text:
# Experimental verification of a fully inseparable tripartite continuous-variable state

## Abstract

A continuous-variable tripartite entangled state is experimentally generated by combining three independent squeezed vacuum states and the variances of its relative positions and total momentum are measured. We show that the measured values violate the separability criteria based on the sum of these quantities and prove the full inseparability of the generated state.

The remarkable proposal of quantum teleportation Bennett et al. (1993) demonstrates that the quantum correlations of a shared entangled state enable two parties to reliably exchange quantum information. So far, several experiments on quantum communication with discrete-variable states have been carried out. In the domain of continuous variables (CVs), the unconditional quantum teleportation of arbitrary coherent states Furusawa et al. (1998); Bowen et al. (2003); Zhang et al. (2003) and quantum dense coding Li et al. (2002) have been demonstrated. These successful experiments show the advantage of CV bipartite entanglement for the implementation of quantum protocols; that is, the simplicity of its generation and manipulation and the applicability of efficient homodyne techniques to its detection.

CV entanglement may also be applicable to quantum protocols involving more than two parties. For example, tripartite entanglement (the entanglement shared by three parties) enables one to construct a quantum teleportation network van Loock and Braunstein (2000), to build an optimal one to two telecloner van Loock and Braunstein (2001a), or to perform controlled dense coding Zhang et al. (2002). CV tripartite entanglement can be generated in a similar way as in the case of CV bipartite entanglement. It only requires combining three modes using linear optics, where at least one of these modes is in a squeezed state van Loock and Braunstein (2000). In fact, as pointed out in Ref. Giedke et al. (2001), CV tripartite entanglement has already been generated in the CV quantum teleportation experiment of Ref. Furusawa et al. (1998), although no further investigation was made there. On the other hand, the separability properties of tripartite states are more complicated than in the bipartite case; three-mode Gaussian states are classified into five different classes Giedke et al. (2001). In order to exploit the tripartite entanglement for three-party quantum protocols such as that from Ref. van Loock and Braunstein (2000), the state involved has to be fully inseparable (class 1 in Ref. Giedke et al. (2001)). Although the output state that emerges from the beam splitters with one or more squeezed input states is in principle fully inseparable for any nonzero squeezing van Loock and Braunstein (2000), inevitable losses in the real experiment may destroy the genuine tripartite entanglement and convert the state into a partially or fully separable one. This would make a true tripartite quantum protocol fail. In other words, the success of a true tripartite quantum protocol (e.g., a coherent-state quantum teleportation network with fidelities better than one half) is a sufficient criterion for the full inseparability of the state involved van Loock and Braunstein (2000). It should be noted here that the success of a tripartite quantum protocol between two parties with the help of the third party (e.g., via a momentum detection of the mode three) does not guarantee that the third party is inseparable from the rest. In the example of the protocol of Ref. van Loock and Braunstein (2000), the full inseparability can be proven only when the protocol succeeds between at least two different pairs or, more generally, when the positions and momenta of all three parties are part of the protocol.

Though the full inseparability is unambiguously verified this way, an alternative verification scheme that does not rely on a full quantum protocol is desirable. In the bipartite case, the inseparability may also be verified simply by measuring the variances of relative position and total momentum Duan et al. (2000); Simon (2000). Recently, a similar scheme to verify the full inseparability of CV tripartite entangled states was proposed van Loock and Furusawa (2002), based on the variances of appropriate linear combinations in position and momentum. In this letter, we generate a tripartite entangled state by combining three independent squeezed vacuum states and demonstrate its full inseparability by applying the scheme of Ref. van Loock and Furusawa (2002).

Let us introduce the position and momentum quadrature-phase amplitude operators $`\widehat{x}`$ and $`\widehat{p}`$ corresponding to the real and imaginary part of an electromagnetic field mode’s annihilation operator, respectively: $`\widehat{a}=\widehat{x}+i\widehat{p}`$ (units-free with $`\mathrm{}=\frac{1}{2}`$, $`[\widehat{x},\widehat{p}]=\frac{i}{2}`$). The simplest way to generate a tripartite entangled state is to send a single-mode squeezed vacuum state $`|x=0`$ (idealized by an eigenstate corresponding to infinite squeezing) into a series of two beam splitters van Loock and Braunstein (2000). In this case, the inputs of the two unused ports are vacuum states. This is practically easy to implement, but when applied to a quantum protocol, the performance would be of only limited quality due to the two vacuum input states. For example, in the teleportation network, the maximum fidelity between any pair is then $`1/\sqrt{2}`$ in the limit of infinite squeezing (excluding additional local squeezers). In order to approach unit fidelity (perfect teleportation), one needs to send squeezed states into all input ports. An example is the CV counterpart of the Greenberger-Horne-Zeilinger (GHZ) state Greenberger et al. (1990), $`𝑑x|x_1|x_2|x_3`$. This CV GHZ state can be generated by sending a momentum-squeezed vacuum state $`|p=0_1`$ and two position-squeezed vacuum states $`|x=0_2`$ and $`|x=0_3`$ into a “tritter” Braunstein (1998), which consists of two beam splitters with transmittance/reflectivity of 1/2 and 1/1. In order to show this, we define a beam splitter operator $`\widehat{B}_{ij}(\theta )`$ which transforms two input modes $`\widehat{a}_{i,j}`$ as

$`\widehat{B}_{ij}^{}(\theta )\left(\begin{array}{c}\widehat{a}_i\\ \widehat{a}_j\end{array}\right)\widehat{B}_{ij}(\theta )=\left(\begin{array}{c}\widehat{a}_i\mathrm{cos}\theta +\widehat{a}_j\mathrm{sin}\theta \\ \widehat{a}_i\mathrm{sin}\theta \widehat{a}_j\mathrm{cos}\theta \end{array}\right).`$ (5)
The transmittance $`T`$ and the reflectivity $`R`$ of the beam splitter are expressed by $`T=\mathrm{cos}^2\theta `$ and $`R=\mathrm{sin}^2\theta `$, respectively. Applying first $`\widehat{B}_{12}(\mathrm{cos}^11/\sqrt{3})`$ and then $`\widehat{B}_{23}(\pi /4)`$ to the input state $`|p=0_1|x=0_2|x=0_3`$ yields $`𝑑x|x_1|x_2|x_3`$. This CV GHZ state is a simultaneous eigenstate of zero total momentum ($`p_1+p_2+p_3=0`$) and zero relative positions ($`x_ix_j=0`$) and exhibits maximal entanglement.

In the real experiment, only finite squeezing is available. Thus the output state is no longer the ideal CV GHZ state and it can never be maximally entangled. Accordingly, total momentum and relative positions have finite variances: $`[\mathrm{\Delta }(\widehat{p}_1+\widehat{p}_2+\widehat{p}_3)]^2>0`$ and $`[\mathrm{\Delta }(\widehat{x}_i\widehat{x}_j)]^2>0`$. This becomes clear when we express the operators for the three output modes in the Heisenberg picture van Loock and Braunstein (2000):

$`\widehat{x}_1`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{3}}}e^{+r_1}\widehat{x}_1^{(0)}+\sqrt{{\displaystyle \frac{2}{3}}}e^{r_2}\widehat{x}_2^{(0)},`$
$`\widehat{p}_1`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{3}}}e^{r_1}\widehat{p}_1^{(0)}+\sqrt{{\displaystyle \frac{2}{3}}}e^{+r_2}\widehat{p}_2^{(0)},`$
$`\widehat{x}_2`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{3}}}e^{+r_1}\widehat{x}_1^{(0)}{\displaystyle \frac{1}{\sqrt{6}}}e^{r_2}\widehat{x}_2^{(0)}+{\displaystyle \frac{1}{\sqrt{2}}}e^{r_3}\widehat{x}_3^{(0)},`$
$`\widehat{p}_2`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{3}}}e^{r_1}\widehat{p}_1^{(0)}{\displaystyle \frac{1}{\sqrt{6}}}e^{+r_2}\widehat{p}_2^{(0)}+{\displaystyle \frac{1}{\sqrt{2}}}e^{+r_3}\widehat{p}_3^{(0)},`$
$`\widehat{x}_3`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{3}}}e^{+r_1}\widehat{x}_1^{(0)}{\displaystyle \frac{1}{\sqrt{6}}}e^{r_2}\widehat{x}_2^{(0)}{\displaystyle \frac{1}{\sqrt{2}}}e^{r_3}\widehat{x}_3^{(0)},`$
$`\widehat{p}_3`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{3}}}e^{r_1}\widehat{p}_1^{(0)}{\displaystyle \frac{1}{\sqrt{6}}}e^{+r_2}\widehat{p}_2^{(0)}{\displaystyle \frac{1}{\sqrt{2}}}e^{+r_3}\widehat{p}_3^{(0)}.`$ (6)
Here a superscript $`(0)`$ denotes initial vacuum modes, and $`r_1`$, $`r_2`$, and $`r_3`$ are the squeezing parameters. In addition to the finite squeezing, the inevitable losses in the experiment further degrade the entanglement. It is important to stabilize the relative phase of the three input modes in order to properly adjust the squeezing directions. The phase fluctuations in this stabilization lead to an extra degradation of the entanglement. As a result, the output state does not necessarily exhibit genuine tripartite entanglement: it may be fully or partially separable. Therefore, we need to experimentally verify the full inseparability of the state.

A feasible scheme for this purpose is to check the following set of inequalities van Loock and Furusawa (2002):

$`\mathrm{I}.[\mathrm{\Delta }(\widehat{x}_1\widehat{x}_2)]^2+[\mathrm{\Delta }(\widehat{p}_1+\widehat{p}_2+g_3\widehat{p}_3)]^21,`$
$`\mathrm{II}.[\mathrm{\Delta }(\widehat{x}_2\widehat{x}_3)]^2+[\mathrm{\Delta }(g_1\widehat{p}_1+\widehat{p}_2+\widehat{p}_3)]^21,`$
$`\mathrm{III}.[\mathrm{\Delta }(\widehat{x}_3\widehat{x}_1)]^2+[\mathrm{\Delta }(\widehat{p}_1+g_2\widehat{p}_2+\widehat{p}_3)]^21.`$ (7)
Here, the $`g_i`$ are arbitrary real parameters. Note that the variances of the vacuum state are $`(\mathrm{\Delta }\widehat{x}_i^{(0)})^2=(\mathrm{\Delta }\widehat{p}_i^{(0)})^2=\frac{1}{4}`$. The violation of inequality $`\mathrm{I}.`$ is a sufficient condition for the inseparability of modes 1 and 2 and is a criterion for the success of a quantum protocol between parties 1 and 2. Note that inequality $`\mathrm{I}.`$ alone does not impose any restriction on the separability of mode 3 from the others. In other words, the success of a quantum protocol between parties 1 and 2 with the help of party 3 (by conveying classical information about a measurement of $`\widehat{p}_3`$ van Loock and Braunstein (2000)) does not prove the inseparability of the third party from the rest. Thus, we need to check the violation of at least two of the three inequalities (7) to verify the full inseparability of the tripartite entangled state van Loock and Furusawa (2002).

From Eq. (6) we find that the optimum gain $`g_i^{\mathrm{opt}}`$ to minimize the l.h.s. of Ineq. (7) depends on the squeezing parameters, namely

