text:
# Terawatt-scale sub-10-fs laser technology – key to generation of GW-level attosecond pulses in X-ray free electron laser

## 1 Introduction

At the start of this century, we have seen a revolution in synchrotron radiation source intensities. This revolution stemmed from the technique of free electron laser (FEL) based on self-amplified spontaneous emission (SASE), combined with recent progress in accelerator technologies, developed in connection with high-energy linear colliders. In 2001, the VUV FEL at the TESLA Test Facility at DESY (Hamburg, Germany) has successfully demonstrated saturation from 82 nm to 125 nm with GW-level peak power and pulse duration down to 40 fs . It is the first result from this device that Wabnitz et al. reported in . They illuminated xenon clusters with high-intensity ($`10^{14}`$ W/cm2) VUV FEL pulses and observed an unexpectedly strong absorption of the VUV radiation. Such a highly nonlinear optical interaction between light and matter at VUV wavelength range has never been seen before and these fascinating results show the potential of this new class of light sources for scientific research. While modern third generation synchrotron light sources are reaching their fundamental performance limit, recent successes in the development of the VUV FEL at DESY have paved the way for the construction of the novel type of light source which will combine most of the positive aspects of both a laser and a synchrotron. Starting in 2004, the phase 2 extension of TTF will deliver FEL radiation down to the soft X-ray spectral range with minimum wavelength of about 6 nm in the first harmonic and reaching into ”water window” in the second harmonic .

Recently the German government, encouraged by these results, approved funding a hard X-ray SASE FEL user facility – the European X-Ray Laser Project XFEL. The US Department of Energy (DOE) has given SLAC the goahead for engineering design of the Linac Coherent Light Source (LCLS) device to be constructed at SLAC. These devices should produce 100 fs X-ray pulses with over 10 GW of peak power . These new X-ray sources will be able to produce intensities of the order of $`10^{18}\mathrm{W}/\mathrm{cm}^2`$. The main difference between the two projects is the linear accelerator, an existing room temperature linac for LCLS at SLAC, and a future superconducting linac for European XFEL. The XFEL based on superconducting accelerator technology will make possible not only a jump in a peak brilliance by ten orders of magnitude, but also an increase in average brilliance by five orders of magnitude compared to modern 3rd generation synchrotron radiation facilities. The LCLS and European XFEL projects are scheduled to start operation in 2008 and 2012, respectively.

The motivation for the development of XFELs was recently described in detail in . The discussion in the scientific community over the past decade has produced many ideas for novel applications of the X-ray laser. Brilliance, coherence, and timing down to the femtosecond regime are the three properties which have the highest potential for new science to be explored with an XFEL. In its initial configuration the XFEL pulse duration is about 100 femtoseconds. Even though this is a few hundreds times shorter than in third generation light sources, it can probably be further reduced to about 10 femtoseconds . A novel way to generate sub-10 fs x-ray pulses – the slotted spoiler method (P. Emma, 2003) has been proposed recently. This method is based on spoiling the beam phase density in a part of the electron bunch so that this part will not lase, while preserving lasing in a short length of the bunch. The FEL performance of the spoiled beam approach was computed using the time-dependent GENESIS simulation. It has been shown that it is possible to produce X-ray pulses with duration of 3-4 fs FWHM for nominal LCLS bunch compression parameters .

Femtosecond-resolution experiments with X-rays can possibly show directly how matter is formed out of atoms. In fact, X-ray pulse duration even shorter than one femtosecond may be useful for many scientific applications. The reason is that phenomena inside atoms occur on sub-femtosecond timescale. Generating single attosecond $`0.1`$ nm X-ray pulses is one of the biggest challenges in physics. The use of such a tool will enable to trace process inside the atoms for the first time. If there is any place where we have a chance to test the main principles of quantum mechanics in the pure way, it is there.

The interest in the science with attosecond pulses is growing rapidly in the very large laser community. This community is familiar with attosecond pulses of light at photon energies up to about 0.5 keV (3 nm). This is achieved by focusing a fs laser into a gas target creating radiation of high harmonics of fundamental laser frequency. The key to these developments was the invention of laser systems delivering pulses in the range of 5 fs with pulse energies higher than a fraction of mJ. This approach produced the first measurable XUV pulses in the 200 as regime . In principle, table-top ultra-fast X-ray sources have the right duration to provide us with a view of subatomic transformation processes. However, their power and photon energy are by far low. The XFEL is ideally suited the purpose of this emerging field of science. Recently an approach for the generation of attosecond pulses combining fs quantum laser and harmonic cascade (HC) FEL scheme  was proposed in . The HC FEL scheme has the potential to produce coherent light down to wavelengths of a few nm in an undulator sequence . The analysis presented in  shows that this technique has potential to produce 100 as long radiation pulses with MW-level of output power down to 1 nm wavelength.

The X-ray SASE FEL holds a great promise as a source of radiation for generating high power, single attosecond pulses. What ultimately limits the XFEL pulse duration? Since the temporal and spectral characteristics of the radiation field are related to each other through Fourier transform, the bandwidth of the XFEL and the pulse duration cannot vary independently of each other. There is a minimum duration-bandwidth product (uncertainty principle). The shortest possible X-ray pulse duration generated by XFEL is limited by the intrinsic bandwidth of the SASE process. In the case of the European XFEL and the LCLS, the FWHM bandwidth near saturation (at 0.1 nm) is about 0.1%, indicating a 300-as coherence time determined by the bandwidth product. Recently a scheme to achieve pulse durations down to 400-600 attoseconds at a wavelength of 0.1 nm has been proposed . It uses a statistical properties of SASE FEL high harmonic radiation. The selection of a single 10-GW level attosecond pulses is achieved by using a special trigger in data acquisition system. A promising scheme for attophysics experiments using this approach has been studied and could be implemented in the XFEL design .

In this paper we propose a new method allowing to reduce the pulse length of the X-ray SASE FEL to the shortest conceptual limit of about 300 as. It is based on the application of a sub-10-fs laser for slice energy modulation of the electron beam, and application of a crystal monochromator for the selection of single attosecond pulses with GW-level output power.

## 2 The principle of attosecond techniques based on the use of XFEL combined with fs quantum laser

A basic scheme of the attosecond X-ray source is shown in Fig. 1. An ultrashort laser pulse is used to modulate the energy of electrons within the femtosecond slice of the electron bunch at the seed laser frequency. The seed laser pulse will be timed to overlap with the central area of the electron bunch. It serves as a seed for a modulator which consists of a short (a few periods) undulator. Following the energy modulator the beam enters the X-ray undulator. The process of amplification of radiation in this undulator develops in the same way as in a conventional X-ray SASE FEL: fluctuations of the electron beam current serve as the input signal . The proposed scheme for the generation of attosecond pulses is based on frequency-chirping the SASE radiation pulse. When an electron beam traverses an undulator, it emits radiation at the resonance wavelength $`\lambda =\lambda _\mathrm{w}(1+K^2/2)/(2\gamma ^2)`$. Here $`\lambda _\mathrm{w}`$ is the undulator period, $`mc^2\gamma `$ is the electron beam energy, and $`K`$ is the undulator parameter. The laser-driven sinusoidal energy chirp produces a correlated frequency chirp of the resonant radiation $`\delta \omega /\omega 2\delta \gamma /\gamma `$. After the undulator, the radiation is passed through a crystal monochromator which reflects a narrow bandwidth. Since the radiation frequency is correlated to the longitudinal position within the beam, a short temporal radiation pulse is transmitted through the monochromator.

Recent technological advances in ultrafast optics have permitted the generation of optical pulses comprising only a few oscillation cycles of the electric and magnetic fields. The pulses are delivered in a diffraction-limited beam . The combination of a X-ray SASE FEL and a few-cycle laser field techniques is very promising. Our concept of an attosecond X-ray facility is based on the use of a few-cycle optical pulse from a Ti:sapphire laser system. This optical pulse is used for the modulation of the energy of the electrons within a slice of the electron bunch at a wavelength of 800 nm. Due to the extreme temporal confinement, moderate optical pulse energies of the order of a few mJ can result in an electron energy modulation amplitude larger than 30-40 MeV. In few-cycle laser fields high intensities can be ”switched on” nonadiabatically within a few optical periods. As a result, a central peak electron energy modulation is larger than other peaks. This relative energy difference is used for the selection of SASE radiation pulses with a single spike in the time domain by means of a crystal monochromator. A schematic, illustrating these processes, is shown in Fig. 3. Many different output fields can be realized by using different spectral windowing. For instance, it is possible to generate a sequence of 300-as X-ray pulses, separated by $`T_0`$ (or $`2T_0`$), where $`T_0`$ is the Ti:sapphire laser oscillation period. Such operation of the attosecond X-ray source is illustrated in Fig. 3.

The discussion in this paper is focused on the parameters for the European XFEL operating in the wavelength range around 0.1 nm . Optimization of the attosecond SASE FEL has been performed with the three-dimensional, time dependent code FAST  taking into account all physical effects influencing the SASE FEL operation (diffraction effects, energy spread, emittance, slippage effect, etc.). In our scheme the separation of the frequency offset from the central frequency by a monochromator is used to distinguish the 300-as pulse from the 100 fs intense SASE pulse. The monochromatization is straightforward: for the 0.1 nm wavelength range, Bragg diffraction is the main tool used for such purposes. In this case, one has to take care that the short pulse duration is preserved. Transmission through the monochromator will produce some intrinsic spreading of the pulse, and the minimum pulse duration which may be selected by this method is limited by the uncertainty principle. The number of possible reflections which provide the required spectral width is rather limited. We are discussing here only Ge crystals, which have the largest relative bandwidth. This is an important feature which ensures the preservation of the single-spike pulse duration. In its simplest configuration the monochromator consists of Ge crystal diffracting from the (111) lattice planes. We show that it is possible to produce X-ray pulses with FWHM duration of 300 as. In some experimental situations this simplest configuration of monochromator is not optimal. In particular, our study has shown that the maximum contrast of the attosecond X-ray pulses does not exceed 80% and is due to the long tail of the intrinsic crystal reflectivity curve. The obvious and technically possible solution of the problem of contrast increase might be to use a premonochromator. One can align the premonochromator so that the main peak of the spectrum is blocked.

## 3 Generation of attosecond pulses from XFEL

In the following we illustrate the operation of an attosecond SASE FEL for the parameters close to those of the European XFEL operating at the wavelength 0.1 nm . The parameters of the electron beam are: energy 15 GeV, charge 1 nC, rms pulse length 25 $`\mu `$m, rms normalized emittance 1.4 mm-mrad, rms energy spread 1 MeV. Undulator period is 3.4 cm.

### 3.1 Slice modulation of the electron beam

The parameters of the seed laser are: wavelength 800 nm, energy in the laser pulse 2–4 mJ, and FWHM pulse duration 5 fs (see Fig. 4). The laser beam is focused onto the electron beam in a short undulator resonant at the optical wavelength of 800 nm. Parameters of the undulator are: period length 50 cm, peak field 1.6 T, number of periods 2.

Optimal conditions of the focusing correspond to the positioning of the laser beam waist in the center of the undulator. In laser pulses comprising just a few wave cycles, the amplitude envelope and the carrier frequency are not sufficient to characterize and control laser radiation, because the evolution of the light field is also influenced by a shift of the carrier wave with respect to the pulse peak . Recently, the generation of intense, few-cycle laser pulses with a stable carrier envelope phase $`\phi _0`$ was demonstrated . Let us consider the principle question for the design of few-cycle pulse experiments: how does the pulse phase behave during linear propagation? In order to answer this question, we can calculate the evolution of few-cycle pulses in vacuum, which is most conveniently described by a parabolic wave equation and by starting from a Gaussian initial spatial pulse profile . Choosing the initial pulse phase to be $`\phi _0`$ at the beam waist, one reveals that the carrier envelope phase in the far field $`\phi (\mathrm{})=\phi _0\pi /2`$ undergoes a phase shift due to the Guoy phase shift $`\pi /2`$ . Note that the Guoy phase shift and all the other changes experienced by the pulse during propagation do not depend on the initial phase $`\phi _0`$.