$`g_i^{\mathrm{opt}}={\displaystyle \frac{e^{+2r_2}e^{2r_1}}{e^{+2r_2}+\frac{1}{2}e^{2r_1}}},`$ (8)
where $`r_2=r_3`$ (which makes the three-mode state totally symmetric and hence $`g_i^{\mathrm{opt}}`$ independent of $`i`$). In the case of infinite squeezing (CV GHZ state), the optimum gain $`g_i^{\mathrm{opt}}`$ is one, while it is less than one for finite squeezing. Although the smallest values of the l.h.s. of Ineq. (7) are observed when we experimentally adjust $`g_i^{\mathrm{opt}}`$, we employ $`g_i=1`$ for all $`i`$. This makes the experimental verification simpler. Moreover, the measured variances then directly correspond to those of the eigenvalues of the ideal CV GHZ state (relative positions and total momentum).

Figure 1 shows the schematic of the experimental setup to generate three independent squeezed vacuum states. We use a subthreshold degenerate optical parametric oscillator (OPO) with a potassium niobate crystal (length 10mm). Each OPO cavity is a bow-tie type ring cavity which consists of two spherical mirrors (radius of curvature 50mm) and two flat mirrors. The round trip length is 500mm and the waist size in the crystal is 20$`\mu `$m. An output of a Ti:Sapphire laser at 860nm is frequency-doubled in an external cavity with the same configuration as for the OPOs and divided into three beams to pump three OPOs. The pump powers are 56, 71, and 78 mW for OPO 1, 2, and 3, respectively. The squeezed vacuum outputs from these OPOs are combined at two beam splitters to generate the approximate CV GHZ state (see Fig. 3). The visibilities of this combination are 0.968 for the input modes 1 and 2, and 0.948 for 2 and 3. The output modes from the beam splitters are fed into the homodyne detectors 1, 2, and 3 with local oscillator (LO) powers of 1.3, 1.7, and 1.5 mW, and visibilities between the input modes to the homodyne detectors and LOs of 0.979, 0.971, and 0.989, respectively.

We first measure the noise power of each output mode. Figure 2 shows the measurement results on output mode 1. The minimum noise level of 1.14$`\pm `$0.25 dB compared to the corresponding vacuum noise level is observed for the $`x`$ quadrature, while the maximum noise level of 4.69$`\pm `$0.26 dB is observed for the $`p`$ quadrature. Similarly, the minimum noise levels of 0.75$`\pm `$0.27 and 1.21$`\pm `$0.29 dB for the $`x`$ quadrature and the maximum noise levels of 4.12$`\pm `$0.27 and 4.69$`\pm `$0.21 dB for $`p`$ are observed for output modes 2 and 3, respectively. Note that the observed noise levels are always above the corresponding vacuum noise level.

Next we measure the variances of the relative positions and the total momentum from Ineq. (7). Figure 3(a) shows the schematic of the measurement of the variances $`[\mathrm{\Delta }(\widehat{x}_3\widehat{x}_1)]^2`$ and $`[\mathrm{\Delta }(\widehat{x}_2\widehat{x}_3)]^2`$. The outputs of the homodyne detection are electronically subtracted, and the noise power is measured by spectrum analyzers. The variance $`[\mathrm{\Delta }(\widehat{x}_1\widehat{x}_2)]^2`$ is measured in a similar manner. In the case of the variance $`[\mathrm{\Delta }(\widehat{p}_1+\widehat{p}_2+\widehat{p}_3)]^2`$, the noise power of the electronical sum of the homodyne detection outputs is measured as shown in Fig. 3(b).

Figure 4 shows a series of measurement results of (a) $`[\mathrm{\Delta }(\widehat{x}_1\widehat{x}_2)]^2`$, (b) $`[\mathrm{\Delta }(\widehat{x}_2\widehat{x}_3)]^2`$, (c) $`[\mathrm{\Delta }(\widehat{x}_3\widehat{x}_1)]^2`$, and (d) $`[\mathrm{\Delta }(\widehat{p}_1+\widehat{p}_2+\widehat{p}_3)]^2`$, which have the average noise power of $`1.95`$, $`2.04`$, $`1.78`$, and $`1.75`$ dB, respectively, compared to the corresponding vacuum noise level. These results clearly show the nonclassical correlations among the three modes. After repeating the measurement series ten times, we obtain the following measured values for the l.h.s. of Ineq. (7),

$`\mathrm{I}.[\mathrm{\Delta }(\widehat{x}_1\widehat{x}_2)]^2+[\mathrm{\Delta }(\widehat{p}_1+\widehat{p}_2+\widehat{p}_3)]^2`$
$`=0.851\pm 0.062<1,`$
$`\mathrm{II}.[\mathrm{\Delta }(\widehat{x}_2\widehat{x}_3)]^2+[\mathrm{\Delta }(\widehat{p}_1+\widehat{p}_2+\widehat{p}_3)]^2`$
$`=0.840\pm 0.065<1,`$
$`\mathrm{III}.[\mathrm{\Delta }(\widehat{x}_3\widehat{x}_1)]^2+[\mathrm{\Delta }(\widehat{p}_1+\widehat{p}_2+\widehat{p}_3)]^2`$
$`=0.867\pm 0.062<1.`$ (9)
Since violations of all the inequalities are demonstrated, we have proven the full inseparability of the generated tripartite entangled state.

In summary, we have generated a tripartite CV entangled state and verified its full inseparability according to the criteria based on the variances of the relative positions and the total momentum. The violations of all the inequalities verify the presence of genuine tripartite entanglement. Moreover, they ensure that a suitable true tripartite quantum communication protocol using the generated state would succeed. For example, a fidelity greater than one half would be achievable between any pair of parties in a tripartite quantum teleportation network with arbitrary coherent signal states. Let us finally note that the generated tripartite entangled state, though being a Gaussian state with an always positive Wigner function, potentially exhibits nonlocal correlations between the three modes: it can be shown van Loock and Braunstein (2001b) that this state violates three-party Bell-Mermin-Klyshko inequalities imposed by local realism Bell (1964); Mermin (1990); Klyshko (1993). However, instead of the relatively simple homodyne measurements, photon number parity measurements are required to detect this nonlocality Banaszek and Wódkiewicz (1998) which is unfeasible with current technology.

This work was supported by the MEXT and the MPHPT of Japan, and Research Foundation for Opto-Science and Technology. PvL acknowledges financial support from the DFG under the Emmy-Noether programme.


format: text

text:
# Recursion method and one-hole spectral function of the Majumdar-Ghosh model

## 1 Introduction

The one-hole spectral function (or the imaginary part of the one-particle Green’s function) is a central quantity in the theory of strongly correlated systems since it contains the information about the charge dynamics in a correlated background Fulde . It was often studied in the theory of high-temperature superconductors (HTSC) using the several variants of the $`t`$ -$`J`$ or the Hubbard model Dagotto . Other examples are the Kondo peak seen in the spectral function of heavy fermion compounds Hewson , or the features of spin-charge separation for one-dimensional problems Voit . The most important experimental method to measure the one-hole spectral function is angle-resolved photoemission.

A powerful method to calculate the one-hole spectral function for strongly correlated systems consists in the exact diagonalisation of small clusters Dagotto  using the Lanczos method Lanczos . That can be understood like a tridiagonalisation of the Hamilton matrix. The method and its application in electronic structure theory has already a long history. It is also known as recursion method according to pioneering studies of Haydock and coworkers in the seventies Haydock72 ; Haydock80 . That time, however, it was mostly applied to different problems of noninteracting electrons. The numerous spectral function calculations using the different versions of the $`t`$-$`J`$ or the Hubbard model Dagotto  have usually very large finite size effects. Normally, they provide the spectral function only as a sequence of well separated $`\delta `$-peaks. On the other hand, it is already known for long time that the recursion method has the potential to provide the information about the bulk spectral function by a proper termination of the continued fraction representation of the Green’s function HayNex85 ; Ducastelle . That idea was used in the field of many-body dynamics for pure spin systems by Viswanath and Müller VisMul .

Below we would like to present an improvement of the recursion method for the charge dynamics in strongly correlated systems. By finding the correct asymptotic behaviour of the recursion coefficients $`a_n`$ and $`b_n`$ for $`n\mathrm{}`$ which is not influenced by boundary effects we are able to extract the bulk spectral function. The method will be demonstrated for the example of one-hole in the Majumdar-Ghosh (MG) model, i.e. in the 1D spin chain with exchange terms between first $`J`$ and second neighbours $`J^{}`$ at the special ratio $`J^{}/J=1/2`$. It has the advantage that one knows explicitly the ground-state wave function of the spin system Majumdar . The MG model is considered here as a generic model for a 1D spin chain with a spin gap. That means that the spin correlation functions decay exponentially. Therefore, the MG model is a good candidate to present our method since the fast decay of spin correlation functions tends to minimise the boundary effects on the recursion coefficients. We analyse also the recursion coefficients corresponding to a recently presented variational study of the one-hole spectral function using the subspace of one-spinon wave-functions Hayn . It completes our approach by providing additional information on the asymptotic behaviour of recursion coefficients and on the spectral function. The main physical question connected with the MG model is: how are the features of spin-charge separation in 1D modified by a spin gap. We would like to present a rather complete description of the one-hole spectral function of the MG model. The paper is organised as follows: first we recall the recursion method (Sec. 2) and introduce the one-dimensional $`tJJ^{}`$ model (Sec. 3). In Sec. 4 we collect analytical and numerical results in the strong coupling limit $`J,J^{}0`$ and in Sec. V we present the asymptotic behaviour of the recursion coefficients. That analysis is used to reconstruct the spectral function (Sec. 6) and to calculate the bound states and band gaps (Sec. 7).

## 2 The Recursion Method

The recursion method (RM) Haydock72 ; Haydock80  is a powerful tool to calculate the matrix elements of the resolvent ($`z\widehat{H})^1`$ of a Hamiltonian $`\widehat{H}`$. Given an arbitrary starting state $`|u_0`$, the RM generates a new basis in which the Hamiltonian matrix is tridiagonal. One can formulate the RM in two ways, in the Hamiltonian or the Liouvillian representation VisMul . The latter one corresponds to a generalisation of the RM to the Liouvillian space of operators. Usually, it is used for finite temperature problems. In the present case, however, both formulations will be demonstrated to be identical and we concentrate on the Hamiltonian variant.

In the Hamiltonian representation, the RM starts from a vector in Hilbert space, i.e. from some wave function of the system $`|u_0`$. The new basis (orthogonal, but not normalised) is generated according to the Lanczos procedure

$$|u_{n+1}=(\widehat{H}a_n)|u_nb_n^2|u_{n1}$$
(1)
with $`|u_1=0`$, and $`b_0^2=u_0|u_0`$. The coefficients are calculated from

$$a_n=u_n|\widehat{H}|u_n/u_n|u_n,b_n^2=u_{n+1}|u_{n+1}/u_n|u_n.$$
(2)
Having calculated the recursion coefficients $`a_n`$ and $`b_n`$, one can easily find the matrix element of the resolvent with $`|u_0`$:

$`R(z)`$ $`=`$ $`u_0|(z\widehat{H})^1|u_0`$
$`=`$ $`{\displaystyle \frac{b_0^2}{za_0\frac{b_1^2}{za_1\frac{b_2^2}{\mathrm{}}}}}{\displaystyle \frac{b_0^2}{za_0}}{\displaystyle \frac{b_1^2}{za_1}}\mathrm{}.`$
For a starting state corresponding to one hole in the correlated state, the resolvent $`R(z)`$ is identical to the one-particle Green’s function and its imaginary part gives the spectral function we are interested in.

In this form the RM was applied to many body problems (see e.g. Refs. GaglBacc ; Dagotto ; mb  and references therein) several times. In these cases the resolvent matrix elements of the full many body Hamiltonian were calculated. This work can be done numerically for finite systems. Unfortunately, at present, such calculations are only possible for very small clusters. The resulting $`R(z)`$ may be related with the dynamic susceptibility GaglBacc . The problem is that $`R(z)`$ for a finite system is a rational function that may be represented as a set of poles. The corresponding spectral density consists of a set of $`\delta `$-functions which are usually represented using an artificial broadening. This is quite different from that what we can expect for a macroscopic system of interacting particles. The spectrum of a nontrivial interacting system contains usually a continuous part (the “continuum” which is connected with a branch cut in $`R(z)`$) and eventually one or several isolated states (the point spectrum)). A sequence of broadened $`\delta `$-functions is a very poor approximation for the continuous part of the spectrum. Furthermore, such a procedure neglects completely the important difference between the continuum and isolated poles. All the fine structure of the spectrum is completely lost.

Below we will show what kind of information about the macroscopic system may be extracted from the exact diagonalisation studies of small clusters and propose a more physical procedure of smoothing the spectral density. The ED has been performed for the 1D $`tJJ^{}`$ model with 12, 16, 20, and 24 lattice sites according to the method described in Refs. GaglBacc ; Dagotto . From these calculations we extract the information about the bulk related asymptotic behaviour of the recursion coefficients $`a_n`$ and $`b_n`$ for large $`n`$. Additional information about the asymptotics can be obtained using the variational approach presented in Ref. Hayn .

## 3 Majumdar-Ghosh model

In order to demonstrate our method, we specialise to the calculation of the one-hole spectral function in the one-dimensional Majumdar-Ghosh model (the 1D $`tJJ^{}`$ model with $`J^{}/J=1/2`$). The interest to 1D models has been revived after the creation of quasi one dimensional transition metal compounds. The reports on angle-resolved photoemission spectroscopy (ARPES) K96  studies show that the most striking theoretical prediction for 1D strongly correlated systems, the spin-charge separation, can be observed. The interest to this phenomenon has increased after Anderson’s proposal that it explains the properties of 2D cuprate superconductors Anderson . The calculation of one-particle Green’s functions for 1D strongly correlated system is a complex and non-trivial task. The complete answer is absent even for the exactly solvable 1D Hubbard model LiWu . That is why the development of simple approximate approaches remains important. In the present study we consider the $`tJJ^{}`$ model, which describes the basic physics of the electron motion in 1D transition metal compounds

$$\widehat{H}=\widehat{t}+\widehat{J}+\widehat{J}^{},$$
(4)
where

$$\widehat{t}=t\underset{i,\alpha }{}(X_i^{\alpha 0}X_{i+1}^{0\alpha }+h.c.),$$

$$\widehat{J}=\frac{J}{2}\underset{i,\alpha ,\beta }{}X_i^{\alpha \beta }X_{i+1}^{\beta \alpha },\widehat{J}^{}=\frac{J^{}}{2}\underset{i,\alpha ,\beta }{}X_i^{\alpha \beta }X_{i+2}^{\beta \alpha },$$

and we introduce Hubbard projection operators that act in the subspace of on-site states, namely boson-like operators

$$X_i^{\alpha \beta }|\alpha ,i\beta ,i|,\alpha ,\beta =,,$$

and fermion-like ones $`X_i^{\alpha 0}`$, $`i`$ being the site index. To make contact with the standard notation we note that

$$\frac{1}{2}\underset{\alpha ,\beta }{}X_i^{\alpha \beta }X_j^{\beta \alpha }=𝐒_i𝐒_j+\frac{1}{4}n_in_j,$$

for $`ij`$, where $`𝐒_i`$ and $`n_i`$ are spin and density operators, respectively. The exchange term to second neighbours $`J^{}`$ is especially relevant in 1D compounds with a spin-Peierls transition, like CuGeO3. For $`J^{}/J`$ larger or equal 0.2411 the quasi long-range order (algebraically decaying spin-correlation functions) of the 1D antiferromagnetic spin-half Heisenberg model is destroyed and one observes a gap in the spin-excitation spectrum CrtFrst .

Our aim is to calculate the one-particle two-time retarded Green’s function

$$G(k,z)=X_k^{\sigma 0}|X_k^{0\sigma }$$
(5)
for $`z=\omega +ı0^+`$ and the spectral density

$$A(k,\omega )=\frac{1}{\pi }\mathrm{Im}G(k,\omega +ı0^+),$$
(6)
where $`X_k^{\sigma 0}=\sqrt{2/L}_m\mathrm{e}^{ıkm}X_m^{\sigma 0}`$, and $`L`$ is the number of sites, the factor $`\sqrt{2}`$ is introduced for the purpose of normalisation. The above one-particle Green’s function describes the motion of one hole in the correlated state described by the Majumdar-Ghosh wave function Majumdar  (from now on we consider the special case $`J^{}/J=1/2`$). The notation means

$$A|Bi_t^{}^{\mathrm{}}𝑑t\mathrm{e}^{i\omega (tt^{})}\{A(t),B(t^{})\},$$
(7)
with $`\{\mathrm{},\}`$ denoting the anticommutator, and where the expectation value means the thermal average over a grand canonical ensemble: $`\mathrm{}=Q^1Sp[\mathrm{e}^{\beta (H\mu N)}\mathrm{}],Q=Sp\mathrm{e}^{\beta (H\mu N)}`$. Here $`Sp`$ implies the trace of an operator, $`N`$ is the particle number operator, $`\beta =(kT)^1`$ is an inverse temperature, and $`\mu `$ represents the chemical potential. The time dependence of the operator $`A(t)`$ is given by $`A(t)=\mathrm{e}^{it(H\mu N)}A\mathrm{e}^{it(H\mu N)}`$. At zero temperature, $`\mathrm{}`$ goes over into the expectation value with the ground-state wave function $`|\mathrm{\Psi }_0`$. It is not difficult to see that in the given case the one-particle Green’s function (5) may be formulated as a resolvent matrix element

$$G(k,z)=X_k^{\sigma 0}|(z+E_0\widehat{H})^1|X_k^{\sigma 0},$$
(8)
with the state $`|X_k^{\sigma 0}=X_k^{\sigma 0}|\mathrm{\Psi }_0`$ denoted in the same way as the corresponding operator. It means also that the Hamiltonian and Liouvillian formulation of the RM are now equivalent if we take into account the shift of the diagonal recursion coefficients $`a_n`$ by the ground-state energy $`E_0`$ of the pure spin system.

When the ground state is degenerate, as it is the case in the Majumdar-Ghosh model, we have to take $`Sp`$ over the ground-state manifold in the calculation of $`\mathrm{}`$ in (7). We have found that the spectral function calculated for two orthogonal ground states of the model according to the prescription given in Sec. 6 coincides within the accuracy of the method.

## 4 Strong coupling limit

For the following analysis it is important to recall some former results for the one-hole spectral function in the MG model using a variational ansatz Hayn . It was shown that one can give an exact result in the strong coupling limit $`J,J^{}0`$. In this limit only the $`\widehat{t}`$ operator remains in Eq. (4). Note that it is a true many-body Hamiltonian due to the constraint of no double occupancy, it is often called the $`\widehat{t}`$-model. Formally, the $`\widehat{t}`$-model is the limiting case of various models ($`U\mathrm{}`$ of the Hubbard model, $`J0`$ of the $`tJ`$ model, etc.), and its ground state is strongly degenerate. The degeneration is removed by an infinitesimal perturbation that fixes the ground state. It is clear that different spin-models give different answers for the Green’s function (5) and the spectral density (6) in the limit of the $`\widehat{t}`$-model since they may differ by the spin ground state Hayn .

The exact result was found using the basis operator set

$$v_{m,r}=\sqrt{\frac{2}{L}}\underset{\alpha _1,\mathrm{},\alpha _r}{}X_m^{\sigma \alpha _1}X_{m+g}^{\alpha _1\alpha _2}\mathrm{}X_{m+rg}^{\alpha _{r1}\alpha _r}X_{m+r}^{\alpha _r0},$$
(9)
$`g=r/\left|r\right|`$ . The operator $`v_{m,r}`$ may be interpreted as holon (’right end’ $`\mathrm{}X_{m+r}^{\alpha _r0}`$) and spinon (’left end’ $`X_m^{\sigma \alpha _1}\mathrm{}`$) excitations connected by a string of spin ’flips’ (in fact, this string contains also diagonal spin operators but they also create excitations in a quantum antiferromagnet). The double Fourier transform over the $`m`$ and $`r`$ indices gives exact eigenstates of the $`\widehat{t}`$-term for any magnetic background. In the case of finite exchange values, the basis set (9) is not complete and corresponds to the single spinon approximation.

It was shown that the Green’s function (5) for fixed $`k`$ has a formal analogy with the on-site Green’s function of the 1D tight-binding model of non-interacting electrons in a non-orthogonal basis

$`G(k,z)={\displaystyle \frac{1}{L}}{\displaystyle \underset{Q}{}}{\displaystyle \frac{2Z(Q)}{z+ϵ_h(kQ)}}`$ (10)
$`={\displaystyle \frac{1}{\sqrt{z^24t^2}}}\left[1+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}2\mathrm{cos}kn\mathrm{\Omega }_{0n}\tau ^n(z)\right],`$
where $`ϵ_h(k)2t\mathrm{cos}(k)`$.

The spectral function may be represented like

$$A(k,\omega )=_\pi ^\pi \frac{dQ}{\pi }Z(Q)\delta (\omega +ϵ_h(kQ)),$$
(11)
with

$`Z(Q){\displaystyle \frac{1}{2}}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}\mathrm{exp}\left[ı(Q\pi )n\right]\mathrm{\Omega }_{0n},`$ (12)
$`\mathrm{\Omega }_{0r}={\displaystyle \underset{\alpha _1,\mathrm{},\alpha _r,\sigma }{}}X_0^{\sigma \alpha _1}\mathrm{}X_{rg}^{\alpha _{r1}\alpha _r}X_r^{\alpha _r\sigma }`$
$`=(2𝐒_0𝐒_g+{\displaystyle \frac{1}{2}})(2𝐒_g𝐒_{2g}+{\displaystyle \frac{1}{2}})\mathrm{}(2𝐒_{rg}𝐒_r+{\displaystyle \frac{1}{2}}),`$
$`\tau (z){\displaystyle \frac{z\sqrt{z^24t^2}}{2t}}={\displaystyle \frac{t}{z}}{\displaystyle \frac{t^2}{z}}{\displaystyle \frac{t^2}{z}}\mathrm{}.`$
Integrating over the $`\delta `$-function, the spectral function (11) may be given in the form SP

$$A(k,\omega )=\frac{Z(k+\varphi )+Z(k\varphi )}{\pi \sqrt{4t^2\omega ^2}},\mathrm{cos}(\varphi )=\frac{\omega }{2t}.$$

The expressions (11,4) imply that the spectrum is continuous and bounded on the interval $`2t\omega 2t`$ for all the possible spin models in the considered limit. It may be interpreted as holon motion with the dispersion $`ϵ_h(k)`$ with an immovable spinon. The difference in hole spectral densities and Green’s functions arises due to the expectation value $`Z(Q)`$ that depends on the spin part of the Hamiltonian. Some parts of the spectrum may be excluded by a vanishing $`Z(Q)`$.

For the 1D pure (non frustrated) antiferromagnetic Heisenberg model, $`Z(Q)`$ was given in Ref. OgSh ; SP  as $`Z(Q)\mathrm{\Theta }(Q\pi /2)/\sqrt{Q\pi /2}`$. The square-root singularity leads to additional peaks in the spectral function, such that for each $`k`$ value one can distinguish well pronounced spinon and holon peaks (for details see also Ref. Nagaosa97 ). For the Majumdar-Ghosh (MG) wave function, on the other hand, one finds $`Z(Q)=\frac{3}{2}\left(1+\mathrm{cos}Q\right)/\left(5+4\mathrm{cos}2Q\right)`$ (see Ref. Maekawa ) where $`Z(Q)`$ is non zero for all values of $`Q`$. The peak at $`Q=\pi /2`$ becomes very wide. That leads to a large damping of the holon peak Hayn . One finds the following explicit formula for the spectral function of the MG model in the limit $`J,J^{}0`$:

$`A(k,\omega )={\displaystyle \frac{3}{8t\pi \sqrt{1x^2}}}`$
$`{\displaystyle \frac{20+4x\mathrm{cos}k+16(2x^21)\mathrm{cos}2k8D}{25+40(2x^21)\mathrm{cos}2k+8\mathrm{cos}4k+8(8x^48x^2+1)}},`$
$`x{\displaystyle \frac{\omega }{2t}},Dx\mathrm{cos}3k+4x^3\mathrm{cos}k`$ (13)
which is non zero for all $`2t\omega 2t`$.

Now, for the MG model, we compare the continued fraction expansion of the one-particle Green’s function for the macroscopic system (10) with that for a finite ring. The summation of series in (10) gives

$$G(k,z)=\frac{42\left[2\tau (z)+\tau ^3(z)\right]\mathrm{cos}k\tau ^4(z)}{\sqrt{z^24t^2}\left[4+4\mathrm{cos}2k\tau ^2(z)+\tau ^4(z)\right]},$$
(14)
that may be rewritten using the identity $`\tau (z)=t/(zt\tau (z))`$ in the form

$$G(k,z)=\frac{1}{za_0}\frac{b_1^2}{za_1b_2^2t^1\tau (z)\left[1+F(z)\right]},$$
(15)
where

$`a_0={\displaystyle \frac{t}{2}}ϵ,b_1^2=3t^2(1{\displaystyle \frac{3}{4}}ϵ^2),a_1={\displaystyle \frac{9ϵ^36ϵ}{8b_1^2}}t^3,`$
$`b_2^2={\displaystyle \frac{9(87ϵ^2)t^4}{16b_1^2}},ϵ\mathrm{cos}k,`$ (16)
and the function $`F(z)`$:

$`F(z)`$ $`=`$ $`{\displaystyle \frac{ϵ(1ϵ^2)V(z)}{(87ϵ^2)W(z)}},`$
$`V(z)z(18ϵ^427ϵ^2+8)+4ϵ(109ϵ^2)`$
$`3(12ϵ^421ϵ^224)\tau (z),`$
$`W(z)2z^2(43ϵ^2)2zϵ12`$
$`+11ϵ^2\left[2z(43ϵ^2)ϵ\right]\tau (z),`$
vanishes for $`k=0,\pm \frac{\pi }{2},\pm \pi `$. For these $`k`$ values we thus readily obtain $`G(k,z)`$ in the continued fraction form (2) with $`a_n=a_{\mathrm{}}=0,b_n^2=b_{\mathrm{}}^2=t^2,n>2`$ that follows from (4). For other $`k`$ values the coefficients asymptotically tend to the same constant limits $`a_{\mathrm{}},b_{\mathrm{}}`$.

Fig. 1 compares the analytical result for $`a_n,b_n`$ with the values provided by the ED of a ring with $`L=24`$ sites at $`J=0.01`$. One sees that the values of $`a_n`$ with $`n10`$ are not influenced by the boundary. These first values give us bulk-related information. Comparing with other cluster sizes $`L`$, one finds that $`a_n`$ has no finite size effect for $`n(L4)/2`$. For the off-diagonal elements $`b_n`$ one has $`n(L4)/2+1`$. So, we may conclude that ED studies of finite systems are able to provide the information that concerns the macroscopic system and may be used for the analysis of the spectral density. Despite the small number of ’bulk related’ coefficients they definitely reached their asymptotic regime and may be extrapolated to $`n\mathrm{}`$ giving the exact Green’s function and spectral density for infinite chain. The situation is not always so simple, but this example clearly shows a manageable way to extrapolate ED results to the macroscopic system.

Concluding the consideration of the strong coupling limit we recall that it is commonly believed that the main changes in the spectral density at larger finite couplings (i.e. $`J0`$) come from the appearance of the spinon dispersion $`ϵ_s`$, and may be accounted for by the substitution

$$ϵ_h(kQ)ϵ_h(kQ)+ϵ_s(Q)$$
(17)
in Eqs.(10,11) SP ; Nagaosa97 . Below we shall demonstrate that holon and spinon scattering lead to additional changes: the continuous part of the spectrum expands up to infinite energies and for the case of the MG model, additional discrete states appear.

## 5 Termination of the continued fraction

Now we concentrate on results for finite $`J,J^{}`$. In the previous section we have shown that the few first recursion coefficients given by the exact diagonalisation study of a small cluster are (almost) not affected by boundaries. Thus, in order to infer the shape of the spectral density for macroscopic system we should rely on these ’bulk related’ coefficients. To define the bulk related coefficients we observe from Fig. 2a that the coefficients for $`L=12`$ start to deviate from the common line at $`n=4`$, those for $`L=16`$ at $`n=6`$ and for $`L=20`$ at $`n=8`$. The corresponding numbers for $`b_n`$ (Fig. 2b) are 5, 7, and 9. Generalising we find $`a_n:n<(L4)/2`$ and $`b_n:n<(L4)/2+1`$ to be bulk related coefficients. From the inset of Fig. 2a we clearly see that ’bulk related’ coefficients for finite $`J,J^{}`$ show a different tendency (to grow with moderate oscillations) compared to the rest of the sequence for a cluster (strong oscillations around a constant value).

For finite $`J`$ the number of bulk related coefficients is slightly reduced in comparison to Fig. 1. It should be noted that the ’bulk related’ coefficients with large $`n`$ have still a small finite size effect, but that is several orders of magnitude smaller than for the rest of the sequence. We have also observed that in the case of the pure $`tJ`$ model (i.e. $`J^{}=0`$) the situation is worse. The number of recursion coefficients that are not influenced by finite size effects is considerably reduced. One may recall that according to Eq. (2) the coefficients are related to static spin correlation functions of the spin background and speculate that the difference arises due to the algebraically decaying spin correlation functions for the 1D Heisenberg model in difference to the exponential decay for the MG model.

For the calculation of the spectral density we will try to extrapolate the tendency shown by the ’bulk related’ coefficients. In other words, we rewrite the expression (2) in the form

$$G(k,z)=\frac{b_0^2}{za_0}\frac{b_1^2}{za_1}\mathrm{}\frac{b_{n_0}^2}{za_{n_0}T_{n_0}(k,z)},$$
(18)
and try to find a function $`\stackrel{~}{T}_{n_0}`$ (a so called ”terminator”) that is close to $`T_{n_0}`$. The terminator should have such an analytic behaviour which corresponds to the asymptotics of $`a_n`$ and $`b_n`$.

Various ways to construct such approximations are described in the literature on the recursion method (see Refs. HayNex85 ; Ducastelle ; Nex78 ; Magnus85 ). The asymptotic behaviour of continued fraction coefficients is governed by the band structure and the singularities of the spectral density Magnus85 . The problem is well studied for a bounded spectrum: $`\{a_n\}`$ and $`\{b_n\}`$ converge toward limits in the single band case, oscillate endlessly in a predictable way in the multiband case. Damped oscillations are created by isolated singularities. But the growth of coefficients indicates that the spectrum we deal with is unbounded. The infinite growth of diagonal coefficients is observed for the $`tJ_z`$ model that differs from the $`tJ`$ model by omitting spin fluctuations, $`a_n=nJ_z,b_n^2=z_0t^2`$, ($`z_0`$ is the number of nearest neighbours) Izing96 . This growth has a clear physical meaning: a hole moving on the Neel background (the ground state for $`tJ_z`$ model without holes) creates a string of overturned spins; every $`|u_n`$ state (1) is the combination of such strings containing $`n`$ overturned spins, thus its energy $`u_n\left|\widehat{t}+\widehat{J}_z\right|u_nnJ_z`$. The non-diagonal coefficients $`b_n`$ are constant. As a result, the hole is localised in the $`tJ_z`$ model, its spectral function does not depend on quasimomentum $`k`$ and consists of an infinite set of $`\delta `$ -functions. The spectrum is discrete and unbounded IzNote .

When spin fluctuations are added the situation becomes more complicated, now the non-diagonal coefficients also grow. Such a situation was studied in the two-dimensional $`tJ`$ model as well as for spin-fermion models of the CuO2 plane, where the following asymptotics was found within the framework of the self-consistent Born approximation SCBA ; linterm :

$$b_n\lambda _1n+\lambda _2,a_n2\lambda _1n+\lambda _3,\lambda _i=\lambda _i(k),n1.$$
(19)
It indicates the existence of an exponential tail in the spectral density. But the connection between an unbounded spectrum (toward high energies) and the growth of recursion coefficients is far more general. For example, the model spectral density

$$A(\omega )=c^t\omega ^\rho \mathrm{exp}(\omega ^\beta ),c(\beta )=\frac{2\mathrm{\Gamma }(\beta )}{\left[\mathrm{\Gamma }(\beta /2)\right]^2},\omega >0,$$

with the parameters $`t`$, $`\rho `$, and $`\beta `$, gives the following asymptotic behaviour of the recursion coefficients Magnus85 :

$$b_n\left(\frac{n}{c(2\beta )}\right)^{1/\beta },a_n2\left(\frac{n}{c(2\beta )}\right)^{1/\beta },$$

which does not depend on $`t`$ and $`\rho `$. For $`\beta =1`$ we have a linear dependence on $`n`$ for $`\{a_n\}`$ and $`\{b_n\}`$ Haydock80 .

The linear growth holds also for the MG model. Fig. 3 shows the linear contributions to the growth of ’bulk related’ coefficients. We have found $`\lambda _1,\lambda _2`$ using least mean square fit of the $`\{b_n\}`$ sequence and only $`\lambda _3`$ was adjusted for $`\{a_n\}`$. Fig. 3 shows that our assumption about the slope $`da_n/dn2db_n/dn`$ holds.

The continued fraction expansion of the incomplete Gamma function has the same asymptotics (19) and is written as NR

$`\mathrm{\Gamma }(\alpha ,x)={\displaystyle \frac{\mathrm{e}^xx^\alpha }{x+1\alpha }}{\displaystyle \frac{1(1\alpha )}{x+3\alpha }}{\displaystyle \frac{2(2\alpha )}{x+5\alpha }}\mathrm{}`$
$`\mathrm{}{\displaystyle \frac{n(n\alpha )}{x+(2n+1\alpha )\mathrm{}}}.`$ (20)
We shall use this circumstance for the construction of the terminator $`\stackrel{~}{T}_N(k,z)`$ for $`G(k,z)`$ (18).

Let us introduce the model Green’s function

$$\stackrel{~}{G}(k,z)=\frac{\mathrm{\Gamma }(\alpha ,x\mathrm{e}^{ı\pi })\mathrm{e}^{ı\pi \alpha }}{\lambda _1\mathrm{e}^xx^\alpha },$$
(21)
with the spectral density

$$\stackrel{~}{A}(k,\omega )=\frac{1}{\lambda _1\mathrm{\Gamma }(1\alpha )}\left[\frac{(\omega \omega _0)}{\lambda _1}\right]^\alpha \mathrm{exp}\left[\frac{(\omega \omega _0)}{\lambda _1}\right].$$
(22)
Here $`x=(z\omega _0)/\lambda _1`$,