For an attosecond X-ray source it is of great interest to maximize the central peak energy offset, which depends sensitively on the absolute phase of the seed laser pulse $`\phi _0`$. We start with an illustration of the few-cycle-driven energy modulation assuming that the peak electric field appears at the peak of the envelope when the laser pulse passes the undulator center (i.e. $`\phi _0=0`$ at the Gaussian beam waist). The interaction with the laser light in the undulator then produces a time dependent electron energy modulation as shown in Fig. 5. For the laser (FWHM) pulse duration of 5 fs at a laser pulse energy 2-4 mJ, we expect a central peak energy modulation 30-40 MeV.

### 3.2 Monochromator

The width of the spectral distribution of the SASE radiation will be determined by the frequency chirp, provided the frequency chirp is larger than FEL bandwidth. A monochromator may be used to select the pulses of short duration, due to correlation between frequency and longitudinal position in the radiation pulse. For 12 keV photons, we consider Bragg diffraction in crystals as a method of bandwidth selection. Special attention is called to the fact that the relative spectral width for the given Bragg reflection is independent of the wavelength or glancing angle of X-rays and is given merely by properties of the crystal and the reflecting atomic planes. In particular, it implies that the choice of a crystal and reflecting atomic planes determines the spectral resolution. For example, one can consider Si(111) crystals, which have a FWHM bandwidth of $`\mathrm{\Delta }\lambda /\lambda =1.3\times 10^4`$, or Ge(111)crystals, which have a FWHM bandwidth of $`\mathrm{\Delta }\lambda /\lambda =3.4\times 10^4`$. Monochromators at synchrotron beam lines are most commonly fabricated from silicon. The reason is that the semiconductor industry has created a huge demand for defect-free , perfect single crystals. Silicon is by no means the only choice, and in recent years diamond has become a popular alternative, due to the fact it has the highest thermal conductivity and low absorption.

An attosecond X-ray source requires a relatively broadband monochromator. The larger the monochromator bandwidth is, the shorter the minimal pulse duration than can be extracted. We are discussing here only Germanium single crystals which have the largest relative bandwidth. Although Ge is not as perfect as silicon or diamond, sufficiently large perfect Ge crystals are available today. For 12 keV photons Bragg peaks of Ge crystals have reflectivities of approximately 75%. Figure 6 gives an example of a reflectivity curve for a thick absorbing crystal. The drawing of Fig. 6 shows several interesting features. The shape is asymmetric and is due to absorption effect. The tails of the reflectivity curve decrease as $`(\mathrm{\Delta }\lambda /\lambda )^2`$. It should be pointed out that the tail of reflectivity curve plays important role in the operation of the attosecond X-ray source, and this characteristic of spectral window and attosecond pulse contrast are ultimately connected. Good crystal quality is required for high resolving power. Similarly, a good resolving power requires a collimated beam. In fact, the angular spread of an insufficiently collimated beam negatively affects the wavelength resolution just like poor crystal quality. The Ge monochromator angular acceptance is of order 50 $`\mu `$rad for a wavelength 0.1 nm, and is well matched to the natural opening angle, (1 $`\mu `$rad), of an XFEL source. Therefore, a crystal monochromator works better with XFEL radiation than with conventional synchrotron source.

Besides the crystal quality, other factors must be considered in selecting a scheme for a monochromator. The monochromator crystal must be thermally stable and capable of being exposed to XFEL output radiation with limited radiation damage. We have chosen the one Ge (111) crystal scheme of the X-ray monochromator with silicon premonochromator, which withdraws the major heat load from the actual short pulse selection Ge monochromator. We consider Laue diffraction in Si crystals as a method of bandwidth selection in premonochromator. In the pre-monochromator part one can use ten Si crystal plates of $`15\mu `$m thickness and the reflection Si(111). Given the crystal plate is perfect, it reflects 90% of the incident X-rays within a band of $`\mathrm{\Delta }\lambda /\lambda 10^4`$. One can align the Si plates so that the main peak of the spectrum is blocked. The radiation power which reaches the broadband (Ge) monochromator crystal is 10% of the initial value. Only 30% of the offset frequency radiation is absorbed in ($`15\times 10=150\mu `$m thickness) premonochromator totally, and the rest passes through. Another advantage of the premonochromator is the possibility to increase the contrast of output attosecond X-ray pulses.

### 3.3 Output characteristics of attosecond FEL

A complete description of the X-ray FEL can be performed only with three-dimensional, time-dependent numerical simulation code. Since amplification process starts from shot noise, properties of a single-spike selection should be described in statistical term. The statistics of concern are defined over an ensemble of radiation pulses.

In the present scheme an electron beam with slice modulation of the energy passes through the undulator and produces SASE radiation. Since only a small fraction of the bunch is modulated (10 fs versus 200 fs of FWHM electron pulse duration, see Fig.5), the total energy in the radiation pulse remains approximately the same as in the case of nonmodulated electron beam (see Fig. 7), and saturation is achieved at an undulator length of about 120 m. Figure 8 shows typical temporal and spectral structure of the radiation pulse at saturation. In the present numerical example the central part of the electron bunch was subjected to the slice energy modulation, and one can notice its clear signature in the temporal structure. Figure 9 shows an enlarged view of the central part of the radiation pulse. The dotted lines in this figure show the initial energy modulation of the electron beam. The temporal structure of the radiation pulse has a clear physical explanation. The FEL process is strongly suppressed in the regions of the electron bunch with large energy chirp, and only regions of the electron bunch with nearly zero energy chirp produce radiation. From a physical point of view each part of the bunch near the local extremum of the energy modulation can be considered as an isolated electron bunch of short duration. At the chosen parameters of the system its duration is about 300 attosecond which is about of coherence time. Thus, it is not surprising that only a single radiation spike is produced by each area near the local extremum. An asymmetry of averaged power (at symmetric energy modulation) is due to the nonsymmetry of the FEL process with respect to the sign and the value of the energy chirp. In particular, spikes related to the negative of energy offset have a higher amplitude. This is a typical nonlinear effect allowing to prolong interaction of the radiation pulse with the electron pulse at the ”correct” sign of the energy chirp. The temporal structure of the radiation pulse nearly repeats the temporal structure of the energy modulation. Slight deviations from the periodic structure are also due to the nonsymmetry of the FEL process with respect to the sign of the energy chirp. In particular, spikes related to the negative energy offset slip more forward than those related to the positive energy offset.

Let us turn back to the main subject of our study, i.e. to the production of attosecond pulses. The lower plot in Fig. 8 shows the total radiation spectrum of the radiation pulse. At this scale a signature of the slice energy modulation can hardly be seen. In Fig: 10 we present the tails of the spectrum at an enlarged scale. Each of three clearly visible bumps in the averaged spectrum corresponds to a local extremum of the energy offset shown in Fig. 5. The bump marked as $`M_1`$ corresponds to the central peak energy offset. The bump $`M_2`$ corresponds to the neighboring two positive energy offsets. The bump $`M_3`$ comes from the areas of the electron bunch with negative energy offset. The single-shot spectrum (shown as grey line) exhibits an oscillatory behavior near bumps $`M_2`$ and $`M_3`$. That is due to an interference of two radiation wavepackets with close frequencies coming from different parts of the electron bunch. Other maxima can be hardly distinguished, since they are located within the bandwidth of the main spectrum.

Figures 3, 3, and 10 give a clear idea about separation of the attosecond radiation pulses. Positioning of the monochromator to different maxima of the spectrum allows us to select single pulse, or a two pulse sequence of attosecond duration. The calculation involves the following steps. The FEL simulation code produces 3-D arrays for the radiation field in the near zone. This field is recalculated into the far field zone, and is subjected to the Fourier transform. The latter result is convoluted with the reflectivity function of Ge(111) monochromator (see Fig. 6), and is subjected to inverse Fourier transform giving temporal structure of the radiation pulse behind the monochromator.

By selecting the frequency offset of the monochromator to the position marked as $`M_1`$ in Fig. 10, we select single pulses. Their properties are illustrated with Fig. 11. An analysis of single pulses shows that their pulse duration is about 300 as, the average power has GW-level, and the radiation pulse energy is about a $`\mu `$J. The larger width of the averaged curve is partially due to shot-to-shot fluctuations of the position of the radiation pulse (a fraction of coherence time). Note that shot-to-shot fluctuations of the radiation energy after monochromator are suppressed significantly due to ultrashort duration of the lasing fraction of the electron bunch . An advantage of single-pulse selection is the small background from the main radiation pulse due to a large offset from the resonant frequency.

By positioning of the monochromator central frequency to the spectrum bumps $`M_2`$ or $`M_3`$ one can select a two pulse sequence as illustrated in Figs. 12 and 13. Two pulses are separated by two or one oscillation period of optical laser depending on the choice of the monochromator tuning. Note that due to the statistical nature of the SASE process the time jitter between two pulses is about 200 as, a fraction of the coherence time. One should not wonder that pulse amplitudes differ visibly for the case of pulse separation by one laser oscillation period (see Fig. 12). As mentioned above, this is a typical nonlinear effect related to the sensitivity of the FEL process to the sign and the value of the energy chirp. Although the energy modulation amplitude is the same in both maxima, the shape of the energy chirp is asymmetric.

## 4 Discussion

Successful operation of the attosecond XFEL requires the fulfillment of several requirements. The requirement that the SASE FEL bandwidth is much less than the separation of the few-cycle-driven frequency offset is of critical importance for the performance of the attosecond XFEL. In this case a crystal monochromator can be used to distinguish the attosecond pulses from the intense SASE pulses. Obviously, this requirement is easier to achieve for high power optical laser systems. For 800 nm laser radiation and for 0.1 nm output radiation, for example, the peak power of few-cycle laser pulse must be larger than 500-700 GW. This condition can be satisfied by a terawatt-scale sub-10 fs Ti:sapphire laser system which seems feasible.

Our scheme of attosecond X-ray source is based on the assumption that the beam density modulation does not appreciably change as the beam propagates through the energy modulator undulator. When the resonance condition takes place, the electrons with different arrival phases acquire different values of the energy increments (positive or negative), which result in the modulation of the longitudinal velocity of the electrons with the laser frequency. Since the velocity modulation is transformed into a density modulation of the electron beam when, passing the undulator, an additional wakefield exists because of a variation in amplitude density modulation. It is interesting to estimate the amount of bunching produced during the 800 nm undulator pass. An undulator is a sequence of bending magnets where particles with different energies have different path length, $`\mathrm{\Delta }z=R_{56}\delta E/E`$. The net compaction factor of the undulator is given by $`R_{56}=2\lambda _0N_\mathrm{w}`$, where $`\lambda _0=800`$ nm is the resonance wavelength and $`N_\mathrm{w}=2`$ is the number of undulator periods. An induced correlated energy spread at the exit of (800 nm) undulator is about 0.3%. Therefore, a rough estimate for the induced bunching is $`\delta a(\pi R_{56}/\lambda _0)(\delta E/E)3\times 10^2`$. Since this value is much less than unity, we can conclude that density modulation in the 800 nm undulator due to few-cycle-driven energy spread should not be a serious limitation in our case.

The next problem is that of synchronization. Frequency chirp in the XFEL is seeded by positioning a fs optical pulse on the electron bunch. Even when femtosecond pulses from laser system are synchronized to the photoinjector master clock with phase-locking technique, the synchronization of the optical seed laser with the electron pulses to an accuracy of 100 fs is not yet achievable. A more serious problem is the timing jitter of electron and seed laser pulses. The jitter of electron pulses originates in the photoinjector laser system (laser pulse jitter) and in the magnetic bunch compressors (from electron bunch energy jitter). Due to this uncertainty, not every fs optical pulse will produce an attosecond X-ray pulse. Random production of attosecond X-ray pulses needs to be controlled. A basic question at this point is how attosecond X-ray pulses will be identified. Separation of attosecond pulse frequency from the central frequency can be used to distinguish the 300 as pulses from the intense 100 fs pulses. Appearing of X-ray pulses at the frequency offset will indicate that the seed optical pulse is overlapped with the central part of the electron bunch.

Analysis of parameters of an attosecond X-ray source shows that its repetition rate is clearly limited by the value of repetition rate achievable with terawatt-scale sub-10 fs Ti:sapphire laser system having 1-kHz repetition rate . The single-pass scheme considered here is the simplest one. The laser beam, which is essentially unaltered in the electron beam modulation process, is then disposed of. This is not optimum for a couple of reasons. The idea is roughly the following. The attosecond pulse repetition rate can be significantly increased if the laser pulse can be reused, because the laser pulse suffers little loss in energy after each interaction with electron beam. The solution of this problem is a multipass approach (based on reflective optical elements) in which a laser pulse is made to pass through the modulator undulator a finite number of times before being thrown away. In this way, the attosecond pulse repetition rate is increased by increasing the number of laser pulses used.

## 5 Conclusion

Operation of the proposed scheme was illustrated for the parameters of the European XFEL. Although the present work is concerned primarily for use in the wavelength range around 0.1 nm, its applicability is not restricted to this range, for example 0.15 nm LCLS facility is a suitable candidate for application of attosecond techniques described here. It is important that proposed attosecond scheme is based on the nominal XFEL parameters, and operates in a ”parasitic” mode not interfering with the main mode of the XFEL operation. It can be realized with minimum additional efforts. The machine design should foresee the space for installation of modulator undulator and a viewport for input optical system. Many of the components of the required laser system can be achieved with technology which is currently being developed for applications other than the attosecond X-ray source. As a result, a laser system could be developed over the next few years and can meet the XFEL requirements well in advance of XFEL construction schedule.

## Acknowledgments

We thank G. Grübel, D. Novikov, E. Weckert for many useful discussions. We thank R. Brinkmann, J.R. Schneider, A. Schwarz, and D. Trines for interest in this work.


format: text

text:
# Classical phase-space descriptions of continuous-variable teleportation

## Abstract

The nonnegative Wigner function of all quantum states involved in teleportation of Gaussian states using the standard continuous-variable teleportation protocol means that there is a local realistic phase-space description of the process. This includes the coherent states teleported up to now in experiments. We extend the phase-space description to teleportation of non-Gaussian states using the standard protocol and conclude that teleportation of non-Gaussian pure states with fidelity of $`2/3`$ is a “gold standard” for this kind of teleportation.

Quantum teleportation is a process that can transfer an arbitrary quantum state from a system held by one party, usually called Alice, to a system held by a second party, usually called Bob. The process requires a pair of systems, shared by Alice and Bob, in an entangled state—the entangled resource—and a “small” amount of classical information transmitted from Alice to Bob. Originally proposed for qubit states bennett , teleportation protocols were later extended to states of a system described by continuous phase-space variables, such as a massive particle or a mode of an optical field vaidman ; braunsteinkimble . This continuous-variable teleportation protocol was implemented in an experiment that teleported a coherent state of an optical-frequency electromagnetic mode with fidelity $`0.58\pm 0.02`$ furasawa . Two recent experiments have improved the experimental fidelity of the teleported coherent state to values of $`0.64\pm 0.02`$ bowen  and $`0.61\pm 0.02`$ zhang .

In the standard continuous-variable teleportation protocol braunsteinkimble , Alice and Bob share an entangled Gaussian state of two modes, $`A`$ and $`B`$, which have annihilation operators $`a`$ and $`b`$; this entangled resource is ideally a two-mode squeezed state schumaker . The state to be teleported is the pure purestate  state $`\rho =|\psi \psi |`$ of a mode $`V`$ in Alice’s possession, which has annihilation operator $`v`$. The protocol consists of (i) Alice’s measuring two (commuting) joint quadrature components of modes $`V`$ and $`A`$, specifically the Hermitian real and imaginary parts of the operator $`v+a^{}`$, (ii) Alice’s communicating the (complex) result $`\xi `$ to Bob, and (iii) Bob’s displacing mode $`B`$ by $`\xi `$. The efficacy of the protocol is quantified by the fidelity between the output state of mode $`B`$ and the input state $`|\psi `$, averaged over the possible measurement results.

Experiments to date have teleported only coherent states. It is generally believed, though not proved, that teleporting coherent states with average fidelity $`>1/2`$ requires an entangled resource braunstein00 . It has thus been argued that teleportation of coherent states with fidelities above $`1/2`$ constitutes truly quantum teleportation braunstein00 ; braunstein01 . Using a variety of arguments, other workers have contended that $`=2/3`$ is the appropriate boundary between classical and quantum teleportation ralph ; grangier . While acknowledging the need for an entangled resource for teleporting coherent states with fidelity $`>1/2`$ (indeed, we provide additional evidence), we add a fresh perspective by investigating whether the entangled resource is used in a way that can be accounted for by a classical phase-space description. When such a description exists, it provides a local realistic hidden-variable model for the teleportation protocol.

Our investigation is motivated by the fact that the coherent states and the Gaussian entangled resource used in the experiments have nonnegative Wigner functions wigner , which are phase-space probability distributions that provide a classical description of measurements of the quadrature components. A nonnegative Wigner function does not give a classical description of measurements other than those of quadrature components; specific such measurements on an entangled Gaussian state cannot be given a local realistic description and thus violate Bell inequalities banaszek . Since the standard protocol uses only quadrature measurements, however, we conclude that for teleporting coherent states—or any Gaussian state—using the standard teleportation protocol, the nonnegative Wigner function of the three modes gives a classical, local realistic description for all fidelities. This means that all the experiments to date—and any such experiment that teleports coherent states, no matter what fidelity is achieved—can be accounted for in terms of purely classical correlations, with no need for a quantum-mechanical explanation.

To find situations where the Wigner function does not provide a classical phase-space description of the standard protocol, we must look to teleportation of non-Gaussian states, which (for pure states) have Wigner functions that take on negative values HP . To accommodate non-Gaussian states, we extend our hidden-variable model by allowing (i) Alice to substitute a randomly displaced state with nonnegative Wigner function in place of the non-Gaussian state and (ii) Alice and Bob to cheat by teleporting this new “smeared-out” state with perfect fidelity. Teleportation of non-Gaussian pure states with fidelity $`2/3`$ cannot be accommodated within this extended hidden-variable model, thus making a fidelity of 2/3 a “gold standard” for teleportation of non-Gaussian pure states.

We begin with a brief Wigner-function-based review of the teleportation protocol. The state $`\rho _{AB}`$ of modes $`A`$ and $`B`$ has Wigner function $`W_{AB}(\alpha ,\beta )`$, which is a quasidistribution for the c-number complex amplitudes $`\alpha `$ and $`\beta `$ corresponding to the annihilation operators $`a`$ and $`b`$ garzoll . In the standard protocol, $`W_{AB}(\alpha ,\beta )`$ is a Gaussian, but for the present, we allow it to be a general Wigner function. The (pure purestate ) state $`\rho =|\psi \psi |`$ of mode $`V`$ has Wigner function $`W_\rho (\nu )`$, where $`\nu `$ is the c-number complex amplitude corresponding to annihilation operator $`v`$. The overall Wigner function of the three modes is $`W_\rho (\nu )W_{AB}(\alpha ,\beta )`$. The state of mode $`B`$ after a measurement of $`v+a^{}`$ that yields result $`\xi `$ has Wigner function

$$W^{}(\beta |\xi )=\frac{1}{p(\xi )}d^2\nu d^2\alpha \delta (\nu +\alpha ^{}\xi )W_\rho (\nu )W_{AB}(\alpha ,\beta ),$$
(1)
where

$$p(\xi )=d^2\nu d^2\alpha d^2\beta \delta (\nu +\alpha ^{}\xi )W_\rho (\nu )W_{AB}(\alpha ,\beta )$$
(2)
is the probability to obtain result $`\xi `$.

Having received result $`\xi `$ from Alice, Bob displaces the complex amplitude of mode $`B`$ by $`\xi `$, yielding a state $`\rho _{\mathrm{out}}(\xi )`$ with Wigner function $`W_{\mathrm{out}}(\beta |\xi )=W^{}(\beta \xi |\xi )`$. The fidelity of this output state and the input state is $`F(\xi )=\psi |\rho _{\mathrm{out}}(\xi )|\psi `$. We are interested in the average of this fidelity over all measurement results, $`=d^2\xi p(\xi )F(\xi )=\psi |\overline{\rho }_{\mathrm{out}}|\psi `$, where $`\overline{\rho }_{\mathrm{out}}=d^2\xi p(\xi )\rho _{\mathrm{out}}(\xi )`$ is the average output state, having Wigner function

$$W_{\overline{\rho }_{\mathrm{out}}}(\beta )=d^2\xi p(\xi )W_{\mathrm{out}}(\beta |\xi )=d^2\nu G(\nu )W_\rho (\beta \nu ).$$
(3)
Here

$$G(\nu )=d^2\alpha d^2\beta \delta (\beta +\alpha ^{}\nu )W_{AB}(\alpha ,\beta )$$
(4)
is the (nonnegative) probability to obtain result $`\nu `$ in a measurement of $`b+a^{}`$ on modes $`A`$ and $`B`$. Equation (3) shows that the average output state is a mixture of displaced input states,

$$\overline{\rho }_{\mathrm{out}}=d^2\nu G(\nu )D(\nu )\rho D^{}(\nu ),$$
(5)
where $`D(\nu )`$ is the displacement operator.

We can now write the average output fidelity in two complementary forms,

$``$ $`=`$ $`{\displaystyle d^2\nu G(\nu )|C_\rho (\nu )|^2}`$ (6)
$`=`$ $`\pi {\displaystyle d^2\beta d^2\nu G(\beta \nu )W_\rho (\beta )W_\rho (\nu )},`$
where

$$C_\rho (\nu )=\psi |D(\nu )|\psi =d^2\mu W_\rho (\mu )e^{\nu \mu ^{}\nu ^{}\mu },$$
(7)
the symmetrically ordered (Wigner-Weyl) characteristic function of the input state cahillglauber , is the Fourier transform of the Wigner function. The first form in Eq. (6) comes directly from Eq. (5), and the second from from writing the fidelity as an overlap of the Wigner functions for the input and average output states. The effect of the initial state of modes $`A`$ and $`B`$ on the average fidelity is contained wholly in the marginal distribution $`G(\nu )`$. High-fidelity teleportation occurs when $`G(\nu )`$ is very narrow, i.e., when the quadrature components contained in $`b+a^{}`$ are sharp, expressing a particular kind of correlation between modes $`A`$ and $`B`$. Using the Fourier transform (7), we can derive two other, equivalent forms for the average fidelity,

$``$ $`=`$ $`\pi {\displaystyle d^2\beta d^2\nu \stackrel{~}{G}(\beta \nu )W_\rho (\beta )W_\rho (\nu )}`$ (8)
$`=`$ $`{\displaystyle d^2\nu \stackrel{~}{G}(\nu )|C_\rho (\nu )|^2},`$
where

$$\stackrel{~}{G}(\nu )=d^2\mu G(\mu )e^{\nu \mu ^{}\nu ^{}\mu }$$
(9)
is the Fourier transform of $`G(\mu )`$.

Before proceeding to the standard protocol and our hidden-variable models, we pause here to demonstrate the one technical result we need. We wish to find the maximum value of the integral

$$I=d^2\alpha d^2\beta e^{t|\alpha \beta |^2/2}W_{AB}(\alpha ,\beta ),t0,$$
(10)
over the Wigner function $`W_{AB}(\alpha ,\beta )`$ of a joint state $`\rho _{AB}`$ of modes $`A`$ and $`B`$. Introducing annihilation operators $`c=(a+b)/\sqrt{2}`$ and $`d=(ab)/\sqrt{2}`$, with corresponding c-number variables $`\gamma `$ and $`\delta `$, we can rewrite $`I`$ as