$$\omega _0\lambda _32\lambda _2\lambda _1,$$
(23)
is the lower bound of the continuous spectrum, $`\alpha =2(M\lambda _2/\lambda _1)`$, $`M`$ is the integer part of $`\lambda _2/\lambda _1`$. The continued fraction expansion coefficients for $`\stackrel{~}{G}(k,z)`$ are

$`\stackrel{~}{a}_m=\lambda _1(2m+1\alpha )+\omega _0,m0,\stackrel{~}{b}_0=1,`$ (24)
$`\stackrel{~}{b}_m=\lambda _1\sqrt{m(m\alpha )}\lambda _1(m\alpha /2),m>0.`$
Comparing with the definition (19) we find

$$b_n\stackrel{~}{b}_{M+n},a_n\stackrel{~}{a}_{M+n}.$$

Then we slightly generalise the recipe of Ref. HayNex85  and approximate linterm

$`\stackrel{~}{T}_{M+n_0}(k,z)`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{q}_{M+n_01}(k,z)\stackrel{~}{G}(k,z)\stackrel{~}{p}_{M+n_0}(k,z)}{\stackrel{~}{q}_{M+n_02}(k,z)\stackrel{~}{G}(k,z)\stackrel{~}{p}_{M+n_01}(k,z)}}`$ (25)
$``$ $`T_{n_0}(k,z).`$
Here $`\stackrel{~}{p}_n(k,z),\stackrel{~}{q}_n(k,z)`$ are polynomials of first and second kind orthogonal with respect to the spectral density (22). They are calculated by the recursion relation Haydock80 ; Nex84

$`\stackrel{~}{p}_1=\stackrel{~}{q}_1=0,\stackrel{~}{p}_0=1,\stackrel{~}{q}_0=\stackrel{~}{b}_0^2,`$ (26)
$`\stackrel{~}{p}_n=(z\stackrel{~}{a}_{n1})\stackrel{~}{p}_{n1}\stackrel{~}{b}_{n1}^2\stackrel{~}{p}_{n2},`$
$`\stackrel{~}{q}_{n1}=(z\stackrel{~}{a}_{n1})\stackrel{~}{q}_{n2}\stackrel{~}{b}_{n1}^2\stackrel{~}{q}_{n3}.`$
$`\stackrel{~}{p}_n(k,z)`$ are proportional to the generalised Laguerre polynomials. In our calculations it holds always: $`\lambda _2/\lambda _1>0`$ and $`1<\alpha <0`$.

We thus obtain $`G(k,z)`$ in the whole complex energy plane. It has the correct analytic properties and coincides with the retarded Green’s function for $`z`$ in the upper half plane and with the advanced Green’s function for $`z`$ in the lower half plane. It has a branch cut on the real axis for $`\omega _0<z<\mathrm{}`$ that corresponds to the continuum part of the spectrum. In addition, it may have isolated poles for $`z<\omega _0`$ which represent quasiparticle excitations.

Figs. 2 and 3 show that the recursion coefficients oscillate around their asymptotics. The influence of these oscillations on the spectral density was neglected so far. For a bounded spectrum this is the indication of a multiband spectrum. The problem was well studied in Refs. Ducastelle ; Magnus85 . In the case of one gap, i.e. a spectrum lying in the disjoint intervals $`E_1\omega E_2`$ and $`E_3\omega E_4`$ ($`E_2<E_3`$) any pair of neighbouring coefficients $`(y,x)=(b_n^2,a_n)`$ or $`(b_{n+1}^2,a_n)`$ obeys the law Ducastelle

$$\left(x^2+A_1x+A_2+2y\right)^2=X(A_1x),$$
(27)
with

$`A_1{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}E_i,A_2{\displaystyle \frac{1}{2}}{\displaystyle \underset{i<j}{}}E_iE_j{\displaystyle \frac{A_1^2}{2}},`$
$`X(z){\displaystyle \underset{i}{}}\left(zE_i\right).`$
This relation may be represented in phase space. In the asymptotic limit, all pairs $`(b_n^2,a_n)`$ and $`(b_{n+1}^2,a_n)`$ should lie on a single closed curve (27). The best fit provides the band edges directly. After establishing the $`E_i`$ values, the relation (27) may be used for the extrapolation of coefficients beyond the known values. Fig. 4 shows the plot for the single spinon approximation (which leads to a bounded spectrum) at $`k=0`$ and $`J=0.4`$ . Several known values and their extrapolation form a single curve. An analogous analysis for different $`k`$ and $`J`$ has shown that in the single spinon approximation we have one gap in the continuous part of the spectrum. The asymptotic regime is reached for $`n>n_020`$ for $`k=0,\pi `$. For the more structured spectral density at other $`k`$ values the asymptotics is reached after $`n_070`$. The bounds of the coefficients are

$$\frac{E_1+E_4}{2}g<a_n<\frac{E_1+E_4}{2}+g,g\frac{E_3E_2}{2},$$

$$\frac{E_4E_1}{4}\frac{g}{2}b_n\frac{E_4E_1}{4}+\frac{g}{2},$$

and the half gap value $`g0.074`$ equals the amplitude of oscillations. This relation between the amplitude of the coefficient oscillations and the gap does not depend on the total bandwidth $`E_4E_1`$ . One may expect that it holds also for an unbounded spectrum. We have modelled such a spectrum by the creation of a gap in the spectral density (22). Then we have calculated the coefficients from the gapped spectral density and we have found that they oscillate around their asymptotics (19) with an amplitude equal to half of the gap. As we have mentioned above, the internal singularities lead to damped oscillations of the coefficients. The restricted number of exact ’bulk-related’ ED coefficients makes it impossible at present time to analyse their oscillations in detail and to extract quantitatively the values and positions of the gaps.

## 6 Reconstruction of spectral density

One may summarise the following recipe of spectral density reconstruction from the Lanczos recursion coefficients: i) Establish the number of recursion coefficients that are not strongly affected by the boundaries. ii) Investigate their asymptotics and find a model function that has a continued fraction expansion with similar asymptotics, or equivalently find a reliable procedure to extrapolate the coefficients. iii) Calculate the Green’s function $`G(k,z)`$ using expressions (18,25) (then $`z`$ may be real) or directly (2) (with complex $`z`$, the numerical efforts depend on the distance to the real axis) up to convergence. Then the spectral density results from the definition (6).

As we have already mentioned, the number of bulk-related coefficients for the Majumdar-Ghosh model (4) is $`n_0(L4)/2`$, where $`L`$ is the number of sites. We have performed the ED for rings of the length up to 24 sites, i.e. we have about 10 pairs of coefficients for every $`k`$ value. The last 6 pairs were used for the determination of the parameters of the linear law (19) which are necessary to terminate the continued fraction. Fig. 5 shows the spectral density $`A(k,z),z=\omega +ı\eta `$ for $`k=\pi /2`$, $`J=0.4`$, $`\eta =0.01`$. On panel (a) the spectral density calculated according to the above recipe (solid line) is compared with the spectral density for the $`L=24`$ sites ring (dotted line) calculated with the same broadening $`\eta `$. The qualitative difference of both curves is evident. Panel (b) of the figure shows the result of the single spinon approximation. With a broadening of $`\eta =0.01`$ the overall shape of the curve is similar to the ED result and is characterised by two peaks at $`\omega \pm 2t`$. The small feature near $`\omega 0.4`$ corresponds to a bound state in the gap, but all bound states and gaps are smeared out for a broadening of $`\eta =0.01`$ (compare with Fig. 6b). The essential difference with respect to the ED spectrum (for the single spinon approximation the spectrum is bounded) manifests itself as a difference of the relative weight of the peaks, the higher energy peak of ED being smoothed by the exponential tail of the spectral density that tends up to $`\omega \mathrm{}`$.

Fig. 6 shows the spectral density $`A(k,\omega +ı\eta )`$ for various $`k`$-values, $`J=0.4`$. Panel (a) corresponds to the ED with $`\eta =0.01`$. The interpretation of the various peaks as different collective excitations were already broadly discussed in Ref. Hayn : the lower edge of the continuum corresponds to the spinon excitation with a small dispersion of the order $`J`$. For $`k`$ different from $`\pi /2`$ one can also observe the holon peak with a much larger dispersion (proportional $`t`$) and with a large damping Hayn .

For the single spinon approximation the law (27) makes it possible to reveal the gap in the continuum part of the spectrum at every $`k`$ and the existence of quasiparticle states outside the continuum spectrum and within the gap. To show these features a very small broadening of $`\eta =0.0001`$ was used in panel (b).

## 7 Bound states and band gaps

The existence of a quasiparticle state in the Majumdar-Ghosh model (i.e. a bound state below the continuum) was proved by ED in Ref. Maekawa  by an analysis of the system size dependence of the pole strength $`Z_h`$ corresponding to the lowest eigenvalue. An analogous way was used to show the existence of a quasiparticle in the single spinon approximation Hayn  and in the diagrammatic approach Brenig01 . It was found that the single spinon approximation gives $`Z_h`$ approximately twice as big as that of the ED. That puzzle can be resolved when one takes into account the linear growth of the recursion coefficients which does not show up in the single spinon approximation. The linear growth may have an influence on the weight $`Z_h`$. From a physical point of view the quasiparticle in this model represents the bound state of holon and spinon. Scattering processes with the contribution of high energy multispinon excited states that are not taken into account in the single spinon approximation enlarge the probability of decay of the bound state.

Fig. 7 shows the relative weights of the hole wave function decomposition over recursion vectors $`|u_n(k)`$ for the three lowest eigenenergies (thick solid, thin solid and dashed lines correspondingly) at the band minimum $`k=\pi /2`$ . 90 recursive coefficients were taken: (a) ’bulk-related’ ED coefficients, extrapolated by linear low; (b) single spinon approximation. A qualitatively different behaviour is seen for the spin-polaron quasiparticle wave function (thick line) and the wave functions that belong to the continuum. The different large-$`n`$ behaviour of recursion coefficients leads to different radius of the polaron state. For comparison, the dotted line on panel (b) shows the change in the wave function caused by the addition of the linear growth (19) to the single spinon approximation coefficients. The linear growth leads to a larger radius of the polaron state and to a decrease of the weight $`Z_h`$ by nearly a factor of two.

The physical origin of the gap in the continuum part of the spectrum that is evident in the single spinon approximation may be understood as follows. Let us recall that in the single spinon approximation we diagonalise the Hamiltonian (4) within the subspace spanned by the operator basis (9). (See the details in Sec. V.B and Appendix C of Ref. Hayn .) The Green’s function (5) for every $`k`$ is found from the solution of a generalised eigenvalue problem that formally resembles an effective tight-binding model (that is different for every $`k`$ and is determined by the spin background) with Hamiltonian $`E_{r,r^{}}(k)`$ in a non-orthogonal basis with overlap matrix $`S_{r,r^{}}(k)`$. The index $`r`$ denotes the distance between ’spinon’ and ’holon’ in (9). The spectral density (6) formally coincides with the density of states of the effective model. In the limit $`J,J^{}0`$ both matrices $`E`$ and $`S`$ depend only on $`rr^{}`$ and the effective model is ’translationally invariant’ if we consider $`r`$ and $`r^{}`$ as site indices of the effective model. And we obtain the results (10,11) making the Fourier transformation over $`rr^{}`$. This is the consequence of the degeneracy of all spin states which makes it impossible for the holon to scatter on inhomogeneities of the spin state. By switching on finite $`J`$ and $`J^{}`$ the matrix $`E_{r,r^{}}(k)`$ acquires two contributions. One of them corresponds to the appearance of the spinon coherent motion with the dispersion $`ϵ_s(Q)=2J\mathrm{cos}Q`$. This contribution is also ’translationally invariant’ and leads to changes in the spectral density (11) described by (17). Another contribution comes from the commutation of the ’holon end’ of the operator (9) with the spin part of the Hamiltonian (4). It describes two processes, the loss of magnetic energy due to the presence of the holon, and the holon-spinon scattering. The former participates in the bound state formation. The latter has a contribution that is different for odd and even $`rr^{}`$ , i.e. we observe a ’period doubling’ in the effective Hamiltonian. This is intimately connected with the dimerised nature of the Majumdar-Ghosh wave function. It is interesting that this feature does not lead to a real period doubling in the system, i.e. the Green’s function (5) and the spectral density (6) have the whole Brillouin zone periodicity, as we clearly see in Fig. 6.

## 8 Summary

We propose a new way how to extract the information about the infinite system from the exact diagonalisation of small clusters. It is based on the consideration of Lanczos recursion coefficients that are provided by ED. We have found that the comparison of the results for several cluster sizes allows to obtain the set of the coefficients that are not affected by the finite size effects. These coefficients contain the needed information about the macroscopic system. When it is possible to infer their asymptotic behaviour (that may strongly differ for the ’bulk related part’ compared to the rest of the set for a finite system) we propose to restore the Green’s function and spectral density using the terminator technique HayNex85 . This may considerably improve the shape of the spectral density compared to the direct ED result which gives the spectral function in the form of a set of delta functions.

## Acknowledgements

This work was supported by DFG (436 UKR 113/49/41) and by a NATO Collaborative Linkage Grant (PST.CLG. 976416). R.O.K. thanks for hospitality the IFW Dresden, where the main part of this work has been carried out. For several numerical calculations we used the Cambridge Recursion Library Nex84 .


text:
# SHOCKS AND PARTICLE ACCELERATION IN SUPERNOVA REMNANTS: OBSERVATIONAL FEATURES

## 1. Introduction

Supernova remnant shocks are considered to be the prime source of cosmic rays for energies at least up to the “knee” at $`3\times 10^{15}`$ eV and probably up to the “ankle” at $`10^{18}`$ eV (above which the cosmic rays are thought to be of extra-galactic origin). One of the main reasons is that supernova remnants (SNRs) are the only galactic sources that are able to provide the energy necessary to maintain the observed cosmic ray density, $`10^{48}`$ erg/yr, if about 10% of the average supernova kinetic energy is dumped into cosmic ray acceleration, and assuming a supernova rate of 1 every 50 - 100 yrs (e.g. Longair, 1994). Another reason is that radio synchrotron emission from SNRs have spectral indexes implying a spectral energy slope of $`2.2`$, consistent with the locally observed cosmic ray spectral index, which has a slightly steeper index due to propagation effects.

For a long time the observational evidence that SNRs are indeed the source of cosmic rays consisted solely of the observation of radio synchrotron emission, caused by the presence of relativistic electrons. In some cases SNRs may just light up the relativistic galactic background electrons by means of compressing the local magnetic field (van der Laan, 1962), but for the very bright, young, SNRs, the most likely explanation is that the electrons are accelerated by SNR shocks. A good example is the youngest and brightest known galactic SNR, Cassiopeia A (Cas A).

Although for these reasons it is likely that SNRs accelerate electrons to relativistic energies, at least up to several GeV, radio synchrotron emission does not prove that SNRs are capable of accelerating cosmic rays up to the “knee”, and, moreover, it only proves that electrons are being accelerated, whereas the cosmic rays observed at earth are dominated by ions. However, over the last decade evidence has been accumulating that cosmic rays are indeed being accelerated up to the “knee”. This progress has been made possible by the coming of age of X-ray astronomy, with missions like ASCA, BeppoSAX, RXTE, and, more recently Chandra and XMM-Newton. In addition, a lot has been learned from $`\gamma `$-ray observations with space based instruments like CGRO-EGRET and ground based Cherenkov telescopes such as WhippleHEGRA, and CANGAROO.

Although the detection of ion cosmic rays in SNRs can be considered one of the main goals of cosmic ray physics, another important, but still poorly understood, aspect of cosmic ray physics is the injection mechanism at low energies. This process may be intimately connected to shock heating and temperature equilibration by collisionless shocks.

## 2. Shock structure and electron-ion tempererature equilibration

The particle collision length scales in supernova remnant shocks are much larger than the typical size of the shock structure. For instance, the young supernova remnant Cas A has a shock velocity of about 5000 km/s (Vink et al., 1998; Delaney & Rudnick, 2003) and a typical ion density of 5 cm-3, implying that a proton entering the shock has a relative energy of $`8`$ keV. The deflection time scale for a proton, assuming proton proton interactions, is $`\tau _{\mathrm{pp}}800`$ yr (e.g. Huba, 2002), more than the age of Cas A, which is about $`320`$ yr! So other, collective rather than two body, interactions, have to be at work to heat the gas behind a supernova remnant shocks. They are therefore referred to as collisionless shocks.

Although SNRs have collisionless shocks, the Hugoniot relations, which are based on the conservation of mass, momentum, and energy, are still valid. If we neglect, for the moment, the energy deposition in cosmic rays, and radiative losses, the Hugoniot relations in the limit of high Mach number shocks state (e.g. Zeldovich & Raizer, 1966; McKee & Hollenbach, 1980):

$`{\displaystyle \frac{\rho _{2,i}}{\rho _{1,i}}}=\chi ={\displaystyle \frac{\gamma +1}{\gamma 1}}=4`$ (1)
$`kT_{2,i}={\displaystyle \frac{2(\gamma 1)}{\gamma +1)^2}}m_iv_s^2={\displaystyle \frac{3}{16}}m_iv_s^2,`$ (2)
where subscript $`1`$ and $`2`$ refer to the pre-shock and post-shock conditions for particle species $`i`$; with $`v_s`$ the shock velocity. The numerical values were obtained for $`\gamma =5/3`$, the specific heat ratio for a monatomic gas. Note, however, that for a relativistic gas (i.e. the pressure is dominated by relativistic cosmic rays) $`\gamma =4/3`$, giving a compression ratio of $`\chi =7`$. The post-shock gas moves away from the shock front with a relative velocity of $`u_2=v_s/\chi `$.

The shock jump conditions hold for each particle species $`i`$, and raises the question whether the, still poorly understood, mechanism for shock heating the particles will also tend to rapidly equilibrate the temperatures of the various particles involved. In addition, if a large neutral fraction exists in the pre-shock material, or the shock is moving through a medium with a high metal content (e.g. ejecta), a large fraction of the electrons behind the shock may originate from ionizations. These so-called secondary electrons have to be heated by primary electrons, resulting in an overall lower electron temperature than indicated by Eq. (2), see Itoh (1984).

The earliest indications that the electron and ion temperatures are indeed not equilibrated are the relatively low electron temperatures inside SNRs, which, according to X-ray observations, in no object seem to exceed 5 keV, whereas a typical shock velocity of $`4000`$ km/s should give rise to a mean plasma temperature of 19 keV.

More direct evidence for only modest equilibration of electrons and ions behind SNR shocks comes from optical and UV spectroscopy. H$`\alpha `$ emission from non-radiative shocks in a number of SNRs is characterized spectroscopically by a narrow component and a broad component (Ghavamian et al., 2001). The narrow component is thought to come from neutral hydrogen briefly excited after being overtaken by the shock. The broad component is the result of charge exchange between neutral and ionized hydrogen. The ionized hydrogen is heated according to Eq. (2), and after charge exchange, which leaves the atom in an excited state, the thermal motion of the heated hydrogen gives rise to Doppler broadening. The width of the broad line component can be used to determine the shock velocity, if the equilibration fraction is known (see Eq. 2). The equilibration fraction is determined from the ratio between the narrow and broad component (Ghavamian et al., 2001).

Several SNRs with broadened H$`\alpha `$ emission are known (e.g. Smith et al., 1991; Ghavamian et al., 2001). In particular the work by Ghavamian et al. (2001, 2002) provides clear evidence that the youngest SNRs, i.e. the ones with large shock velocities such as SN 1006 and Tycho, have a large ratio between the electron and proton temperatures ($`T_\mathrm{e}`$/$`T_\mathrm{p}<0.07`$). Additional proof for low equilibration behind the shock front of SN 1006 comes from UV spectra obtained with the Hopkins Ultraviolet Telescope, which show broad C IV, N V and O VI lines, indicating temperature non-equilibration of these elements (Laming et al., 1996).

Very recently Vink et al. (2003) determined the amount of electron-oxygen equilibration from X-ray observations by XMM-Newton. The observation pointed at a bright, but narrow (0.4 FWHM) knot in the northeast of SN 1006, which made it a good target for the XMM-Newton reflective grating spectrometer (RGS). The size of the knot translates into a spectral broadening of $`\mathrm{\Delta }\lambda 0.05`$ Å (FWHM) or a resolution power of $`\lambda /\mathrm{\Delta }\lambda 430`$ for the dominant O VII He$`\alpha `$ line emission. However, SN 1006 itself is very extended (30) and contributes to the spectrum as well. What made it, nevertheless, a good target for measuring the thermal Doppler broadening is that this knot is compact and lies at the edge of the remnant, making it very likely that Doppler broadening is due to thermal motions rather than internal bulk velocities.

Modeling of the O VII line emission, which consists of three bright lines, revealed both a very low ionization parameter, $`\mathrm{log}(n_\mathrm{e}t)9.2`$ (cgs units), and a substantial line broadening with $`\sigma _E=3.4\pm 0.5`$ eV at 540 eV (Fig. 2). The detection of broadening is statistically very significant ($`>6.5\sigma `$ level), and corresponds to a temperature for the oxygen ions of $`kT=528\pm 150`$ keV, or, according to Eq. (2), a shock velocity of $`4100`$ km/s. This is comparable, but somewhat faster, than the shock velocity inferred from the H$`\alpha `$ emission ($`3000`$ km/s). This is not surprising as the X-ray emission comes from further downstream the shock front, where the plasma may have been shocked earlier, when the shock velocity was higher than at present. The electron temperature of the knot has been measured with the CCD instruments of XMM-Newton, and is $`T_\mathrm{e}`$$`=1.5\pm 0.2`$ keV. Calculations of the plasma history show that the difference between electron and oxygen temperature support only a small equilibration fraction near the shock front ($``$5%), and subsequent equilibration proceeds through Coulomb interactions. A conclusion that, at least for the case of SN 1006, is based on optical, UV and X-ray spectroscopy.

## 3. Cosmic ray acceleration and X-ray synchrotron emission from shell type SNRs

The most likely mechanism for cosmic ray production in SNRs is first order Fermi acceleration, according to which an energetic particle repeatedly scatters upstream and downstream due to the presence of plasma waves. For each scattering the energy of the particle in the local frame is conserved, resulting in an increase in energy after each crossing of the shock front, i.e. after each change of local frame (see e.g. Bell, 1978; Blandford & Ostriker, 1978). As the flow in the post shock region has a tendency to sweep the particles away from the shock front with a probability of $`4u_2/v_s`$, i.e. independent of the particle energy, the particle distribution is an inverse power law in momentum. The predicted power law index is $`\mathrm{\Gamma }=(2u_2u_1)/(u_1u_2)=(\chi +2)/(\chi 1)`$, which corresponds to $`\mathrm{\Gamma }=2`$ for a shock compression ratio of $`\chi =4`$.

The earliest evidence that SNRs are indeed sites of cosmic ray acceleration was their non-thermal radio emission, which for the brightest remnants could only be explained by assuming recent electron cosmic ray acceleration. The SNR radio emission is the result of synchrotron radiation caused by relativistic electrons with energies in the MeV to GeV range. Most shell-type SNRs have radio spectral indexes of $`\alpha 0.6`$ (Stephenson & Green, 2002), but some remnant have steeper spectra, like Cas A, which has $`\alpha =0.77`$, steeper than predicted by simple first order Fermi acceleration theory, $`\alpha =(\mathrm{\Gamma }1)/2=0.5`$. This is probably the result of the back reaction of the cosmic ray pressure on the shock structure. For a substantial cosmic ray pressure the flow ahead of the shock will be slowed down, as seen from the shock frame. As a result a shock pre-cursor forms and the density jump is smoothed over larger length scales (e.g. Eichler, 1979; Bell, 1987). Moreover the density jump experienced by relativistic particles depends on the particle’s gyroradius, i.e. particle energy. Very energetic particles will experience the full density jump from the undisturbed medium to the fully shocked plasma, whereas less energetic particles will only experience the difference between the shock precursor and the main shock, which may be less than the canonical factor 4 (see Eq. 2). As a result, the radio spectral index should be steeper at low radio frequencies than at high frequencies. This theory has been tested and confirmed for the Tycho and Kepler SNRs by Reynolds & Ellison (1992). Very recently Jones et al. (2003) reported that infrared observations of Cas A also indicate a flattening of the synchrotron spectrum at high frequencies.

An extreme case of cosmic ray modification of the shock structure is one in which the pressure is completely dominated by cosmic rays. In that case the shock structure is smooth instead of characterized by a distinct pressure jump. Malkov (1997) has shown that this will result in a particle index $`\mathrm{\Gamma }=1.5`$. It is not clear how physical this solution of the couple gas/cosmic ray model is, but there is no observational evidence that some SNR shocks are characterized by this extreme case of cosmic ray acceleration, which should be distinguishable by a rather flat radio index of $`\alpha =0.25`$ and large energy losses due to escaping high energy cosmic rays (for $`\mathrm{\Gamma }2`$ the highest energy cosmic rays carry most of the energy). $`\mathrm{\Gamma }=1.5`$ is also the spectral index expected for a relativistic gas with $`\chi =7`$, but Malkov (1997) has shown that for cosmic ray dominated shocks $`\mathrm{\Gamma }=1.5`$ even for $`\chi >7`$.

Although SNRs have the available energy to mark them as the primary sites of cosmic ray acceleration, one long standing questions is if they are capable to accelerate cosmic rays at least up to the observed break in the cosmic ray spectrum at $`3\times 10^{15}`$ eV (Lagage & Cesarsky, 1983). A partial answer has been provided by the recent discovery that some shell-type SNRs, which usually are dominated by thermal X-ray emission, in fact also emit X-ray synchrotron radiation. This was first discovered in SN 1006.

For a long time the featureless X-ray spectrum of SN 1006 was difficult to explain by thermal radiation, as thermal radiation produces line emission. An interpretation invoking carbon rich, hot plasma, which suppresses line emission above 0.5 keV (Hamilton et al., 1986), became untenable after ASCA observations showed that SN 1006 does display line emission typical for hot plasma, but that the overall X-ray emission is dominated by continuum emission coming from the northeastern and southwestern rims of the remnant (Koyama et al., 1995). Fig. 2 illustrates this with XMM-Newton observations. This strongly suggested that synchrotron radiation is the dominant source of X-ray emission. The spectral index of $`3`$ (Vink et al., 2000; Allen et al., 2001; Dyer et al., 2001) is steeper than expected from Fermi shock acceleration, indicating that the photons were emitted by electrons with energies close to the maximum electron cosmic ray energy.

The energy cut-off itself depends on the magnetic field strength and is approximately, for a photon cut-off energy of 1 keV,

$$E_\mathrm{e}=230/\sqrt{B_\mu }\mathrm{TeV},$$
(3)
where $`B_\mu `$ is the magnetic field strength in $`\mu `$G. As the interstellar magnetic field strength is $`6\mu `$G, and a shock compression factor of 4 may enhance it up to 20 $`\mu `$G, the implied electron cosmic ray cut-off energy is around 50 TeV. However, there is some dispute about the actual magnetic field strengths in SNRs, as will be addressed in the next section.

After SN 1006 X-ray synchrotron emission has been identified in a few other supernova remnants, notably Cas A (Allen et al., 1997), RXJ 1713.7-3946 (Koyama et al., 1997), RCW 86 (Borkowski et al., 2001; Bamba et al., 2000), and G266.2-1.2 (Slane et al., 2001). I will return to some of those objects below. Some of them have X-ray emission dominated by synchrotron radiation. However, the X-ray emission from Cas A is dominated by line emission, but a non-thermal tail of emission is seen at hard X-rays, which has been attributed to either synchrotron radiation or non-thermal bremsstrahlung (Allen et al., 1997; Favata et al., 1997; Bleeker et al., 2001; Vink & Laming, 2003).

## 4. Cosmic ray acceleration and magnetic fields in SNRs

The detection of synchrotron radiation in SN 1006 and other remnants has important consequences. On the one hand it is further proof that SNRs accelerate electrons, as the high electron cut-off indicates that the electrons can only have been accelerated recently. On the other hand, $`100`$ TeV seems to be the maximum attainable electron energy (Reynolds & Keohane, 1999), which is one order of magnitude lower than the energy at the “knee” of the cosmic ray spectrum.

The severity of the problem depends, however, on the nature of the cut-off of the electron cosmic ray spectrum. The ion and electron acceleration process is identical, apart from the initial cosmic ray injection, which means that if the energy cut-off is determined by the acceleration efficiency alone (i.e. the spectrum is age limited), the maximum electron energy is similar to the maximum ion energy. The electron energy cut-off is much less of a problem, if the magnetic fields in SNRs are greatly enhanced, and the spectra are loss limited, as losses affect predominantly the electron cosmic ray spectra. The relative importance of loss limitations to age limitations depends on the mean magnetic field strength of the medium through which the electrons move.

An additional effect of higher magnetic fields is that the acceleration is more efficient and age limitations of the ion cosmic ray spectrum can result in maximum energies beyond the “knee” energy. For that reason Biermann & Cassinelli (1993) proposed that most of the cosmic rays around the “knee” are produced in SNRs moving through the stellar winds of their progenitors, which may have enhanced the local magnetic field with respect to the average galactic magnetic field. An alternative idea is that the plasma waves generated by the cosmic rays result in non-linear behavior, enhancing the background magnetic field (Bell & Lucek, 2001). In other words, once cosmic rays are present they may significantly enhance further cosmic ray acceleration.

Vink & Laming (2003) proposed a way of measuring the magnetic field strength near the shock front of Cas A. Independently Bamba et al. (2003) proposed the same method for SN 1006. As an example I will concentrate on Cas A.

The main idea is that the narrow rim of X-ray continuum emission surrounding Cas A is indeed X-ray synchrotron emission (Gotthelf et al., 2001). Near the shock front, electrons are continuously being accelerated, but as soon as electrons move away from the shock front, radiative losses make that the electron population cuts off at progressively lower energies, until it does no longer produce X-ray emission. The plasma moves downstream of the shock front with a velocity of $`\frac{1}{4}v_s`$. So the timescale, $`\tau `$, for radiative losses is coupled to a length scale by $`\mathrm{\Delta }r=\frac{1}{4}v_s\tau `$.

The shock velocity of Cas A has been measured to be $`5000`$ km/s (Vink et al., 1998; Delaney & Rudnick, 2003), whereas the width of the rim, i.e. the length scale $`\mathrm{\Delta }r`$, is between 1.5′′ to 4′′. For a distance of Cas A of 3.4 kpc (Reed et al., 1995), this corresponds to loss times of 18 yr to 50 yr. The loss times are inversely proportional to $`B^2E_\mathrm{e}`$, whereas the photon energy corresponding to the electron energy ($`E_\mathrm{e}`$) scales as $`BE_\mathrm{e}^2`$. Combining the measured loss times and the fact that X-ray synchrotron radiation is observed around 5 keV allows the determination of the typical magnetic field strength near the shock front, which is about $`10^4`$ G (Fig. 3), more than an order of magnitude higher than the average galactic field strength. Note that this approach assumes that the electrons at those high energies are still coupled to the plasma, i.e. the mean free particle path, $`\lambda _{mfp}`$, should be smaller than the rim width. In other words, the diffusion coefficient needs to be small. As this is a prerequisite for acceleration at those energies, this seems a reasonable assumption. In order to be self-consistent the mean free path in a radial direction should be larger than the gyroradius, $`r_g`$, for the derived magnetic. The mean free path is often parameterized as $`\lambda _{mfp}=\eta r_g`$, with $`\eta 1`$ (c.f. Reynolds, 1998). The width of the rim implies $`\lambda _{mfp}<10^{17}`$ cm2, or $`\eta <100`$ for $`B=0.1`$ mG and $`E`$50 TeV.

The relatively high magnetic field in Cas A and SN 1006 suggests that rapid cosmic ray acceleration of electrons and ions is possible, and, as the electron population is loss limited, the ion cosmic rays can in principle be accelerated up to, or beyond the “knee”. In fact, close to the shock front even electrons may obtain high energies, as the electron energy indicated by Fig. 3 applies to the energy of the population downstream of the shock.

Vink & Laming (2003) also used another method for estimating the average magnetic field strength inside Cas A, which consisted in comparing an upper limit on the non-thermal bremsstrahlung above 100 keV, assuming it is caused by the low energy part of the electron cosmic ray spectrum, and the overall radio flux (c.f. Cowsik & Sarkar, 1980; Atoyan et al., 2000).

The bremsstrahlung normalization scales with $`n_\mathrm{e}\mathrm{\Sigma }_in_\mathrm{i}Z_i^2`$, whereas the radio synchrotron emission scales with $`n_\mathrm{e}B_{}^{(\mathrm{\Gamma }+1)/2}`$. As $`\mathrm{\Sigma }_in_\mathrm{i}Z_i^2`$ is approximately known, the average relativistic electron density and magnetic field strength can be inferred. The measurements imply $`B>0.5`$ mG, higher than the magnetic field strength near the shock. The reason is presumably magnetic field enhancements due to turbulence associated with the contact discontinuity of shocked ejecta and swept up circumstellar matter. This may explain why the brightest radio structure in Cas A is a shell coinciding with the ejecta shell, instead of with the outer rim.

The high overall magnetic field makes it unlikely that X-ray synchrotron radiation is coming from the bright radio and X-ray shell of Cas A; inside the bright shell the synchrotron loss time for the X-ray emitting electrons is only of the order of 10 yr, much less than the age of Cas A. As some hard X-ray emission seems to be associated with that bright shell (Bleeker et al., 2001), it is likely that another emission mechanism is contributing to the hard X-ray emission. Laming (2001a, b) has proposed that internal shocks produce lower hybrid plasma waves, which are responsible for accelerating electrons up to $`100`$ keV, where the observed X-ray spectrum has a cut-off. This predicted spectrum was found to be consistent with the hard X-ray spectrum obtained with a deep, 500 ks, BeppoSAX observation of Cas A (Vink & Laming, 2003, Fig. 4).

## 5. TeV emission: Evidence for hadronic cosmic ray acceleration or not?

So far our discussion has focused on electron cosmic rays, as they are most easily observed through synchrotron radiation and bremsstrahlung. However, hadronic cosmic rays are arguably the most important part of the cosmic rays, as they constitute the cosmic rays best observed at earth, and the ratio of hadronic to leptonic (electron) cosmic rays is roughly $`100:1`$. The direct observation of ultra-relativistic in SNRs is therefore an important part of the proof that shock acceleration in SNRs is responsible for the observed cosmic ray spectrum up to the “knee”. This goal is coming nearer, or may even have been reached, with the recent advances made with telescopes observing Cherenkov radiation in the earth atmosphere caused by TeV $`\gamma `$-rays.

Although the electromagnetic radiation from relativistic ions is weak, ions colliding with background ions betray themselves by producing pions ($`\pi ^+,\pi ^{},\pi ^0`$), which decay quickly into muons, positrons, electrons and neutrinos, except for $`\pi ^0`$, which decays in 99% of the cases into two photons with energies of 68 GeV each in the $`\pi ^0`$ rest frame. The threshold for pion creation is 290 MeV. Well above the production threshold 30% of the particle energy will be transferred to the pion, and the emerging $`\gamma `$-ray spectrum will have a power law slope similar to the index of the particle spectrum.

Pion decay is the dominant component of the galactic background $`\gamma `$-ray emission above 1 GeV (Hunter et al., 1997). However, the situation for $`\gamma `$-ray emission from SNRs is not so clear. First of all, CGRO-EGRET detected only few of the known SNRs, most of which turn out to be older SNRs, like IC 443 and $`\gamma `$-Cygni (see Torres et al., 2003, for a review). Secondly, inverse Compton scattering of background photons by ultra-relativistic electrons is a very plausible, alternative, explanation for TeV $`\gamma `$-ray emission from young SNRs.

The first shell-type SNR detected above energies of $`1.7`$ TeV is SN 1006 (Tanimori et al., 1998). The discovery was made with the CANGAROO telescope. Less than two years before this detection the synchrotron nature of the X-ray emission had been established, so it was natural to suggest that the same electrons responsible for the X-ray synchrotron radiation were responsible for the $`\gamma `$-ray emission by means of inverse Compton scattering of cosmic microwave background (CMB) photons. As the normalization for inverse Compton scattering scales with $`n_en_{ph,CMB}`$, and synchrotron radiation with $`n_eB^{(\mathrm{\Gamma }+1)/2}`$, with $`n_{ph,CMB}`$ the photon density, one can in principle infer the magnetic field from the ratio of the hard X-ray and TeV $`\gamma `$-ray emission, which turns out to be $`B=6.5\pm 2\mu `$G (Tanimori et al., 1998).

However, this is only valid if the TeV emission is indeed due to inverse Compton scattering. Berezhko et al. (2002) recently argued that from a theoretical point of view a higher magnetic field is preferred, a point of view that seems to be supported by the magnetic field measurement by Bamba et al. (2003). This would mean that the relativistic electron density necessary to explain the X-ray synchrotron spectrum can be lower, which decreases the inferred inverse Compton scattering contribution and the maximum electron energy (Eq. 3). Berezhko et al. (2002) argue that therefore pion decay is the likely origin of the TeV $`\gamma `$-ray emission from SN 1006. An additional argument against dominant inverse Compton scattering above 1 TeV is that the $`\gamma `$-ray emission seems to come from the northeastern side of SN 1006, whereas the synchrotron radiation, caused by the same electrons that upscatter background photos, is observed from both the northeastern and southwestern rims of the remnants (Fig. 2). In the case of pion decay this asymmetry could be explained by the presence of a density enhancement in the vicinity of the northeastern rim.

Two other young remnants have been discovered by Cherenkov telescopes. Cas A was detected at the 4.9$`\sigma `$ level with the HEGRA telescope (Aharonian et al., 2001). For Cas A the case for pion decay as the origin for the TeV emission is even more compelling than for SN 1006, although still not conclusive. As explained in the previous section, the average magnetic field inside Cas A is $`>0.5`$ mG, which means that the electron synchrotron loss time is short, making it unlikely that most electron populations inside Cas A extend all the way up to 100 TeV, even if most of the hard X-ray emission is synchrotron radiation. This makes inverse Compton emission a less likely mechanism for producing TeV emission. However, as Vink & Laming (2003) have shown (see previous section), the magnetic field near the shock front is $`0.1`$ mG, and the relativistic electron density is somewhat higher. Moreover, Cas A itself is a bright infrared source, so that more background photons are present, which have on average a higher energy as well (Atoyan et al., 2000). This gives rise to the additional complication of an anisotropic photon field. Despite these reservations, Vink & Laming (2003) argue that pion decay and not inverse Compton emission is the likely source for the TeV emission, but not by a wide enough margin to be completely confident that the inverse Compton scattering interpretation can be excluded (Fig. 5).

Of the three remnants detected above 1 TeV, RXJ 1713.7-3946 is the most surprising one. The remnant was discovered with ROSAT (Pfeffermann & Aschenbach, 1996), and ASCA spectra showed that the X-ray emission was dominated by non-thermal radiation (Koyama et al., 1997). The detection of TeV emission was made with the CANGAROO telescope (Muraishi et al., 2000).

One of the interesting features of RXJ 1713.7-3946 is that it is probably an old remnant, judged by its size of $`34`$ and the likely association with a molecular cloud at a distance of roughly 6 kpc (Slane et al., 1999), whereas it is still a source of X-ray synchrotron radiation. The cosmic ray spectra in older remnants is likely to be loss limited, but the maximum electron energy is proportional to the shock velocity for loss limited models (Reynolds, 1998), which means that older remnants are not expected to produce considerable X-ray synchrotron radiation.

As pointed out by Slane et al. (1999) some of this problem may be alleviated, if RXJ 1713.7-3946 is a remnant evolving in stellar wind cavity. In that case the shock velocity remains high until it encounters the cavity wall, after which the remnant is believed to brighten rapidly and radiate most of its energy away in a relatively short time. Another surprising feature of this remnant is that no evidence has turned up yet that there is thermal X-ray emission coming from this remnant, unlike SN 1006. This is only possible, if either the electron temperature is very cool (below 0.1 keV), or the density inside the remnant is very low, (cf. SN 1006 which has $`n_\mathrm{e}0.1`$ cm-3).

The fact that the X-ray emission from RXJ 1713.7-3946 is synchrotron dominated suggested, as it did for SN 1006, that the emission above 1 TeV is caused by inverse Compton emission. However, recent observations with CANGAROO showed that the 1-10 TeV $`\gamma `$-ray spectrum, in combination with the X-ray spectrum, is incompatible with inverse Compton emission, and that it is therefore likely that the $`\gamma `$-ray emission is the result of pion decay (Enomoto et al., 2002). The CANGAROO $`\gamma `$-ray spectrum would therefore be the first direct evidence that nucleons are accelerated by SNRs. However, shortly after the announcement that ultra-relativistic ions had finally been observed in SNRs both Reimer & Pohl (2002) and Butt et al. (2002) pointed out that the spectrum of a nearby and possibly associated CGRO-EGRET source, 3EG J1714-3857, was not consistent with a pion decay spectrum from RXJ 1713.7-3946, as the CGRO-EGRET data points lie below the pion decay model. Although an obvious objection is that 3EG J1714-3857 and RXJ 1713.7-3946 do not spatially coincide, it is hard to believe that CGRO-EGRET detected 3EG J1714-3857, whereas a brighter nearby source should have been missed.

## 6. The cosmic ray injection efficiency

One topic not yet addressed in this article is the injection process and injection efficiency of cosmic rays. Our knowledge of the injection process, especially for electrons, is limited, but some progress has been made over the last decade (e.g. Bykov & Uvarov, 1999), and computers are now powerful enough to simulate the microscopic processes with particle in cell methods (Schmitz et al., 2002a, b).

The cosmic ray efficiency is also of major importance for the hydrodynamical modeling of supernova remnants, because, if a major fraction of the shock energy goes into accelerating particles, the plasma may become cosmic ray dominated, possibly increasing the shock compression ratio and lowering the mean gas temperature (Eq. 2). If the acceleration is very efficient, high energy particles may diffuse away from the remnant altogether, making the shock essentially radiative, under which conditions the compression ratio may even exceed a factor of 7 (e.g. Decourchelle et al., 2000). This has in fact been claimed for the Small Magellanic Cloud remnant 1E 0102.2-7219, for which Hughes et al. (2000) has shown that the electron temperature is lower than can be expected, based on the observed shock velocity, even when allowing for slow electron-ion temperature equilibration.

However, it should be pointed out that in most cases, even if very efficient cosmic ray acceleration takes place, it is not likely that the shocks becomes radiative or even acts like a relativistic gas with $`\gamma =4/3`$. The reason is that most cosmic ray spectral indexes are observed to be $`\mathrm{\Gamma }>2`$, for which the total cosmic ray energy333This follows from the integration over the cosmic ray spectrum, $`_{E_0}^{\mathrm{}}\kappa E^\mathrm{\Gamma }E𝑑E`$. is dominated by non-relativistic cosmic rays. This can be illustrated for Cas A, as the magnetic field strength measurement allows for the estimation of the electron cosmic ray density normalization, which turns out to be $`\kappa 5\times 10^{11}`$ (cgs).444See the previous footnote. $`\kappa `$ normalizes the total number of electron cosmic rays per cm-3 (for Cas A $`\mathrm{\Gamma }=2.55`$). If we integrate the number density, assuming a power law spectrum in momentum, we find that the number of particles above $`p_0=0.01m_\mathrm{p}c`$ (somewhat arbitrary, but following Bell, 1987) is $`8\times 10^4(\zeta +1)`$, with $`\zeta `$ the ratio of nucleonic to electron cosmic rays. As the thermal particle density in Cas A is about 20 cm-3 and $`\zeta `$ is generally believed to be in the range 1-100, we find a cosmic ray injection efficiency between $`\varphi =4\times 10^5`$ to $`4\times 10^3`$. Interestingly, this range in injection efficiencies coincides with the transition from thermally dominated pressure to cosmic ray dominated pressure (Bell, 1987), suggesting that perhaps the acceleration process is indeed self-regulating. Note, however, that for the steep cosmic ray spectrum of Cas A the cosmic ray pressure is dominated by non-relativistic nucleonic cosmic rays.

Our understanding would be greatly advanced if we could observe the injection spectrum directly. Vink et al. (1997) suggested that the the X-ray spectrum from certain parts of RCW 86, which are characterized by featureless spectra except for the presence of iron line emission at 6.4 keV, could be explained by a low temperature electron gas with a non-thermal bremsstrahlung tail, which also causes line emission from underionized iron (see also Vink et al., 2002). This non-thermal tail may be caused by the electron cosmic ray injection spectrum. However, Rho et al. (2002) have shown that such a bi-modal electron distribution will quickly evolve toward a Maxwellian distribution, resulting in more line emission. They prefer to explain the observed spectral features by assuming that the continuum is X-ray synchrotron radiation, whereas part of the plasma inside the remnants consists of a hot, low density, pure iron gas (see also Bamba et al., 2000; Borkowski et al., 2001; Bykov, 2002).

## 7. Concluding remarks

I have reviewed the current observational evidence for efficient cosmic ray acceleration by the collisionless shocks in supernova remnants. For a long time the evidence consisted only of the observed non-thermal radio spectra, but the last decade the more important evidence has come from X-ray and $`\gamma `$-ray observations, which indicate that electrons are accelerated up to 100 TeV. However, direct evidence that SNRs accelerate cosmic rays up to, or beyond, the “knee”, and that both electrons and ions are accelerated, is not unambiguous yet. This situation is likely to improve considerably over the next five years, with new space missions like INTEGRAL and GLAST, and the completion of new Cherenkov telescopes, such as HESS.

## Acknowledgements

It wish to express my gratitude to Martin Laming, Johan Bleeker, and Jelle Kaastra, my collaborators on some of the topics discussed here. This work is supported by the NASA through Chandra Postdoctoral Fellowship Award Number PF0-10011 issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-39073.