$$I=d^2\gamma d^2\delta e^{t|\delta |^2}W_{CD}(\gamma ,\delta )=d^2\delta e^{t|\delta |^2}W_D(\delta ),$$
(11)
where $`W_{CD}(\gamma ,\delta )=W_{AB}(\alpha ,\beta )`$ is the Wigner function written in terms of modes $`C`$ and $`D`$, and $`W_D(\delta )`$ is the Wigner function for mode $`D`$ alone. The integral can now be thought of as the expectation value, $`\mathrm{tr}(A_t\rho _D)=\mathrm{tr}(A_t\rho _{AB})`$, of the $`D`$-mode operator $`A_t`$ whose symmetrically ordered associated function cahillglauber  is $`e^{t|\delta |^2}`$, this operator being $`A_t=(1+t/2)^1[(1t/2)/(1+t/2)]^{d^{}d}`$.

The integral $`I`$ being the expectation value of $`A_t`$, $`I`$ is bounded above by the largest eigenvalue of $`A_t`$. Since $`A_t`$ is diagonal in the number-state basis, with eigenvalues that decrease in magnitude with the number of quanta, we have

$$I_{\mathrm{max}}=\text{(largest eigenvalue of }A_t\text{)}=\frac{1}{1+t/2},$$
(12)
with the maximum achieved if and only if $`\rho _{AB}`$ is the vacuum state for mode $`D`$.

For the case that $`\rho _{AB}`$ is a pure product state, $`|\mathrm{\Psi }=|\psi _A|\psi _B`$, which turns out to be the case of interest here, the condition for achieving the maximum becomes $`d|\mathrm{\Psi }=0`$ or, equivalently, $`a|\psi _A|\psi _B=|\psi _Ab|\psi _B`$, from which it follows that $`a|\psi _A=|\psi _A\psi _B|b|\psi _B`$ and $`b|\psi _B=|\psi _B\psi _A|a|\psi _A`$, implying that $`|\psi _A`$ and $`|\psi _B`$ are identical coherent states. Thus the only pure product states that achieve the maximum in Eq. (12) are products of identical coherent states.

We can use Eq. (12) to get one interesting result immediately: The maximum average fidelity for teleporting a coherent state using the standard protocol, but with a separable state for modes $`A`$ and $`B`$, is $`1/2`$. To show this, suppose first that modes $`A`$ and $`B`$ are in a pure product state, with factorizable Wigner function $`W_A(\alpha )W_B(\beta )`$. The characteristic function for any coherent state satisfies $`|C_{\mathrm{coh}}(\nu )|^2=e^{|\nu |^2}`$, so we can use Eq. (4) and the first form in Eq. (6) to write the average fidelity as

$$=d^2\alpha d^2\beta e^{|\alpha \beta |^2}W_A(\alpha ^{})W_B(\beta ).$$
(13)
Here $`W_A(\alpha ^{})`$ is the Wigner function for the time-reversed, parity-inverted state of mode $`A`$. The $`t=2`$ general bound (12) implies $`1/2`$, with equality if and only if mode $`A`$ is in a coherent state $`|\alpha `$ and mode $`B`$ is in the time-reversed, parity-inverted coherent state $`|\mathrm{-}\alpha ^{}`$.

Now suppose modes $`A`$ and $`B`$ are initially in a separable state, thus having a pure product-state ensemble decomposition. The fidelity is the average over the pure product-state ensemble, which shows that the fidelity is still bounded above by $`1/2`$, with equality if and only if the separable state is a mixture of product states of the form $`|\alpha |\mathrm{-}\alpha ^{}`$. This does not show that $`1/2`$ is the maximum fidelity for coherent-state teleportation in the absence of entanglement, since the result applies only to the standard protocol, but it is an additional piece of evidence, distinct from the results reported in Ref. braunstein01 .

We now take up again our analysis of the standard teleportation protocol, assuming that modes $`A`$ and $`B`$ are in a Gaussian state with Wigner function

$$W_{AB}(\alpha ,\beta )=\frac{4(𝚌^2𝚜^2)}{\pi ^2}e^{2𝚌(|\alpha |^2+|\beta |^2)2𝚜(\alpha \beta +\alpha ^{}\beta ^{})},$$
(14)
where $`𝚌`$ and $`𝚜`$ satisfy $`|𝚜|<𝚌\sqrt{1+𝚜^2}`$. This state is pure if and only if $`𝚌=\sqrt{1+𝚜^2}`$, in which case the state becomes a two-mode squeezed state with $`𝚌=\mathrm{cosh}2r`$ and $`𝚜=\mathrm{sinh}2r`$, where $`r`$ is the squeeze parameter schumaker . The state (14) is separable if and only if $`𝚌1|𝚜|`$ tutorial .

For the Wigner function (14), the distribution (4) is a Gaussian,

$$G(\nu )=\frac{𝚌+𝚜}{\pi }e^{(𝚌+𝚜)|\nu |^2}=\frac{2}{\pi t}e^{2|\nu |^2/t},$$
(15)
where $`t2/(𝚌+𝚜)`$ is the single parameter needed to characterize the fidelity that can be achieved with this entanglement resource. The Wigner function (3) of the average output state is the $`(s=t)`$-ordered quasidistribution, $`W_\rho ^{(s)}(\nu )`$ cahillglauber , of the input state:

$$W_{\overline{\rho }_{\mathrm{out}}}(\beta )=\frac{2}{\pi t}d^2\nu e^{2|\beta \nu |^2/t}W_\rho (\nu )=W_\rho ^{(s)}(\nu ).$$
(16)
For $`t=0`$, $`G(\nu )`$ is a $`\delta `$-function, and the output state is identical to the input state (perfect teleportation). For $`t=1`$, the Wigner function of the average output state is the Husimi $`Q`$ distribution of the input state, i.e., $`W_{\overline{\rho }_{\mathrm{out}}}(\beta )=W_\rho ^{(1)}(\beta )=Q_\rho (\beta )=\beta |\rho |\beta /\pi `$.

For $`0t<2`$ ($`𝚌>1𝚜`$), the state (14) is entangled, with the right sort of correlations for this protocol, these correlations decreasing as $`t`$ increases. At $`t=2`$, the state passes through the separability boundary $`𝚌+𝚜=1`$, and for $`t2`$, the state either is separable ($`𝚌1|𝚜|`$) or, though entangled ($`𝚌>1+𝚜`$), has the wrong sort of correlations for this protocol.

The average fidelity of Eqs. (6) and (8) now becomes

$`_\rho (t)`$ $`=`$ $`{\displaystyle \frac{2}{\pi t}}{\displaystyle d^2\nu e^{2|\nu |^2/t}|C_\rho (\nu )|^2}`$ (17)
$`=`$ $`{\displaystyle \frac{2}{t}}{\displaystyle d^2\beta d^2\nu e^{2|\beta \nu |^2/t}W_\rho (\beta )W_\rho (\nu )}`$
$`=`$ $`{\displaystyle d^2\beta d^2\nu e^{t|\beta \nu |^2/2}W_\rho (\beta )W_\rho (\nu )}`$
$`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle d^2\nu e^{t|\nu |^2/2}|C_\rho (\nu )|^2}.`$
The second form is the overlap of the Wigner function and the $`s`$-ordered quasidistribution for the input state. The first two derivatives of the last form show that $`_\rho (t)`$ is a strictly decreasing, strictly concave function of $`t`$. These forms also show that the average fidelity obeys the scaling relation $`_\rho (t)=2_\rho (4/t)/t`$, which again draws attention to the separability boundary at $`t=2`$.

Given that $`|C_{\mathrm{coh}}(\nu )|^2=e^{|\nu |^2}`$, the average fidelity for teleporting a coherent state is $`_{\mathrm{coh}}(t)=(1+t/2)^1`$. For number states $`|n`$, whose Wigner functions take on negative values (except for $`n=0`$), we have obtained an analytic formula for a generating function

$`(\lambda ,t)`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\lambda ^n_{|nn|}(t)`$
$`=`$ $`{\displaystyle \frac{1}{\sqrt{(1+t/2)^22\lambda (1+t^2/4)+\lambda ^2(1t/2)^2}}}.`$
The resulting fidelity for teleporting a number state is

$$_{|nn|}(t)=\frac{(1t/2)^n}{(1+t/2)^{n+1}}P_n\left(\frac{1+t^2/4}{1t^2/4}\right),$$
(19)
where $`P_n(x)`$ is a Legendre polynomial. This gives $`_{|nn|}(1)=2P_n(5/3)/3^{n+1}`$ and $`_{|nn|}(2)=(2n)!/2^{2n+1}(n!)^2`$. Another example of a state with a Wigner function that takes on negative values is the superposition $`|\psi =(|0+|1)/\sqrt{2}`$, which has teleportation fidelity $`_{|\psi \psi |}(t)=(1+3t/4+t^2/4)/(1+t/2)^3`$.

We now return to our Wigner-function-based discussion of hidden-variable models for teleportation. For Gaussian input states and for the two-mode entangled resource (14), all the Wigner functions are nonnegative wigner , so they provide a classical phase-space description—and hence a local hidden-variable description—for this kind of teleportation, no matter what fidelity is achieved. The hidden variables are the quadrature components of all the modes, and the overall Wigner function is a probability distribution for these hidden variables.

All non-Gaussian input pure states have Wigner functions that take on negative values (Hudson-Piquet theorem HP ) and thus cannot be incorporated in the simple hidden-variable model. To see what can be achieved within a classical phase-space description, suppose that before performing the teleportation protocol, Alice “kicks” the input state $`\rho `$ randomly in phase space. The random kick is described by a Gaussian so that the average state after the kick is

$$\rho ^{}=\frac{2}{\pi t}d^2\nu e^{2|\nu |^2/t}D(\nu )\rho D^{}(\nu ).$$
(20)
We choose the kicking strength $`t`$ to be the minimum value necessary to make $`\rho ^{}`$ have a nonnegative Wigner function, thus giving $`\rho ^{}`$ a classical phase-space description and allowing it to be incorporated it within our hidden-variable model. For all non-Gaussian states, this minimum kicking strength is one vacuum unit, i.e., $`t=1`$ lutkenhaus , implying that $`\rho ^{}`$ is the state whose Wigner function, $`W_\rho ^{}(\nu )=W_\rho ^{(1)}(\nu )=Q_\rho (\nu )`$, is the $`Q`$ function of the original state $`\rho `$. Further suppose that Alice and Bob cheat by teleporting $`\rho ^{}`$ with perfect fidelity. Then the fidelity of the overall process is the overlap of the Wigner and $`Q`$ functions of $`\rho `$, i.e., the $`t=1`$ fidelity (17) of the standard protocol. Notice that bigger kicks ($`t>1`$) would give smaller fidelity, making clear why we choose the smallest kicking strength consistent with giving $`\rho ^{}`$ a nonnegative Wigner function.

These considerations, coupled with wanting to know the maximum teleportation fidelity for a given entangled resource $`t`$, motivate us to find the maximum value of the average fidelity $`_\rho (t)`$ over all input pure states $`\rho `$, $`t=1`$ being the value relevant for our hidden-variable model. The task can be restated as finding the pure state $`\rho `$ that maximizes the overlap of the Wigner function and the $`s`$-ordered quasidistribution $`W_\rho ^{(s)}(\nu )`$. We apply the bound (12) to the third form of the average fidelity (17), in this case maximizing over pure product states $`\rho \rho `$. The resulting maximum is $`_{\mathrm{max}}(t)=(1+t/2)^1`$, achieved if and only if $`\rho `$ is a coherent state.

Each non-Gaussian input state $`\rho `$ has its own threshhold fidelity, $`_\rho (1)<_{\mathrm{max}}(1)=2/3`$, below which its teleportation can be accommodated within our extended phase-space hidden-variable model and above which it cannot. Thus teleportation of a non-Gaussian state with fidelity exceeding $`_\rho (1)`$ is required to rule out an explanation in terms of classical phase-space correlations. A fidelity of $`2/3`$ emerges as a “gold standard” for continuous-variable teleportation in the sense that teleportation of any non-Gaussian pure state with $`2/3`$ cannot be fitted within our extended hidden-variable model. This conclusion applies only to our phase-space-based hidden-variable model; we have not shown that there is no local hidden-variable model that can accommodate teleportation fidelities of 2/3 or above.

This work was supported in part by US Office of Naval Research Grant No. N00014-03-1-0426 (CMC) and by KBN Grant No. 2PO3B 02123 (KW).


text:
# Generalised zeta-function regularization for scalar one-loop effective action

Dipartimento di Fisica, Università di Trento
and Istituto Nazionale di Fisica Nucleare,
Gruppo Collegato di Trento, Italia

Abstract:

The one-loop effective action for a scalar field defined in the ultrastatic space-time where non standard logarithmic terms in the asymptotic heat-kernel expansion are present, is investigated by a generalisation of zeta-function regularisation. It is shown that additional divergences may appear at one-loop level. The one-loop renormalisability of the model is discussed and the one-loop renormalisation group equations are derived.

Within the so-called one-loop approximation in quantum field theory, the Euclidean one-loop effective action may be expressed in terms of the sum of the classical action and a contribution depending on a functional determinant of an elliptic differential operator, the so called fluctuation operator. The ultraviolet one-loop divergences which are eventually present, have to be regularised by means of suitable techniques (for recent reviews, see ). Very powerful is the use of zeta-function regularisation, which can be useful defined as a Mellin-like tranform of the heat-kernel trace of the small disurbance operator $`L`$. In standard cases, the small $`t`$ asymptotic expansion of the opearator $`Tr\mathrm{exp}(tL)`$ has only powers in $`t`$. In such a case the corresponding zeta function is regular at the origin and zeta-function regularisation give finite results. The situation dramatically changes when the heat-kernel asymptotics contains also logarithmic terms. In such a case the zeta-function has a pole at the origin and a revision of zeta-function regularisation is needed.

One may have logarithmic terms in the heat-kernel trace in the case of non smooth manifolds, for example when one considers the Laplace operator on higher dimensional cones , but also in 4-dimensional spacetimes with a 3-dimensional, non-compact, hyperbolic spatial section of finite volume . More recently the presence of logarithmic terms in self-interacting scalar field theory defined on manifolds with non-commutative coordinates have also been pointed out .

Here we shall discuss the modification of the zeta-function formalism due to the presence of logarithmic terms in the heat-kernel asymptotics.

Let us suppose to deal with a general expression of the kind

$`Tre^{tL}{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}B_jt^{j2}+{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}P_j\mathrm{ln}tt^{j2}.`$ (1)
The Mellin transform gives

$`\zeta (s|L)={\displaystyle \frac{1}{\mathrm{\Gamma }(s)}}{\displaystyle _0^{\mathrm{}}}𝑑tt^{s1}Tre^{tL}.`$ (2)
. Making use of the small $`t`$ asymptotics, one has

$`\zeta (s|L)={\displaystyle \frac{1}{\mathrm{\Gamma }(s)}}\left({\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{B_j(L)}{s+j2}}{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{P_j(L)}{(s+j2)^2}}+J(s)\right),`$ (3)
the function $`J(s)`$ being analytic. We see that in contrast with the standard situation, here the zeta function has also double poles and, if $`P_2`$ is non vanishing, it is no longer regular at the origin, but it has a simple pole with residue $`P_2`$. Another important consequence for physics is that, due to the presence of logarithmic terms, the heat-kernel coefficients $`B_n`$, with respect to scale transformations, transforms in a non homogeneous manner and in particular $`B_2`$ is not scale invariant.

The consequences of the presence of the pole at the origin on the one-loop effective action can be investigated by using the following regularisation for the functional determinant

$`\mathrm{ln}det\left(L/\mu ^2\right)_\epsilon ={\displaystyle _0^{\mathrm{}}}𝑑t{\displaystyle \frac{t^{\epsilon 1}}{\mathrm{\Gamma }(1+\epsilon )}}Tre^{tL/\mu ^2}={\displaystyle \frac{\zeta (\epsilon |L/\mu ^2)}{\epsilon }}={\displaystyle \frac{\omega (\epsilon |L/\mu ^2)}{\epsilon ^2}},`$ (4)
where we have conveniently introduced the new kind of zeta function $`\omega `$, regular at the origin, by means of the relation

$`\omega (s|L)=s\zeta (s|L),\omega (s|L/\mu ^2)=s\mu ^{2s}\zeta (s|L)=\mu ^{2s}\omega (s|L).`$ (5)
We may expand $`\omega `$ in Taylor series around $`\epsilon =0`$, obtaining in this way

$`\mathrm{ln}det\left(L/\mu ^2\right)_\epsilon ={\displaystyle \frac{1}{\epsilon ^2}}\omega (0|L/\mu ^2){\displaystyle \frac{1}{\epsilon }}\omega ^{}(0|L/\mu ^2){\displaystyle \frac{1}{2}}\omega ^{\prime \prime }(0|L/\mu ^2)+O(\epsilon ).`$ (6)
As a consequence, the one loop-divergences are governed by the two coefficients $`\omega (0|L/\mu ^2)`$and $`\omega ^{}(0|L/\mu ^2)`$, while the non trivial finite part is given by $`\frac{1}{2}\omega ^{\prime \prime }(0|L/\mu ^2)`$. This suggests a generalisation of the zeta-function regularisation for a functional determinant associated with an elliptic operator $`L`$, namely

$`\mathrm{ln}detL={\displaystyle \frac{1}{2}}\omega ^{\prime \prime }(0|L).`$ (7)
Of course, this reduces to the usual zeta-function regularisation when $`\zeta (s|L)`$ is regular at the origin.

The two coefficients governing the one-loop divergences can be computed making use of the meromorphic structure

$`\omega (s|L)={\displaystyle \frac{1}{\mathrm{\Gamma }(s)}}{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{sB_j}{s+j2}}{\displaystyle \frac{1}{\mathrm{\Gamma }(s)}}{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{sP_j}{(s+j2)^2}}+{\displaystyle \frac{s}{\mathrm{\Gamma }(s)}}J(s).`$ (8)
One has

$`\omega (s|L)=P_2(L)+[B_2(L)\gamma P_2(L)]s+O(s^2),`$ (9)
thus, using Eq. (5), one obtains

$`\omega (0|L/\mu ^2)=P_2(L),\omega ^{}(0|L/\mu ^2)=B_2(L)(\gamma +\mathrm{ln}\mu ^2)P_2(L),`$ (10)
$`\gamma `$ being the Euler-Mascheroni constant.

The model is one-loop renormalisable, if the dependence of $`B_2`$ and $`P_2`$ on the background field has the same algebric structure of the classical action and the divergences may be reabsorbed by the redefinition of mass and coupling constants.

The derivation of one-loop renormalisation group equations may be obtained by assuming that the mass and all coupling constants appearing in the classical action are depending on $`\mu `$ and requiring

$`\mu {\displaystyle \frac{d}{d\mu }}W_R=0.`$ (11)
In this way we get

$`W_R`$ $`=`$ $`S{\displaystyle \frac{1}{4}}\omega ^{\prime \prime }(0|L/\mu ^2)`$ (12)
$`=`$ $`S{\displaystyle \frac{1}{4}}\left[\omega ^{\prime \prime }(0|L)+2\mathrm{ln}\mu ^2\omega ^{}(0|L)+(\mathrm{ln}\mu ^2)^2\omega (0|L)\right].`$
Making use of Eqs. (10)-(12), we finally get (at one-loop level)

$`\mu {\displaystyle \frac{d}{d\mu }}S=\omega ^{}(0|L)+\mathrm{ln}\mu ^2\omega (0|L)=B_2(L)\left(\gamma +\mathrm{ln}\mu ^2\right)P_2(L).`$ (13)
If the theory is renormalisable, the action and the heat coefficients have the same structure in terms of the fields. More precisely, if the action has the form

$`S={\displaystyle 𝑑V\underset{\alpha }{}\lambda _\alpha (\mu )F_\alpha },`$ (14)
then

$`B_2(L)={\displaystyle 𝑑V\underset{\alpha }{}k_\alpha (\mu )F_\alpha },P_2(L)={\displaystyle 𝑑V\underset{\alpha }{}h_\alpha (\mu )F_\alpha },`$ (15)
where $`F_\alpha (1,\varphi ^2/2,\varphi ^4/24,\mathrm{})`$ are the independent building blocks, $`\lambda _\alpha (\mathrm{\Lambda },m^2,\lambda ,\mathrm{})`$ the collection of all coupling constants, including the ones concerning the gravitational action, while $`k_\alpha `$ and $`h_\alpha `$ are constants, which can be directly read off from the form of the heat coefficients.

From Eqs. (13)-(15) one obtains the differential equations for the beta functions in the form

$`\beta _\alpha \mu {\displaystyle \frac{d\lambda _\alpha }{d\mu }}=k_\alpha (\gamma +\mathrm{ln}\mu ^2)h_\alpha .`$ (16)
One recovers the usual result when all $`h_\alpha =0`$, namely $`P_2=0`$.

In this contributio we have presented new results concerning the one-loop effective action and the renormalisation group equations in non standard situations, in which the presence of logarithmic terms in the heat-kernel expansion give contributions which cannot be neglected. In particular, we have derived the contributions to beta fuctions and pointed out the explicit dependence on the renormalisation scale parameter.

Logarithmic terms in the heat-kernel expansion can be found, for example, in Friedmann-Robertson-Walker cosmology, in space-times with non trivial topology and also in non-local theories .


text:
# Thermodynamics and Phase Diagram of High Temperature Superconductors.

## Abstract

Thermodynamic quantities are derived for superconducting and pseudogap regimes by taking into account both amplitude and phase fluctuations of the pairing field. In the normal (pseudogap) state of the underdoped cuprates, two domains have to be distinguished: near the superconducting region, phase correlations are important up to the temperature $`T_\varphi `$. Above $`T_\varphi `$, the pseudogap region is only determined by amplitudes, and phases are uncorrelated. Our calculations show excellent quantitative agreement with specific heat and magnetic susceptibility experiments on cuprates. We find that the mean field temperature $`T_0`$ has a similar doping dependence as the pseudogap temperature $`T^{}`$, whereas the pseudogap energy scale is given by the average amplitude above $`T_c`$.

One of the most intriguing problems in high temperature superconductivity is the presence of a region above the critical temperature $`T_c`$ and below a temperature $`T^{}`$ where observable quantities deviate from Fermi liquid behaviour. This region is called pseudogap region loram94 ; timusk  because it contains effects similar to superconductivity like a partial suppression of electronic density of states.

The origin of such a pseudogap above $`T_c`$ is unclear. There are four major approaches concerning its theoretical understanding: the first is based on the formation of incoherent Cooper pairs above $`T_c`$. Phase order emery  or Bose condensation randeria  would then establish superconductivity at $`T_c`$. The second assumes that the pseudogap is induced by anti-ferromagnetic fluctuations schmalian . The third approach is based on spin-charge separation where spins bind together to form spin-singlets and the energy needed to split them apart leads to the formation of a ”spin-gap” lee . The fourth assumes the existence of a quantum critical point sachdev  but the latter has never been observed. However, these approaches seem to be unable to describe specific heat and magnetic susceptibility.

The main aim of this article is to show that various experimental observations can indeed be interpreted in terms of fluctuations of the pairing field $`\psi =|\psi |e^{i\varphi }`$, and that two temperature regions have to be distinguished (see Fig. 5): for a relatively small temperature interval $`T_c<T<T_\varphi `$ the phase of $`\psi `$ is still correlated in space over some correlation length $`\xi `$ (the Kosterlitz-Thouless correlation length in $`2D`$). Thus, in this regime, observables are governed by correlated phase fluctuations described by the XY-model. For $`T_\varphi <T<T^{}`$, phases of $`\psi `$ are essentially uncorrelated ($`\xi `$ is on the order of the lattice constant), but $`|\psi |`$ is still fluctuating and non-zero, signaling local pair fluctuations. This explains the wide hump seen in specific heat experiments loram94 , the depression of the spin susceptibility takigawa  and the persistence of the pseudogap for $`T<T^{}`$.

Our method has a major difference with respect to the Emery and Kivelson phase fluctuations scenario emery  of the pseudogap regime: we show that phase fluctuations influence the pseudogap only up to a temperature $`T_\varphi `$ which is much smaller than $`T^{}`$. Above $`T_\varphi `$, observables are thus only determined by the amplitude of the pairing field.

The picture of two different regimes above $`T_c`$ is also supported by other experiments: a) Demsar et al demsar , interpreting the real-time measurements of the quasiparticle relaxation dynamics, find a temperature interval of only a few K in which pair fluctuations associated with their collective phase are important, whereas the pseudogap persists to much higher temperatures. b) Hall effect measurements matthey  in underdoped GdBa2Cu3O7-x show a characteristic temperature $`T^{}`$ between $`T_c`$ and $`T^{}`$, at which the temperature dependence of $`\mathrm{cot}(\mathrm{\Theta }_H)`$ deviates from $`T^2`$ and the Hall coefficient has a peak. A possible explanation consists in considering vortex excitations as scattering centers modifying the Hall angle $`\mathrm{\Theta }_H`$ and Hall coefficient, which would again suggest identifying $`T^{}`$ with our $`T_\varphi `$ below which correlated vlatex ortices exist.

We base our calculations on a $`d`$-wave attractive Hubbard model

$$H=t\underset{i,j\sigma }{}c_{i\sigma }^{}c_{j\sigma }U\underset{i}{}Q_d^{}(i)Q_d(i)$$
(1)
with a hopping $`t`$ between nearest neighbour sites $`i`$ and $`j`$ on a square lattice. The interaction favours the formation of onsite $`d`$-wave pairs since $`Q_d^{}(i)=_jD_{ij}Q_{ij}^{}`$ where $`D_{ij}=1,(1)`$ for $`i`$ being the nearest neighbour site of $`j`$ in horizontal (vertical) direction. $`Q_{ij}^{}=(c_i^{}c_j^{}c_i^{}c_j^{})/\sqrt{2}`$ is the singlet pair operator of neighbouring sites. Decoupling the interaction with the help of a Stratonovich-Hubbard transformation, the partition function $`Z=\text{Tr}e^{\beta H}`$ is then

$$Z=Z_nD^2\psi 𝒯e^{_0^\beta 𝑑\tau _i\left(\frac{1}{U}|\psi |^2+\psi Q_d^{}(i)+\text{hc}\right)}_{H_n}$$

where $`\psi =\psi (i,\tau )=|\psi (i,\tau )|e^{i\varphi (i,\tau )}`$, and $`H_n`$ is the non-interacting part. The trace over the fermionic operators can be evaluated yielding

$$Z=D^2\psi e^{_0^\beta 𝑑\tau \left[_i\frac{1}{U}|\psi |^2+\text{Tr}\mathrm{ln}G\right]}.$$
(2)
Here $`G`$ is a Nambu matrix of one-electron Green functions for fermions interacting with a given, space and time dependent pairing field $`\psi (i,\tau )`$.

Expanding (2) in powers of $`\stackrel{}{}\psi `$, $`Z`$ can be written as a functional integral involving an action $`S[\psi ]`$ for a field $`\psi `$ that changes slowly in space and that can be taken time-independent:

$$S[\psi ]=d^dr\left[S_0(|\psi |)+S_1(\stackrel{}{}\psi )\right]$$
(3)
where $`S_0`$ is a local function of $`|\psi (𝐫)|`$, and $`S_1=c|\stackrel{}{}\psi |^2/2`$ where $`c`$ is a constant. $`d`$ is the dimension.

Now we would like to compute observables such as energy, specific heat and spin susceptibility. Our main strategy is to treat separately amplitude and phase fluctuations. In this spirit, two different approaches are possible: 1) The amplitude is fixed and determined by a suitable variational equation. 2) The energy is expanded around the average amplitude.

Variational Method: we can neglect amplitude correlations since simulations show that they are weak between different sites $`i,j`$: $`|\psi _i||\psi _j||\psi |^20`$ (i.e. the amplitude is always positive and shows no critical behaviour). Rewriting the free energy $`F`$ in terms of a constant amplitude $`|\overline{\psi }|`$ yields

$$F=\frac{1}{\beta }\mathrm{log}D\varphi e^{\beta \left[VS_0(|\overline{\psi }|)\mathrm{log}(|\overline{\psi }|)V/\beta +{\scriptscriptstyle d^drS_1}\right]}$$

where the Jacobian of the polar transformation is put into the exponential, and $`V`$ is the volume. Equating the derivative of $`F`$ with respect to $`|\overline{\psi }|`$ to zero leads to the self-consistent equation:

$$\frac{S_0(|\overline{\psi }|)}{|\overline{\psi }|}\frac{1}{\beta |\overline{\psi }|}+c|\overline{\psi }||\stackrel{}{}e^{i\varphi }|^2_{XY}=0.$$
(4)
The first term, called the amplitude contribution, leads to the BCS gap equation BCS  if other contributions are neglected. The second comes from the Jacobian and implies that the amplitude is never zero. The third term is the expectation value of the energy in the XY model depending on the constant coupling $`c|\overline{\psi }|^2/2`$. This contribution characterizes the influence of the phase fluctuations. Solutions of equation (4) are reliable for all temperatures except for $`T<<T_c`$. They are accurate at $`T_c`$ only in the underdoped regime.

Average Value Method: now, we keep the local coupling in $`S_1`$ between amplitude and phase by introducing the statistical Ginzburg-Landau (GL) model (see curty ). Averages like $`|\psi |`$ and $`S_1`$ are computed within the GL model with all fluctuations and correlations. Then fermionic observables are expanded with respect to the average amplitude.

Expanding the energy density obtained from action (3) around the average amplitude (see Fig. 1) yields:

$$E=S_S/VS_0(|\psi |)+S_1_S+𝒪(\delta |\psi |^2)$$
(5)
where the square of amplitude fluctuations is neglected in first approximation. Averages are computed using a normalised GL action:

$$S_{GL}[\psi ]=k_BV_0d^dr\left(U_{GL}+|\stackrel{}{}\psi |^2/2\right)$$
(6)
where $`U_{GL}=\eta ^2[(T/T_01)|\stackrel{~}{\psi }(𝐫)|^2+|\stackrel{~}{\psi }(𝐫)|^4/2]`$ comes from the expansion of $`S_0`$ in powers $`\beta |\psi |`$. $`T_0`$ is the mean field pairing temperature, and $`\stackrel{~}{\psi }=\psi /\left|\psi _0\right|`$ where $`\psi _0=\psi (T=0)`$. $`S_{GL}`$ is normalised with a lattice spacing $`\epsilon `$. $`\eta :=\epsilon /\xi _0`$, where $`\xi _0`$ is the mean field correlation length at zero temperature. $`V_0`$ is the zero temperature phase stiffness. The energy becomes

$$ES_0(|\psi |_{GL})+S_1_{GL}$$
(7)
where $`S_0(|\psi |)=(|\psi |^2/U)\frac{2}{\beta V}_q\mathrm{log}[2\mathrm{cosh}\beta E_q/2]`$ corresponds the BCS free energy for which the gap value is determined by the GL average. The $`d`$-wave quasi-particle energy is $`E_q=[(\epsilon _q\mu )^2+|\psi |^2\mathrm{cos}^2(2\theta )]^{1/2}`$ where $`\mu `$ is the chemical potential. and $`\theta `$ is the angle in $`k`$ space with respect to $`k_x`$ direction.

For computer simulations of the statistical ensemble $`\{\psi \}`$ under the action $`S_{GL}`$, we use a standard Monte Carlo procedure to update amplitude $`|\stackrel{~}{\psi }|`$ and a Wolff wolff  algorithm for the phase $`\varphi `$ as for the real $`\mathrm{\Phi }^4`$ model brower . Typically $`10^4`$ sweeps are needed to obtain good statistics for lattice size: $`N=40^2`$. The value of $`\eta `$ depends on the corse-graining procedure and is fixed to 3 in this method. Changes in $`\eta `$ modify slightly the shape of the average amplitude. The link between $`S_0`$ and the GL action is made by fixing ratio of $`|\psi _0|`$ and $`T_0`$. Since we are interested in the temperature domain from $`T_c`$ to $`T^{}`$, the ratio $`|\psi _0|/T_0`$ is fixed at $`T=T_c`$ by requiring that the specific heat jump recovers its BCS value in the overdoped regime. We found that this ratio must be $`|\psi _0|/T_0=3.06`$ in the $`s`$-wave case, and $`|\psi _0|/T_0=3.72`$ for $`d`$-wave symmetry.

Both approaches are valid below and above the critical temperature $`T_c`$ which is the temperature where the phase stiffness becomes zero.

The specific heat $`C`$ is the sum of $`C_0`$ and $`C_1`$, resp. amplitude and gradient contributions. Defining $`\gamma =C/(\gamma _nT)`$ where $`\gamma _n`$ the Sommerfeld constant, the reduced specific heat is

$$\gamma =\gamma _0\left(|\psi |\right)+\gamma _1$$
(8)
where $`\gamma _1`$ is divided by $`T_c`$ instead of $`T`$ since $`S_1`$ is classical and does not satisfy the third law of thermodynamics.

The amplitude contribution $`\gamma _0`$ is 1 at high temperature. The gradient contribution is normalised as:

$$\frac{C_1}{\gamma _n}=\frac{k_B}{\xi _0^3\gamma _n}\frac{C_1^{(s)}}{Nk_B}$$
(9)
where $`V=N\xi _0^3`$, and $`C_1^{(s)}/(Nk_B)`$ is the specific heat per number of lattice sites coming from the simulations. Experiments give $`\gamma _n`$ 26 mJ K-1 mol-1 = 252 J K-1 m-3. For the fit of Fig. 2, using the reasonable value $`\xi _016`$ Å, we get the dimensionless constant $`\alpha =k_B/(\xi _0^3\gamma _n)13.5`$. In Fig. 2 the experimental specific heat of YBa2Cu3O6.73 loram94  is fitted using the variational method reproducing the double peak structure: a sharp peak at $`T_c`$ coming mainly from correlated phase fluctuations and a wide hump rounded by amplitude fluctuations (see Fig. 1 for typical distributions). The crossover temperature $`T_1`$ corresponds to the temperature where $`\gamma _1`$ is less than approximately 2% of the normal specific heat. In the amplitude equation (4), a 2-dimensional density of states $`D(\epsilon )=1/W`$ is used with $`W=5000\text{K}`$, $`\mu =0.25W`$ and $`U=959\text{K}`$. These values give $`T_0200\text{K}`$ and $`\psi _02.14T_0`$ in agreement with experiments kugler . The other parameters are $`V_0=72\text{K}`$ and $`\eta =5`$.

The average value method is used to reproduce specific heat measured for different dopings in Fig. 3. For underdoped systems $`x<0.80`$ we use simulations in $`D=2`$. For the more overdoped, $`x0.80`$, simulations are done in $`D=3`$. Parameters $`V_0`$ and $`T_0`$ extracted from the best fits are shown in the phase diagram of Fig. 5. $`\alpha `$ ranges from 7 to 15.

The magnetic susceptibility $`\chi `$ has two contributions: the paramagnetic spin susceptibility $`\chi _p`$ and the orbital diamagnetic susceptibility $`\chi _d`$. $`\chi _p`$ has been measured by Takigawa et al takigawa  on powder YBa2Cu3O6.63 using Cu and O NMR experiments. The direct contribution coming from phases is negligible because NMR probes essentially the presence of pairs which is related to the amplitude. Therefore, $`\chi _p`$ is given by the amplitude contribution:

$$\chi _p=\frac{\chi _0}{2T}_0^{2\pi }\frac{d\theta }{2\pi }_0^{\mathrm{}}𝑑\epsilon \mathrm{cosh}^2\frac{\sqrt{\epsilon ^2+|\psi |_{GL}^2\mathrm{cos}^2(2\theta )}}{2T}$$
(10)
$`\chi _0`$ is the Pauli spin susceptibility. In Fig. 4, we compare the result of equation (10) and Takigawa’s measurements on powder YBa2Cu3O6.63. The best fit yields $`T_0=159.2\text{K}`$, $`V_0=59.3\text{K}`$. The competition between amplitude and thermal energy enters into the spin susceptibility (10) by the ratio $`|\psi |/T`$. The temperature $`T_\varphi `$ where phases start to influence the thermodynamics is the temperature where $`|\psi |`$ deviates from the average amplitude computed for random phases. This temperature is found above $`T_c`$ at $`T_\varphi 90\text{K}`$. The orbital diamagnetic susceptibility $`\chi _d`$ in YBa2Cu3O6.60 shows fluctuations effects up to about 15 K above $`T_c`$. For $`T>T_\varphi `$, phase fluctuations are so strong that $`\chi _d`$ vanishes sewer .

The Green function can be calculated in CPA (Coherent Potential Approximation) for electrons that are scattered from spatially uncorrelated ”pair impurities”. We extended the $`s`$-wave CPA gyorffy  by using a $`d`$-wave amplitude and uncorrelated phases. We have computed the density of states at half-filling $`\mu =0.5W`$ and maximum amplitude $`|\psi |=0.8W`$. The result is a pseudogap of width $`\mathrm{\Delta }0.2W`$ although phases are completely uncorrelated.
Phase diagram: Values of $`T_0`$, $`V_0`$ and $`T_\varphi `$ extracted from specific heat fits (Fig. 3) are reported in Fig. 5. $`T_1`$ is the temperature where $`C_1/\gamma _n=2\%`$. $`T_\varphi `$ is computed in the same way as shown in Fig. 4 using average amplitudes from specific heat fits. Here, it is the temperature where the amplitude differs by $`2\%`$ from the amplitude in a random phase field. The energy scale $`E_g:=|\psi |_{T=200\text{K}}`$ of the pseudogap is defined as the amplitude at $`T=200\text{K}`$. $`E_g`$ shows the same doping dependence as the one found by Loram and Tallon loram2 . However it is not related to any hidden critical point. Phase correlations above $`T_c`$ grow rapidely in the underdoped regime following the $`T_\varphi `$ line, an have a similar doping dependence as Nernst effect results wang . The gradient specific heat from $`S_1`$ disappears more rapidely like in the Hall effect matthey .
Discussion: amplitude and phase fluctuations are the key for understanding the pseudogap regime of underdoped high temperature superconductors. Phase correlations disappear completely near a temperature $`T_\varphi `$ above $`T_c`$, and therefore, for $`T>T_\varphi `$, the pseudogap region is dominated by amplitude fluctuations. The mean field temperature $`T_0`$ has a similar doping dependence as $`T^{}`$, signaling that the pseudogap region is due to independent fluctuating pairs. Comparisons with measured specific heat on underdoped YBCO reproduce the double peak structure: a sharp peak at $`T_c`$ coming mainly from phase correlations and a separate wide hump below $`T^{}`$ rounded by fluctuations. The spin susceptibility, related to the amplitude, recovers its normal behaviour near $`T^{}`$ whereas the orbital magnetic susceptibility, related to phases, disappears near $`T_\varphi `$. These considerations are independent from the underlying pairing mechanism, and any microscopic theory inducing pairing should lead to similar conclusions.
All these findings provide additional evidence for the fact that superconductivity and pseudogap have the same origin. The former is primarily related to phases of the pairing field, which are ordered below the transition temperature and whose correlations survive over a limited temperature region above $`T_c`$. The pseudogap regime of underdoped materials then extends to much higher temperatures due to the persisting amplitude fluctuations of the pairing field.
We thank J. Loram, B. Janko, C. Moca and C. Wirth. This work has been supported by the Swiss National Science Foundation.


text:
# Modelling Correlations in Portfolio Credit Risk

(January 19, 2004)

## Abstract

The risk of a credit portfolio depends crucially on correlations between the probability of default (PD) in different economic sectors. Often, PD correlations have to be estimated from relatively short time series of default rates, and the resulting estimation error hinders the detection of a signal. We present statistical evidence that PD correlations are well described by a (one-)factorial model. We suggest a method of parameter estimation which avoids in a controlled way the underestimation of correlation risk. Empirical evidence is presented that, in the framework of the CreditRisk+ model with integrated correlations, this method leads to an increased reliability of the economic capital estimate.

Managing portfolio credit risk in a bank requires a sound and stable estimation of the loss distribution with a special emphasis on the high quantiles denoted as Credit Value-at-Risk (CreditVaR). The difference between the CreditVaR and the expected loss has to be covered by the economic capital, a scarce resource of each bank. From a risk management perspective, the definition of industry sectors allows to diversify credit risk. The degree to which this diversification is successful depends on the strength of correlations between the sectors. Moreover, the correlations between sector PDs crucially influence the CreditVaR and hence the economic capital.

In large banks, the concentration risk in industry sectors is a key risk driver. Recently, several approaches for describing and modelling concentration risk were discussed Bue ; NagBa01 ; Frey+01 . In CreditRisk+ CSFB , concentration risk is modelled as a multiplicative random effect on the PD per counterpart in a given sector. In the original version of CreditRisk+, the loss distribution is calculated for independent sector variables. Correlations between PD fluctuations in different sectors can be integrated into CreditRisk+ with the method of Bürgisser et al. Bue . For the calculation of the CreditVaR it is important whether input parameters like the correlation coefficients between sector PDs are known or must be estimated. In the latter case, this estimation leads to an additional variability of the target estimate, in our case the portfolio loss. In this way, uncertainty in the estimation of PD correlations translates itself into uncertainty of the economic capital of a bank.

The estimation of cross–correlations is difficult due to the ”curse of dimensionality”: if the length $`T`$ of the available time series is comparable to the number $`K`$ of industry sectors, the number of estimated correlation coefficients is of the same order as the number of input parameters with the result of large estimation errors. A way out of this dilemma is the use of a factor model with a reduced dimensionality of the parameter space. We present evidence that the PD correlations for $`K=20`$ industry sectors are well captured by a one–factor model. Surprisingly, even the parameter estimation for the one–factor model is subject to large statistical fluctuations and gives rise to a considerable uncertainty in the CreditVaR. We discuss these fluctuations in detail and suggest a bootstrap method which allows to find an upper limit for the parameters. We assess the impact of different conservative estimates with respect to the CreditVaR of a realistic portfolio.

Description of data set
As the economic activity and the probability of default in a given industry sector is not directly observable, we approximate it by the insolvency rate in that sector over the last $`T`$ years. The probability of insolvency $`PD_{kt}`$ of sector $`k`$ in year $`t`$ is calculated as the ratio of the number of insolvencies in that sector to the total number of companies in the sector

$$\widehat{PD}_{kt}=\frac{_{A\mathrm{sector}\mathrm{k}\mathrm{in}\mathrm{year}\mathrm{t}}I_{\{A\mathrm{fails}\}}}{_{A\mathrm{sector}\mathrm{k}\mathrm{in}\mathrm{year}\mathrm{t}}}.$$
(1)
With the help of insolvency rates, the default probability for a given company $`A`$ can be factorized into an individual expected PD $`p_A`$ and the sector specific relative PD movement $`X_k`$ with expectation $`X_k=1`$ according to

$$P(A\mathrm{fails})=\mathrm{p}_\mathrm{A}\mathrm{X}_\mathrm{k}\mathrm{with}\mathrm{X}_{\mathrm{kt}}=\frac{\widehat{\mathrm{PD}}_{\mathrm{kt}}}{\frac{1}{\mathrm{T}}_\mathrm{t}\widehat{\mathrm{PD}}_{\mathrm{kt}}}.$$
(2)
For this study, we use sector specific default histories as supplied by the federal statistical office of Germany. We analyze default rates for a segmentation of the economy into 20 sectors and estimate the sample covariance matrix $`\mathrm{\Sigma }^{\mathrm{emp}}`$ and sample correlation matrix Cemp as

$$\mathrm{\Sigma }_{ij}^{\mathrm{emp}}=\frac{1}{T1}\underset{t=1}{\overset{T}{}}(X_{it}1)(X_{jt}1),C_{ij}^{\mathrm{emp}}=\mathrm{\Sigma }_{ij}^{\mathrm{emp}}/\sigma _{X_i}\sigma _{X_j}$$
(3)
with $`\sigma _{X_i}^2=\mathrm{\Sigma }_{ii}^{\mathrm{emp}}`$.

Test for independent sectors
We first ask whether the sample correlation matrix of the PD time series is compatible with the hypothesis of zero correlations. Ideas for testing this hypothesis for covariance matrices date back to the seventies JOH , and were recently generalized to situations where the number of time series is larger than the sample size LED . Here, we use an adaption of the tests JOH ; LED  to test for the equivalence of correlation matrix to the unit matrix. The test statistics

$$R=\frac{1}{K}\mathrm{tr}\left[𝖢^2\right]1,$$
(4)
for a correlation matrix C is both $`K`$– and $`T`$–consistent with the $`T`$–limiting distribution $`(T1)KR/2\stackrel{D}{}\chi _{K(K1)/2}^2`$ Rosenow03 . The prefactor $`T1`$ rather than $`T`$ is chosen to improve the finite $`T`$ properties of the test. For our example with $`T=7`$ and $`K=20`$, we find $`R=5.805`$, whereas the critical value for $`\alpha =0.05`$ is $`R_{\mathrm{crit}}=3.719`$. Hence, the independence of sector PDs must not be assumed and a model describing sector correlations is needed.

Description of one–factor model
We diagonalize the empirical cross correlation matrix Cemp and rank order its eigenvalues $`\lambda _{i,\mathrm{emp}}<\lambda _{i+1,\mathrm{emp}}`$. As we are interested in modelling correlations rather than covariances, we normalize the $`X_{it}`$ such that they have the same, namely the average variance $`\sigma _X^2=(1/K)_{i=1}^K\sigma _{X_i}^2`$ and subtract the mean

$$\stackrel{~}{X}_{it}=(X_{it}1)\frac{\sigma _X}{\sigma _{X_i}}.$$
(5)
We use the components of the eigenvector $`𝐮_{\mathrm{emp}}^{(K)}`$ corresponding to the largest eigenvalue $`\lambda _{K,\mathrm{emp}}=10.38`$ to define a factor time series

$$Y_t=\underset{i=1}{\overset{K}{}}u_{i,\mathrm{emp}}^{(K)}\stackrel{~}{X}_{it}.$$
(6)
As compared to averaging the sector variables without prior information, the definition of the factor time series from the eigenvector with largest eigenvalue makes sure that the factor explains a maximum amount of correlations. The idea arises from factor analysis, see e.g. Mar . In the context of stock returns, a time series defined according to the prescription of Eq. 6 was found to agree well with a value weighted stock index Ple . We expect that the factor time series Eq. 6 describes economy wide changes of relative PD, possibly weighted by the economic relevance of the individual sectors.

We model the correlations between relative PD movements by a one–factor model

$$\stackrel{~}{X}_{it}=b_iY_t+ϵ_{it}.$$
(7)
The coefficients $`\{b_i\}`$ are found by performing a linear regression. To see whether a one–factor model fully describes the correlations between the $`\{\stackrel{~}{X}_{it}\}`$, we apply the test Eq. 4 to the correlation matrix of the residuals $`\{ϵ_{it}\}`$. Taking into account that the regression reduces the effective length of the residual time series from $`T`$ to $`T1`$, we find $`R=4.409`$ slightly below the threshold $`R_{\mathrm{crit}}=4.463`$. As the assumption of uncorrelated residuals is not rejected, no further factors are needed for the description of correlations.

The point estimator can now be calculated under the assumption that the residua $`\{ϵ_{i,t}\}`$ are iid observations from uncorrelated random variables $`ϵ_i`$ $`i=1,\mathrm{},K`$, i.e. $`ϵ_iϵ_j=0`$ for $`ij`$. Defining the factor variance $`\sigma _Y^2=\frac{1}{T1}_{t=1}^TY_t^2`$, one finds the point estimator for the cross correlation matrix as

$$C_{ij}^{\mathrm{point}}=\delta _{ij}+(1\delta _{ij})b_ib_j\sigma _Y^2/\sigma _X^2.$$
(8)
The largest eigenvalue of Cpoint is found to be $`\lambda _{K,\mathrm{point}}=10.66`$ in good agreement with the original largest eigenvalue. In addition, the corresponding eigenvector $`𝐮_{\mathrm{point}}^{(K)}`$ is found to be very close to the original eigenvector (Fig.1).

Fluctuations in empirical correlation matrices – a toy model

In this section, we use the results of Monte Carlo simulations to study the relation between the cross correlation matrix Cmodel resulting from infinitely long model time series and matrices Csim numerically calculated from finite time series of length $`T`$. We find that the $`\{𝖢^{\mathrm{sim}}\}`$ differ from Cmodel both in a systematic way, for example a shift of the largest eigenvalue towards larger values, and a random way, i.e. an individual member of the simulated ensemble deviates significantly from the average Laloux+99 ; Plerou+99 .

To simplify simulations, we rewrite Eq. 7 as

$$\stackrel{~}{X}_{it}=\alpha \beta _iF_t+c_i\eta _{it}.$$
(9)
The random variables are rescaled according to $`c_i\eta _{it}=ϵ_{it}`$ and $`\alpha \beta _iF_t=b_iY_t`$ such that their variances are $`\mathrm{var}(F)=\sigma _X^2`$ and $`\mathrm{cov}(\eta _i\eta _j)=\sigma _X^2\delta _{ij}`$. In addition, the $`\{\beta _i\}`$ obey the normalization condition $`_{i=1}^K\beta _i^2=1`$, which makes them comparable to the components of the eigenvector $`𝐮^{(K)}`$. The model parameters are subject to the constraint $`c_i^2=1\alpha ^2\beta _i^2`$ in order to enforce $`\mathrm{var}(\stackrel{~}{X}_{it})=\sigma _X^2`$.

The model is completely defined by i) the parameter $`\alpha `$ determining the largest eigenvalue, ii) the parameters $`\{\beta _i\}`$ and iii) by the probability distribution of the random variables $`F`$ and $`\{\eta _i\}`$. As the empirical PD movements are commonly assumed to follow a gamma distribution, we model the random variables $`F`$ and $`\eta _i`$ by gamma distributions as well gammadistribution . In addition, we use normal distributions for the random variables and find that the deviation from a simulation with gamma distributed variables is smaller than $`3\%`$. As the simulation of gaussian random variables is computationally much more efficient than the simulation of gamma distributed variables, we use normally distributed variables for the computationally demanding selfconsistent calculations described in the next section.

The infinite time series correlation matrix of the model is given by

$$C_{ij}^{\mathrm{model}}=\delta _{ij}+(1\delta _{ij})\alpha ^2\beta _i\beta _j.$$
(10)
In this section, we study the outcome of model simulations for the particularly simple hypothetical case $`u_{i,\mathrm{model}}^{(K)}=\beta _i1/\sqrt{K}`$ in order to gain qualitative insight into the occurring fluctuations. For a given value of $`\alpha `$ corresponding to $`\lambda _{K,\mathrm{model}}=10.38`$, we perform 10000 Monte Carlo simulations of Eq. 9. We compute the pdf for the largest eigenvalue $`\lambda _{K,\mathrm{sim}}`$ of Csim and the corresponding eigenvector $`𝐮_{\mathrm{sim}}^{(K)}`$. We find that both quantities have broad distributions (see Fig. 2). The distribution of eigenvalues has an average $`\lambda _{K,\mathrm{sim}}=10.72`$, which is significantly larger than the true eigenvalue $`\lambda _{K,\mathrm{model}}=10.38`$. In addition, one finds simulated eigenvalues as low as $`\lambda _{K,\mathrm{sim}}=5`$. We quantify the systematic shift of eigenvalues by the average $`\mathrm{\Delta }\lambda =\lambda _{K,\mathrm{sim}}\lambda _{K,\mathrm{model}}`$. The magnitude of eigenvalue fluctuations is described by the standard deviation

$$\sigma _\lambda =\sqrt{\lambda _{K,\mathrm{sim}}^2\lambda _{K,\mathrm{sim}}^2}.$$
(11)
For the distribution shown in Fig. 2 we find $`\sigma _\lambda =2.42`$.

There are significant fluctuations of eigenvector components as well. For theoretical eigenvector components $`u_{i\mathrm{model}}^{(K)}=0.224i`$ one even finds negative empirical components indicating spurious anticorrelations, which would lead to dangerous hedges in credit portfolios. Specifically, we calculate the standard deviation

$$\sigma _{u_i}=\sqrt{(u_{i,\mathrm{sim}}^{(K)})^2u_{i,\mathrm{sim}}^{(K)}^2}$$
(12)
and find $`\sigma _{u_i}=0.083`$. Since the $`u_{i,\mathrm{model}}^{(K)}`$ do not vary across $`i`$ we only need to estimate one $`\sigma _{u_i}`$.

As a conclusion, even if the generating process for relative PD movements is a simple one–factor model, the empirically found parameters can deviate significantly from the theoretical ones. We advocate the point of view that the empirical Cemp has to be viewed as a member of such a fluctuating ensemble in that its eigenvalues and eigenvectors can deviate significantly from the unknown “true” correlation matrix of PD movements Laloux+99 ; Plerou+99 . Then, the statistical properties of the ensemble $`\{𝖢^{\mathrm{sim}}\}`$ can be used to derive error bars for both the largest eigenvalue and the components of the corresponding eigenvector.

Conservative estimates
How can we use these results to make a reliable estimate for the correlation matrix of relative PD movements? A bank needs to act in a conservative manner to prevent insolvency. Using the empirical correlation matrix, the bank risks that the correlations are ”accidently” low. The most conservative approach would be to assume all correlations to be $`1`$, i.e. $`u_i^{(K)}=\frac{1}{\sqrt{K}}i`$. But now the model would effectively be a one-sector model. Any possibility to measure concentration risk in certain industry sectors would be prevented. The model would not encourage diversifying the business across sectors.

As a controlled mediation we introduce ”cases” of add-ons of $`x=1,2,3`$ standard deviations to the fluctuating quantities such that the predicted risk for a portfolio is increased. This means correcting the eigenvalue towards larger values and the eigenvector components towards the value $`u_i^{(K)}1/\sqrt{K}`$ indicating the same correlation strength for all sectors and the absence of hedge possibilities.

Specifically, we let $`u_{i,\mathrm{case}}^{(K)}=1/\sqrt{K}`$ if

$$|u_{i,\mathrm{emp}}^{(K)}1/\sqrt{K}|<x\sigma _{u_i}$$
(13)
and $`u_{i,\mathrm{case}}^{(K)}=u_{i,\mathrm{emp}}^{(K)}\pm x\sigma _{u_i}`$ otherwise. The sign is chosen such that the overall risk increases, i.e. such that $`u_{i,\mathrm{case}}^{(K)}`$ falls between the empirical value and $`1/\sqrt{K}`$. After applying these corrections, the eigenvector is normalized. For this calculation, we fix the parameter $`\alpha `$ such that the simulated largest eigenvalue $`\lambda _{K,\mathrm{sim}}`$ is equal to the empirically observed one. We calculate the $`\{\sigma _{u_i}\}`$ selfconsistently, i.e. we calculate $`\sigma _{u_i}`$ for the $`u_{i,\mathrm{case}}^{(K)}`$ which solves Eq.13. The results are shown in Fig.3. We see that for increasing $`x=1,2,3`$, the model eigenvector comes closer to the null hypothesis of an eigenvector with identical components. While the empirical eigenvector has significant negative components indicating anticorrelations between some of the sectors, the negative components in $`𝐮_{1\sigma }^{(K)}`$ are already strongly reduced and completely gone in $`𝐮_{2\sigma }^{(K)}`$.

Similarly, we add a fluctuation margin to the model eigenvalue eigenvalueshift  such that

$$\lambda _{K,\mathrm{case}}=\lambda _{K,\mathrm{emp}}+x\sigma _\lambda .$$
(14)
Here, $`x`$ specifies the width of the confidence interval for the estimation of $`\lambda _{K,\mathrm{model}}`$. We perform this calculation selfconsistently, i.e. we calculate $`\mathrm{\Delta }\lambda `$ and $`\sigma _\lambda `$ for the $`\lambda _{K,\mathrm{case}}`$ which solves Eq.14. We find $`\lambda _{K,1\sigma }=11.17`$, $`\lambda _{K,2\sigma }=12.93`$, and $`\lambda _{K,3\sigma }=15.42`$.

Economic implications of the different correlation matrices
So far, we have described five different estimates for the cross correlation matrix, i.e. Cemp, Cpoint, C$`{}_{1\sigma }{}^{}{}_{}{}^{\mathrm{model}}`$, C$`{}_{2\sigma }{}^{}{}_{}{}^{\mathrm{model}}`$, and C$`{}_{3\sigma }{}^{}{}_{}{}^{\mathrm{model}}`$. To judge the economic implications of these estimates, we study the differences in the loss distribution resulting from these correlation estimations. The CreditVaR is a key quantity in banking when it comes to risk management. Reduced by the expected loss, it quantifies the capital needed to prevent insolvency for a given level of security. As capital is a resource it must be considered in the pricing of credit and trading products. Therefor, we quantify the impact of the different correlation estimates by calculating their influence on CreditVaR.

The portfolio we study is realistic – although fictitious – for an international bank. It consists of 4934 risk units distributed asymmetrically over 20 sectors with 20 to 500 counterparts per sector. The total exposure is 70 bn Euro with a largest exposure of 1.5 bn Euro and a smallest exposure of 0.25 mn Euro. The counterpart specific default probability varies between 0.03% and 7%, the expected loss for the total portfolio is 373.3 mn Euro. Table I shows the CreditVaR calculated by using CreditRisk+ and the method of Bürgisser et al. Bue  for integrating correlations.

We note that the use of a one–factor model changes the CreditVaR only by two percent as compared to the sample cross correlation matrix. Thus, the assumption of a one–factor description and the increase of estimation confidence achieved with this assumption yields portfolio risk estimation compatible with the parameter free estimation.

Our aim is to estimate a quantile of a probability distribution – namely the CreditVaR of the portfolio loss distribution. In the presence of an unknown parameter, it is a well established statistical result (see LEH ) that the use of the point estimate for the parameter – derived by a model or not – leads to an underestimation of the quantile estimate. To account for this additional estimation insecurity, we add a volatility $`\sigma `$ to the parameter estimate, i.e. the correlation matrix. When applying a one–$`\sigma `$ estimate, the CreditVaR increases by 400 mn Euro, for the two–$`\sigma `$ estimate there is another increase by 300 mn Euro, and using the three–$`\sigma `$ estimate the CreditVaR increases by yet another 200 mn Euro. To put these numbers in perspective, we note that the CreditVaR without including correlations is found to be 2.27 bn Euro, and that the assumption of full correlations among all sectors leads to a CreditVaR of 3.952 bn Euro. Because negative PD correlations are not plausible from an economic point of view, the use of the two–$`\sigma `$ estimate guarantees a sufficient forecast reliability on the one hand and allows for some guidance for economical decision on the other hand.

In summary, we have shown that correlations between empirical default rates for economic sectors are statistically significant and must be taken into account. We have described these correlations with a one–factor model and found that this description reproduces well the empirical correlations. However, when using the model to generate short time series and calculating their correlation matrix, one typically observes large statistical fluctuations in the correlation structure. Due to these fluctuations, the parameter estimation for a one–factor model is plagued by large uncertainties. When estimating the model parameters in such a way that the empirically observed ones appear as a worst case scenario, the reliability of the estimate is increased in a systematic way, leading to a moderately increased CreditVaR.

Acknowledgement: We would like to thank A. Müller-Groeling for initiating this project. We thank S. Lösch, C. von Lieres, and A. Wilch for useful discussions.