Fengping Li
email: fengping.li@strath.ac.uk

Wenlong He

Adrian W Cross

Craig R Donaldson

Liang Zhang

Alan D R Phelps

Kevin Ronald

Design and Simulation of a ~390 GHz 7 th Harmonic Gyrotron using a large orbit electron beam

Keywords: Harmonic gyrotron, terahertz, large orbit electron beam

A ~390 GHz harmonic gyrotron based on a cusp electron gun has been designed and numerically modelled. The gyrotron operates at the 7 th harmonic of the electron cyclotron frequency with the beam interacting with a TE 71 waveguide mode. Theoretical as well as numerical simulation results using the 3D Particle-In-Cell (PIC) code MAGIC are presented. The cusp gun generated an axis-encircling, annular shaped electron beam of energy 40 keV, current 1.5 A with a velocity ratio of 3. Smooth cylindrical waveguides have been studied as the interaction cavities and their cavity Q optimized for 390 GHz operation. In the simulations ~600 W of output power at the design frequency has been demonstrated.

Introduction:

Microwave radiation sources are one of the crucial technologies for modern society. High power submillimeter wave sources have attracted significant interest recently due to their potential applications in electron cyclotron resonance heating [START_REF] Wagner | The New Multi-Frequency Electron Cyclotron Resonance Heating System for ASDEX Upgrade[END_REF], plasma diagnostics [START_REF] Saito | Development of a Sub Terahertz High Power Pulse Gyrotron for Collective Thomson Scattering[END_REF], space radar systems [START_REF] Bhartia | Millimeter Wave Engineering and Applications[END_REF], communication systems [START_REF] Siegel | Terahertz Technology[END_REF], materials science and medical imaging and spectroscopy [START_REF] Clery | Brainstorming their way to an imaging revolution[END_REF][START_REF] Torrezen | CW Operation of a Tunable 330/460 GHz Gyrotron for Enhanced Nuclear Magnetic Resonance[END_REF]. The gyrotron is well known as a coherent millimeter wave source capable of high power and high frequency output which is extendable to the Terahertz (THz) frequency region. Recently a gyrotron achieved kilowatts of output power operated at the fundamental electron cyclotron frequency of 1 THz [START_REF] Yu | Generation of 1.5-kW 1-Thz Coherent radiation from a gyrotron with a pulsed magnetic field[END_REF] by using a pulsed solenoid. However, continuous wave (CW) operation at such a frequency is a formidable task because of the large magnetic field (~40 T) that is required. As the output frequency increases, both larger magnetic fields and reduced interaction region size are needed when operating at the fundamental cyclotron frequency. An alternative approach to generate CW high frequency radiation is to work at higher cyclotron harmonics [START_REF] Idehara | A high harmonic Gyrotron with an axis-encircling electron beam and a permanent magnet[END_REF][START_REF] Wang | Stable 1 MW, third-harmonic Gyro-TWT Amplifier[END_REF][START_REF] Wang | Design of a W-Band Second-Harmonic TE 02 Gyro-TWT Amplifier[END_REF][START_REF] Agusu | Design of a 400 GHz Gyrotron for DNP-NMR Spectroscopy[END_REF] which allow the use of larger cavity sizes and smaller magnetic fields by a factor of s, where s is the harmonic number. A CW high harmonic gyrotron operating at a frequency of ~390 GHz at relatively lower B-field was designed and is presented in this paper. In this gyrotron cyclotron resonance takes place between the seventh harmonic of the electron cyclotron frequency (s = 7) and the TE 71 waveguide mode. A large orbit electron beam from a cusp gun is used to reduce the possibility of parasitic interactions. The small-signal theory of the beam-wave interaction was used to calculate the growth rate and the starting current. A smooth-bore cavity was designed with a suitable Q value by optimizing the angle of the output taper and the length of the cavity so that the starting current requirement is met. Finally, the beam-wave interaction of the gyrotron was simulated using the 3D Particle-In-Cell (PIC) code MAGIC and the results are presented.

Cusp electron gun

While operating at harmonics, undesired mode competition as well as parasitic oscillations will be more problematic as the harmonic number increases, therefore, it is advantageous for an electron gun to be able to generate a mode-selective beam. A large orbit electron beam from a cusp gun was selected because for such a beam its harmonic of the cyclotron mode (s) only interacts with a TE mn waveguide mode (m and n are azimuthal and radial mode index respectively) when s = m [START_REF] Mcdermott | Small-signal theory of a large-orbit electron-cyclotron harmonic maser[END_REF]. A cusp electron gun operates by an electron beam passing through a non-adiabatic magnetic field reversal from negative to positive sign.

Invariance of the magnetic moment results in a large orbit, annular electron beam [START_REF] Schmidt | Nonadiabatic Particle Motion in Axialsymmetric Fields[END_REF][START_REF] He | A W-band Gyro-BWO based on a Helically Corrugated Waveguide[END_REF][START_REF] He | Axis-encircling electron beam generation using a smooth magnetic cusp for gyro-devices[END_REF]. Several successful terahertz range large orbit harmonic gyrotrons experiments have been carried out at resonant magnetic fields 10.5 -14 T [START_REF] Bratman | Submillimeter-wave large-orbit Gyrotron[END_REF][START_REF] Bratman | Large-orbit Gyrotron Operation in the Terahertz Frequency Range[END_REF]. The 3D PIC code MAGIC was used to simulate and optimize the design of the cusp electron beam source. The optimized geometry of the cusp gun including magnetic field profile and the resultant beam trajectories are shown in figure 1. An insertion shows in detail the B-field reversal in the cusp region. The B-field profile of the cusp was produced by two coils, a main coil and a reverse coil. The main coil produces a uniform magnetic field at the downstream region which should have a certain length of the flattop region required by the beam-wave interaction in the cavity. Two extra layers of coils were added to both ends of the main solenoid to sharpen its B-field profile so that the overall turns and length of the solenoid were reduced. These extra coils were carefully designed to avoid any bumps in the magnetic field profile while minimizing the rise and fall distance of the magnetic field. The design parameters were optimized resulting in the production of a 40 kV, 1.5 A large orbit electron beam in the simulations. This electron beam has been optimized through adjustment of the magnetic coil arrangement to operate at a velocity ratio from 1.2-3 [START_REF] He | Design, Simulation and Experiment of a Cusp Electron Beam for Millimeter Wave Gyrodevices[END_REF][START_REF] Donaldson | Design and Numerical Optimization of a Cusp-Gun-Based Electron Beam for Millimeter-Wave Gyro-Devices[END_REF] with an inner radius of 0.27 mm and an outer radius of 0.36 mm. This is consistent with the analytical calculated Larmor radius of 0.31 mm. The reflected electrons in figure 1 are low energy electrons emitted during the rise-time of the applied voltage pulse. A velocity pitch spread of ~25% was simulated when reaches 3. Further increase of the beam would result in a significant increase of beam spread in hence the reflection of some beam electrons.

Gyrotron design

Gyrotrons are based on the CRM (Cyclotron Resonance Maser) [START_REF] Chu | The electron cyclotron maser[END_REF] instability and operate through a fast cyclotron wave interaction, which takes place near the frequency cut-off region of the waveguide mode. In a high harmonic gyrotron, the beam (harmonic number s > 1) interacts with a high order TE mn waveguide mode where m = s.

All the possible interactions are shown in figure 2. When the 7 th harmonic of the electron cyclotron frequency is in resonance with the TE 71 mode, interactions below this mode are eliminated because the corresponding harmonic beam lines are far below their corresponding TE m1 waveguide mode curves except for the case of the TE 61 mode and s = 6. The other possible parasitic interactions might be from higher modes, i.e., TE 81 with s = 8 and TE 91 with s = 9, etc. However the most dangerous one is with the neighboring TE 81 interaction as the starting current increases as the mode index increases. But the TE 81 interaction has a starting current higher than the beam current. In this way any other parasitic interactions of higher modes are automatically eliminated. As the gyrotron works near the cut-off frequency region, and a large orbit gyrotron requires s = m, the operating frequency (390 GHz) and harmonic number (s = 7) as well as the azimuthal index of the waveguide mode (m = 7) were chosen prior to the cavity size. The interaction of the waveguide mode and axisencircling electron beam cyclotron mode can be described by the equation below [START_REF] Chu | Gain and Bandwidth of the gyro-TWT and CARM Amplifiers[END_REF][START_REF] Furuno | Operation of a Large-Orbit High-Harmonic Gyro-Traveling-Wave Tube Amplifier[END_REF]:

2 2 2 2 2 2 2 4 2 4 ( )( ) ( ) ( )( ) ( ) (
)

z mn z z mn m mn L w A w m mn k c x m k v s x c J x r r I I r J x β ω ω ⊥ - = - - - -Ω ′ (1)
Where k z is the propagation constant; x mn is the nth root of ( ) 0 m J x ′ = . r w is the radius of the cylindrical cavity;

N e is the axial linear density of electrons. is the electron gyrating frequency, where B is the magnetic field. I is the beam current and 3 0 0 [START_REF] Siegel | Terahertz Technology[END_REF] 17

A I m c e πε = ≈
kA is the Alfven current. The growth rate of this interaction can be obtained by solving the imaginary part of in equation ( 1) as a function of axial wavenumber k z and is shown in figure 3. In the diagram the TE 71 cylindrical waveguide mode is interacting with the 7 th harmonic of the electron cyclotron frequency of a beam with a velocity ratio of 3 at a magnetic field of 2.09 T. The radius of the cylindrical cavity is 1.055 mm, and the cavity length is 50 mm at a reasonable cavity Q of 2000. The cavity length L and the output taper angle are optimized for a reasonably high Q which is a very important parameter in any gyrotron design. The total Q of the resonant cavity is given by 1 1 1

d Q Q Q Ω = + (2)
The diffraction quality factor Q d and the ohmic quality factor Q which are caused by diffraction and ohmic losses respectively, can be calculated using the equations below:

2 1 2 ( ) 4 1 d L Q R R λ π = - (3) 2 2 2 2 2 2 2 1 ( )[ ( )] mn mn w w mn x x c Q r r x m δ ω Ω = + - (4)
Where

2 δ ω µσ =
is the skin depth of the cavity wall, is the conductivity of the cavity material, µ is the absolute permeability. R 1 and R 2 are the reflection coefficients at both ends of the cavity interaction region. The cavity Q was calculated using the code Cascade [START_REF] Ives | Cascade-An advanced computational tool for waveguide components and window design[END_REF]. In figure 4(a) the curve shows the simulated cavity Q as a function of the output taper angle at 390 GHz. In this diagram, the cavity was a 40 mm long smooth cylindrical cavity with a radius of r w = 1.055 mm. The output taper was 10 mm in length and the input section of the cut-off taper had a fixed radius of 0.8 mm and a length of 20 mm. Though in a real situation the tolerance of the cavity dimensions, roughness of the cavity surface [START_REF] Thorpe | RF Conductivity in Copper at 8 mm Wavelengths[END_REF] and the purity of the materials will affect the value of the ohmic losses at high frequencies, the ideal surfaces and oxygen free high conductivity copper were applied to these calculations. From the diagram we can see the cavity Q increases as the angle increases, but in order to sustain enough output power and to decrease the heat load of the cavity [START_REF] Zapevalov | Low-Q cavities for high-power gyrotrons[END_REF], the taper angle needs to be properly selected. Considering all of these factors a 6 degree output taper was chosen in the cavity design. The cutoff taper angle was also studied by changing the length of the taper with a fixed cut-off input radius 0.8 mm. Cascade simulations show that these changes do not affect the cavity Q value very much (less than 20 when the cut-off taper length was varied from 5 mm to 20 mm). shows the simulated cavity Q as a function of cavity length using the code Cascade. In the simulation the input end had a 20 mm taper with an input radius of 0.8 mm. The cavity radius remains the same as before which was 1.055 mm and at the end it had a 6 degree output taper. Simulations using CST Microwave studio showed that there is no significant mode conversion (<10% of power) for this output taper at the operating mode. However, further increase of the output taper length or output taper angle will cause significant mode conversion. For example, when the taper length was increased to 20 mm, 70% of the power was converted to the TM 52 mode. The starting oscillation currents for this cavity at the design mode and competing modes can be calculated by using the equation shown below [START_REF] Mcdermott | Small-signal theory of a large-orbit electron-cyclotron harmonic maser[END_REF]:

2 2 2 2 2 2 1 ( )( )( )(1 ) 16
( ) ( ) ( ) ( / ) w z z A s mn m mn m mn w Lr k I m I Q s c x J x R X J x c r γ ω β ω β - ⊥ = - Ω × ′ Ω (5)
Where:

2 2 2 2 2 ( / ) ( ) {2 2 ( )( ) [( ) 1]( ) ( / ) 2 ( / ) ( )( } ( ) ( / ) ( )( )[ ( ) ] ( ) w n mn w z z mn n mn w mn n mn w w n mn w z z z r s J x c r k c R X x s cx J x c r x c J x c r G X r J x c r ck X G X s ck β β β β β β β β β ω ω β ω ⊥ ⊥ ⊥ ⊥ ⊥ ⊥ ⊥ ⊥ Ω Ω = - + - ′ Ω Ω ′′ Ω + + ′ Ω Ω ′ - Ω (6)
and:

( ) z z X s k c ω β = -Ω , 2 2 2 2 [cos( 2) / (1 )] ( ) ( ) [sin( 2) / (1 )] ( ) qX X odd q G X qX X even q π π - = - (7)
Here The minimum starting current for the 7 th harmonic mode is 1.4 A, at a magnetic field of 2.13 T which are both within the capability of our cusp gun design current of 1.5-1.7 A and the solenoid magnetic field of 2.16 T. This diagram shows clearly that the neighboring modes s = 8 and s = 6 will be eliminated as the starting current for the 8th harmonic requires a starting current of more than 6 A, and the magnetic field requirement for the 6th harmonic is above 2.17 T.

Simulation results

In order to demonstrate the beam-wave interaction can take place in the designed gyrotron cavity, the code MAGIC was used to simulate the harmonic gyrotron interaction. The R-Z plane geometries of the smooth cylindrical structures and electron trajectories are shown in figure . 6. An ideal beam with the parameters from the cusp electron simulation but without velocity and radius spreads was emitted from the left end of the cavity and propagated in a uniform B-field produced by a solenoid in the magic simulation. The right end of the solenoid sits in the middle of the output taper in the Z axis. During the simulation, the cavity length was reduced from 50 mm to 30 mm as the cavity ohmic loss was neglected in the simulation. A schematic diagram of the simulation is shown in figure 6 with the magnetic field profile produced by the solenoid. The simulation result is shown in figure 7. In the simulation the large orbit annular electron beam had a voltage of 40 kV, current of 1.5 A with a beam velocity ratio of 3. The magnetic B-field applied in the cavity was 2.09 T. The magnetic field difference between the MAGIC and analytical calculation is due to the grid density in the gyrotron simulation, when a higher grid density is used, the B-field required in the simulation will be larger and closer to the analytically calculated value. ns arises due to the cavity of 30 mm length being sufficiently long in the absence of wall loss, for the interaction to enter phase trapped saturation resulting in re-absorption of the wave energy -these can be seen in figure 7(c) where the electron distribution in phase angle illustrates the trapping occurring. When 0.1 mm beam thickness and 25% beam spread were introduced in the simulation, it was found that the total output power dropped to ~500 W while a competing TE 61 became slightly dominant at the output. The output power at the designed frequency and mode dropped to ~200 W.

Conclusion

In conclusion a ~390 GHz harmonic gyrotron based on a large orbit electron beam was designed and simulated. The cavity has been analyzed and optimized. Theoretical and numerical simulation results demonstrate that the beam interacts with the single TE 71 mode and a simulated output power of ~600 W was achieved while using a 40 kV, 1.5 A, large orbit electron beam with velocity ratio of 3. It was found from simulation that spread in will cause the reduction of the output power and start-up of mode competition. The effect of beam parameters such as velocity spread, envelop ripple, might have on the starting current and mode excitation [START_REF] Bratman | Theory of gyro devices with thin electron beams (large-orbit gyrotron)[END_REF] will be studied both numerically and experimentally in the future.

Figure 1 .

 1 Figure 1. Simulated cusp gun and resultant beam trajectories with the detailed magnetic field configuration at the cusp point.

Figure 2 .

 2 Figure 2. Dispersion of the TE m1 waveguide modes (m = 1-8) and beam-wave resonances at different harmonics.

  is the relativistic Lorentz factor; m 0 is the rest mass of an electron, V is electron beam voltage.v cβ ⊥ ⊥ = ; L r v ⊥= Ω is the electron Larmor radius.

Figure 3 .

 3 Figure 3. The growth rate at the 7 th harmonic interaction (upper diagram) and the schematic gyrotron geometry (lower diagram).

Figure 4 .

 4 Figure 4. (a) Cavity Q as a function of output taper angle and (b) Cavity Q as a function of cavity length

Figure 4 (

 4 Figure 4(b)shows the simulated cavity Q as a function of cavity length using the code Cascade. In the simulation the input end had a 20 mm taper with an input radius of 0.8 mm. The cavity radius remains the same as before which was 1.055 mm and at the end it had a 6 degree output taper. Simulations using CST Microwave studio showed that there is no significant mode conversion (<10% of power) for this output taper

  figure.5. In this calculation the cavity Q = 2000, and the electron beam velocity ratio was 3.

Figure 5 .

 5 Figure 5. Starting currents for the 7 th harmonic TE 71 mode and its neighboring TE m1 modes.

Figure 6 .

 6 Figure 6. A Schematic diagram of gyrotron in the R-Z plane with simulated electron beam trajectories.

Figure 7 .

 7 Simulation results of the 7 th harmonic gyrotron with a smooth cavity: (a) output power, (b) output spectrum and (c) electron phase angle as a function of axial position (d) microwave E field pattern at resonance. The electron phase angle as a function of axial position Z (figure 7(c)) indicates the interaction occurred at the 7 th harmonic of the electron cyclotron frequency. Very pure radiation output at a frequency of 384 GHz was observed in the simulation as shown in figure 7(b) and ~600 W of output power was achieved in this simulation. The reduction in power observed after ~30

Acknowledgments

The authors would like to thank the Faraday Partnership in High Power RF and the Scottish Universities Physics Alliance (SUPA) and Overseas Research Students Awards Scheme (ORSAS) for supporting this project. The authors would also like to thank "Calabazas Creek Research" for providing the 1 month free license for Cascade.


Christian Pitra

Susana Suárez-Seoane

Carlos A Martín

Wolf-Jürgen Streich

Juan C Alonso

C Gortázar

W.-J Streich

Linking habitat quality with genetic diversity: a lesson from great bustards in Spain

Keywords: Otis tarda, Geographic information systems, Habitat suitability index, mtDNA

The effects of habitat loss and fragmentation on the genetic structure and variability of wild populations have received wide empirical support and theoretical formalization. By contrast, the effects of habitat quality seem largely underinvestigated, partly due to technical difficulties in properly assessing habitat quality. In this study, we combine geographic information system (GIS)based habitat-quality modelling with a landscape genetics approach based on mitochondrial DNA markers to evaluate the possible influence of habitat quality on the levels and distribution of genetic diversity in a range of natural populations (n = 15) of Otis tarda throughout Spain. Ninety-three percent of the population represented by our countrywide sample lives in good-quality habitats, while 4.5% and 2.5% occur respectively in intermediate and poor habitats. Habitat quality was highly correlated with patch size, population size and population density, indicating the reliability and predictive power of the habitat suitability model. Genetic diversity was significantly correlated with habitat quality, size and density of the population, but not with patch size. Three of a total of 20 existing matrilineages from the species' current genetic pool are restricted to poorquality habitats. This study therefore highlights the importance of considering both population genetics and habitat quality in a species of high conservation priority.

Introduction

The great bustard, Otis tarda, is one of the most characteristic avian species of the lowland dry grassland ecosystems in Europe. Currently, the species is threatened and listed as vulnerable A2c+3c+4c version 3.1 under the current International Union for Conservation of Nature (IUCN) Red List of Threatened Species due to their severely fragmented habitat and continuing decline in extent and quality of habitat (www.iucnredlist.org). A study based on nuclear and mitochondrial (mt)DNA analyses showed that in Europe great bustard populations are subdivided into two different evolutionary significant units (ESUs; Moritz 1994) respectively in the European mainland (including Ukraine and south Russia) and Spain (Pitra et al. 2000). The Iberian peninsula is home to ~27,500-30,000 great bustards, approximately 60% of the world population of this species (Palacín and Alonso 2008). For great bustards in Spain, metapopulation structure has recently been shown based on ecological surveys [START_REF] Alonso | The great bustard Otis tarda in Spain: present status, recent trends and an evaluation of earlier censuses[END_REF][START_REF] Alonso | Status and recent trends of the great bustard (Otis tarda) population in the Iberian peninsula[END_REF][START_REF] Alonso | Distribution dynamics of a great bustard metapopulation throughout a decade: influence of conspecific attraction and recruitment[END_REF]) and genetic analyses (Martín et al. 2002;[START_REF] Alonso | Genetic diversity of the great bustard in Iberia and Morocco: risks from current population fragmentation[END_REF]). The birds occupy patches of open, steppe-like landscapes with cerealfallow rotations, a habitat that is particularly under threat because of agriculture intensification and land abandonment linked to the European agricultural policy. Other threats are urbanization, infrastructure expansion and, locally, reforestation. Parallel to an ongoing population decline during the 20th century throughout Europe, the Spanish population of the great bustard has fragmented in several breeding areas distributed across the regions of Madrid, Castilla-La Mancha, Extremadura, Castilla y León, Navarra, Aragón and Andalusia [START_REF] Alonso | Status and recent trends of the great bustard (Otis tarda) population in the Iberian peninsula[END_REF]). Populations at the periphery (the last three mentioned previously) were the most genetically differentiated from the main ones in central Spain [START_REF] Alonso | Genetic diversity of the great bustard in Iberia and Morocco: risks from current population fragmentation[END_REF]). Accordingly, peripheral populations are also found to be more susceptible to decline and extinction than populations at the core of their range (Lane and Alonso 2001; [START_REF] Alonso | The great bustard Otis tarda in Andalucía, southern Spain: status, distribution and trends[END_REF]. These are reasons why the great bustard population in Spain, although now thought to be stable, is apparently concentrating at high-quality areas, disappearing from poor-quality ones [START_REF] Alonso | Distribution dynamics of a great bustard metapopulation throughout a decade: influence of conspecific attraction and recruitment[END_REF]; see also Pinto et al. 2005 for the Portuguese population).

In this context, it becomes particularly important to quantify how variation in ecological factors linked to habitat quality might affect and shape observed levels of genetic variability in wild populations. The effects of habitat loss and fragmentation on the genetic structure and variability of wild populations have received wide empirical support and theoretical formalization [START_REF] Frankham | Inbreeding and conservation: a threshold effect[END_REF]. By contrast, the effects of habitat quality (or suitability) seem largely underinvestigated, partly due to technical difficulties in properly assessing habitat quality (Vergeer et al. 2003). There are a large variety of methods in the literature to measure it (see Johnson 2007 for a review), either assessing directly habitat attributes or recording variables from individuals/populations. Recent developments in geographic information systems (GIS) and habitat modelling provide an ideal methodological framework in which high-resolution analyses of habitat attributes can be performed at the landscape level. Models allow gaining insights into the environmental variables (such as biotic, physical, climatic or human-induced habitat conditions) acting on species survival by comparing current species distribution with the optimal values for each variable (Guisan and Zimmermann 2000;Osborne et al. 2001;Suárez-Seoane et al. 2002). The output of the models is a habitat suitability index (HSI) which ranges from 0 (low quality) to 1 (high quality). It has been widely used since the early 1980s to determine wildlife habitat quality (U.S. Fish andWildlife Service 1980, 1981), as [START_REF] García | Spatial analysis of the habitat quality in a fragmented[END_REF], among other authors, did for a fragmented population of little bustard in NW Spain. This study examines the interactions between habitat quality and genetic diversity in a range of natural populations (n=15) of O. tarda throughout the Spain and considers the implications of the results for the species' conservation. Specifically, we used a GIS-based distribution model previously developed for O. tarda in Spain to infer a spatially explicit HSI capable of identifying, categorizing and ranking landscape quality (Suárez-Seoane et al. 2002). The spatial pattern of the species' genetic diversity was then overlaid to address the following objectives: (1) to investigate at a nationwide scale the spatial patterns of genetic variation, (2) to evaluate the possible influence of habitat quality and some demographic variables (patch size, population size, population density) on genetic variation and (3) to identify those haplotypes and determine the proportion of the overall genetic variation of the species in the study area currently living under poor habitat conditions and deserving high conservation priority.

Material and methods

Bird data and sample collection

Data about numbers and distribution of the species were compiled from censuses carried out in Spain by different teams of observers between mid March and early April (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008), when great bustards gather at their mating areas and thus are easier to count (details and references in the study of [START_REF] Alonso | Status and recent trends of the great bustard (Otis tarda) population in the Iberian peninsula[END_REF]Palacín and Alonso 2008). A team consisted of two observers with experience in counting bustards and knowledge of the area to be surveyed. They drove four-wheel-drive vehicles slowly through most tracks, trying to count all birds and to avoid duplicated counts. Each census was carried out from dawn to dusk, with a midday pause, when birds often sit down and are harder to detect. Sampling regions were chosen to represent the distribution range of the species in Spain [START_REF] Alonso | Status and recent trends of the great bustard (Otis tarda) population in the Iberian peninsula[END_REF]. The corresponding population and patch sizes were delimited according to discontinuities in the distribution of both, birds and areas of steppe-like habitat, considering the movements of a large sample of radio-tagged individuals (more than 700 birds in total, 1991-2009, unpublished data). A total of 302 tissue samples of blood, feathers or unfertilised eggs from various Spanish regions were collected from breeding adults or chicks throughout the great bustard breeding range in Spain. The UTM coordinates were determined for each sample using a GARMIN-12 GPS.

Genetic analysis

DNA extraction, polymerase chain reaction (PCR) amplification of a 657 base pair (bp) fragment of the mitochondrial control region (Domain I=283 bp, Domain II= 374 bp), and DNA sequencing were performed using the protocols outlined by Martín et al. (2002). The PCR products were purified (QIAquick PCR purification kit/ Qiagen), sequenced by the extension-dideoxy-chain termination method (Sanger et al. 1977) with a commercial kit (Dye Terminator Cycle Sequencing Kit, Applied Biosystems/Perkin Elmer) and analysed on a 373 A sequencer (Applied Biosystems). Sequences were aligned manually and all variable sites were confirmed by visual inspection of the chromatograms. Sequences of the same haplotype were identified using the computer program Collapse 1.1 (Posada and Crandall 1998). The resulting 20 different haplotypes were submitted to the GenBank database and included in the subsequent analysis (Table 1). Initial sequence comparisons and measures of genetic variability (haplotype diversity, h; nucleotide diversity, π; and number of nucleotide differences, k) were performed using DNASP 3.0 (Rozas and Rozas 1997). Graphically, haplotype areas were visualized by 95% confidence ellipses for the coordinates of the respective samples. The SPSS 16.0 (SPSS Inc., Chicago, IL, USA) and SYSTAT 12.0 (SYSTAT Software Inc., Richmond, CA, USA) statistical software packages were used for the statistical calculations. The significance level was set to α=0.05.

Habitat modelling

The distribution of great bustards in Spain was modelled at 1×1 km resolution from a set of predictor variables including vegetation, topography and human disturbance estimates. Vegetation was characterized by using a 12month time-series of normalized difference vegetation indices (NDVI), calculated from NOAA-AVHRR satellite data. NDVI separates green vegetation from other surfaces because chlorophyll absorbs red light and reflects near infrared wavelengths. It may vary depending on land use, season and climate. Since regional differences in the timing of seasons and agricultural production can happen across Spain, we applied standardized principal components analysis (PCA) to replace the original 12-monthly variables with 12 independent components that may be interpreted as particular environmental features or events. Topographic variability was estimated in a 5×5 pixel array of 200-m pixels, where altitude was computed to 5-m vertical resolution. Human disturbances were measured as the proportion of 200-m pixels in a 5×5 array containing roads and towns and the distance in km to the nearest 200-m pixel containing those features. Also, proportion and distance to rivers were included in the analysis. Relationships between occurrence bird data and environmental predictors were explored by applying generalized additive models (GAMs). We used S-plus 2000 (Venables and Ripley 1999) and the GRASP (Generalized Regression Analysis and Spatial Prediction) interface (Lehmann et al. 2003) to fit cubic splines with four degrees of freedom for each predictor, using a logit link and binomial error structure. Models were generated by means of a backwards variable selection at a significance level of 0.05 (Pearce and Ferrier 2000) testing, for each selected variable, significant differences between smooth and linear terms. Terms contribution to the model was assessed using a likelihood ratio test and model performance was evaluated through 10-fold crossvalidation. Output model was a HSI ranging from 0 to 1 (Fig. 1), which can be used as an indicator of habitat quality for the species. The model fitted the data well, with an ROC score of 0.96, which reduced by only 0.01 when it was cross-validated, suggesting its robustness and predictability. HSI was produced in IDRISI 32.11 [START_REF] Eastman | Idrisi for Windows[END_REF]. For a more detailed explanation of the data sets and statistical methods used, see Suárez-Seoane et al. (2002). We then calculated the averaged HSI within a moving window of 3×3 km (i.e. 9 km 2 ; the minimum size of a lek) centered sequentially on each pixel corresponding to each sample location. We assigned these HSI mean values to the respective samples. Differences in habitat quality among haplotypes were tested by performing ANOVA.

Results

Spatial pattern in genetic diversity

We found 20 distinct mtDNA haplotypes, defined from 14 variable positions in 302 individuals (Table 1). The four of most frequent haplotypes (nos. 3, 5, 2 and 17 were present in respectively 34%, 15%, 13% and 8% of all individuals investigated. The overall nucleotide diversity (π), i.e. the percentage of the average number of nucleotide differences per site between two sequences (Nei 1978) was π=0.49%, while the overall haplotype diversity (h) was 0.83. The corresponding values within predefined regions are shown in Table 1.

Habitat quality and distribution

The habitat suitability map (Fig. 1) shows a highly diverse mosaic of habitat qualities at countrywide or regional scales, indicating that the HS model was able to discriminate between suitable and unsuitable habitats and that the set of eco-geographical variables allowed us distinguish specific habitats preferred by great bustards from the overall habitat Lowest-quality habitats were also steppe-like areas, but showing a higher diversity of strata with lower herbaceous cover (greenness) in more abrupt terrain, at higher altitude.

Based on the results of predictive habitat modelling, 83.3% of the entire area studied was defined as good-quality habitat for great bustards, with HSI values ranging from 0.62 (Cáceres) to 0.89 (Castilla), 7.5% of the area was of intermediate quality between 0.19 in Toledo West and 0.48 in Madrid South-West, and 9.3% was of poor-quality habitats, with HSI values not significantly different from zero (Córdoba, Sevilla and Huelva) (Table 2). According to our census data, 93.0% of the population represented by our genetic sample lives in good-quality habitats, while 4.5% and 2.5% occur respectively in intermediate and poor habitats.

The HSI differences between the geo-referenced haplotypes are illustrated in Fig. 2. Obviously, the haplotypes 16, 17 and 21 show extremely low habitat quality indices. In order to verify the potential spatial isolation of these three haplotype areas, we calculated spatial 95% confidence ellipses for the pooled haplotypes 16, 17, 21 and the pooled set of the remaining haplotypes based on the sampled individual's geographic coordinates (Fig. 3). In advance, the individuals had been weighted according to the population they belong to.

Correlating genetic variation with demographic variables and habitat quality

We conducted simple correlation analyses to investigate the relationships among nucleotide diversity (π), patch size, population size, population density, and habitat quality (i.e. HSI) (Table 3). Here, we included only populations with at least ten sampled individuals. Habitat quality was highly correlated with patch size (ρ=0.829, P=0.002), population size (ρ=0.909, P<0.001) and population density (ρ=0.773, P=0.005) indicating the reliability and predictive power of the habitat suitability model. Genetic diversity was significantly correlated with habitat quality (ρ=0.620, p=0.042), size (ρ=0.624, p=0.040) and density (ρ=0.756, p=0.007) of the population, but not with patch size (ρ=0.436, p= 0.180). Among these variables, habitat quality was the only extrinsic factor potentially influencing the population development including population size and density. Therefore, partial rank correlation coefficients were calculated controlling for HSI. The resulting partial correlation coefficients regarding genetic diversity were not significant, as opposed to their non-partial equivalents (Table 3). These results suggest the HSI to be one of the underlying variables behind the pairwise correlations between genetic diversity and the population characteristics.

Discussion

Demographic and ecological correlates of genetic variability

Consistent with the classic island biogeography theory [START_REF] Frankham | Relationship of genetic variation to population size in wildlife[END_REF], i.e. that population size and density affect genetic variability, we found significant positive correlations between these demographic parameters and genetic variability. We hypothesized that variations in habitat quality would have an influence on the amount and spatial distribution of genetic variation. Our results indeed confirmed that genetic diversity was significantly correlated with habitat quality. Although that correlation does not imply causation, it suggests that habitat quality may be an important factor determining changes in great bustard numbers, and also causing variations in the genetic diversity and structure of the population studied. We do not believe that this result could be influenced by two caveats concerning the partitioning of genetic diversity in our study. First, the estimate of genetic diversity at each study site is only a crude approximation due to the small sample sizes, implying that estimates of π might change with additional sampling. Second, the fact that we intentionally sampled several individuals in relatively small areas separated from each other by considerable distances, may have exaggerated the similarity within sample regions and the differences between regions. Attempts to use measures of habitat quality to predict current genetic variation in wild populations are in their infancy. Thus far, they rarely appear to be able to explain the partitioning of diversity in natural populations [START_REF] Angers | Canonical correspondence analysis for estimating spatial and environmental effects on microsatellite gene diversity in brook charr (Salvelinus fontinalis)[END_REF][START_REF] Costello | The influence of history and contemporary stream hydrology on the evolution of genetic diversity within species: an examination of microsatellite DNA variation in bull trout, Salvelinus confluentus (Pisces: Salmonidae)[END_REF]. This might reflect our incomplete understanding of the complex nature of the relationship between habitat quality and population dynamics. Our findings of a significant positive correlation between HSI values and demographic data may be attributed to the selection of appropriate environmental variables to build the habitat suitability model for great bustards in Spain. Currently, a wide variety of methods for modelling species habitat suitability is available [START_REF] Franklin | Mapping species distributions[END_REF]. Many of these methods require observations of species presence and absence, which should be widely distributed across space and environmental gradients, to characterize complex response functions (Guisan et al. 2002). Others are based on only-presence data [START_REF] Elith | Novel methods improve prediction of species' distributions from occurrence data[END_REF]Pearce and Boyce 2006). In general terms, according to [START_REF] Brotons | Presenceabsence versus presence-only modelling methods for predicting bird habitat suitability[END_REF], when presence-absences are available, a modelling approach based on a binary response variable is preferable. The cause is that they generally give more accurate predictions than models based on only-presence data. It should be noticed that one of the most reported troubles when dealing with these methods is the quality of absences, which use to be actually 'peudoabsences', mainly because of detectability problems. However, when absences are 'safe' (as it is the case of the species showing a well known distribution pattern, e.g. great bustard in Spain), the application of methods based in presences-absences become more powerful. Among 'presence-absences methods', GAM has been recommended as one of the best choices for species distribution modelling and spatial prediction because it tends to have high prediction accuracy, being a powerful tool able to model non-linear species responses curves (Wood 2006;Meynard and Quinn 2007).

HSI models can be valuable tools, but should always be applied with a critical evaluation of the model assumptions and uncertainties in relation to the research or management problem (Van der Lee et al. 2006). Despite their widespread use, HSI models can be justified only if they are tested against population measures such as density, reproductive success or genetic diversity and structure [START_REF] Fulgione | Landscape fragmentation and habitat suitability in endangered Italian hare (Lepus corsicanus) and European hare (Lepus europaeus) populations[END_REF]). In the future, similar comparative HSI-based ecological and genetic studies with other taxa would help assessing whether there is really a relationship between genetic variation and sensitivity to local habitat conditions, as our results suggest.

Implications for conservation

The recorded genetic diversity of π=0.49% in Spain is similar to that of the great bustards in the Ponto-Caspian steppes (Ukraine and the Lower Volga basin) (π=0.68%; Pitra et al. 2007). Spanish great bustards form a metapopulation of several subpopulations with various degrees of interconnection among them (Martín et al. 2002;[START_REF] Alonso | Genetic diversity of the great bustard in Iberia and Morocco: risks from current population fragmentation[END_REF]. The patch sizes occupied by these subpopulations are significantly correlated with the respective bird numbers. Because 92% of the total number of individuals live in large patches (>250 km 2 ) with good habitat quality, the demo- graphic stability of the metapopulation seems to be secured, as long as habitat quality does not change in these patches. However, the finding that a substantial proportion of the total genetic diversity of Iberian great bustards is currently restricted to low-quality areas was of primary importance. We discuss these patterns and their likely implications for conservation in more detail below. The different Spanish regions have been subjected to various degrees of habitat deterioration mainly due to agricultural transformations within the last decades [START_REF] Alonso | Status and recent trends of the great bustard (Otis tarda) population in the Iberian peninsula[END_REF][START_REF] Alonso | Distribution dynamics of a great bustard metapopulation throughout a decade: influence of conspecific attraction and recruitment[END_REF]. The results of our genetic study reveal that the existing genetic inventory of the Spanish population may be threatened if habitat deterioration proceeds. Based on the spatially nonrandom distribution of the present genetic variation, two risk factors for ongoing genetic erosion were identified: small size of some populations and poor quality of some habitats.

For instance, the populations in Navarra and Aragón are very small (30 and 66 individuals) and show extremely low genetic diversity values (Table 1). These small and genetically impoverished populations might suffer from inbreeding depression and increased probabilities of extinction [START_REF] Frankham | Inbreeding and conservation: a threshold effect[END_REF]. The populations in Andalusia (Huelva, Córdoba and Sevilla) show moderate genetic diversity values (0.22-0.43) but very low HSIs (0.003-0.054). These populations were much larger in the past, and the low HSI values reflect a poor habitat quality caused by agriculture intensification during last decades [START_REF] Alonso | Status and recent trends of the great bustard (Otis tarda) population in the Iberian peninsula[END_REF][START_REF] Alonso | The great bustard Otis tarda in Andalucía, southern Spain: status, distribution and trends[END_REF]. Since genetic diversity in neutral markers is lost at greater rates in smaller than larger populations (Montgomery et al. 2000), this could explain why populations in Andalusia would have reacted slowly to habitat changes. Lane and Alonso (2001) and [START_REF] Alonso | The great bustard Otis tarda in Andalucía, southern Spain: status, distribution and trends[END_REF] have reported local extinctions within the last two decades and population viability analyses under the status quo conditions inferred 85-99% extinction probability within 100 years for local breeding groups in Andalusia. Given the precarious status of the species in this region, our finding that three matrilineages (haplotypes 16, 17 and 21) are restricted to poor-quality habitats in Andalusia is alarming. Assuming the worst case scenario, the extinction of the Andalusian population would cause a loss of three (15%) of a total of 20 existing matrilineages from the species' current genetic pool in Spain. Consequently, strict conservation measures should be directed at populations located in low-quality habitats. Moreover, management actions should focus on the protection of local populations to promote conservation of the current genetic diversity and persistence of the metapopulation structure. The distribution of great bustards in Spain is highly fragmented, and vacant habitat patches may occur for a variety of reasons, including the size and quality of habitat areas as well as the species' very strong fidelity to traditional lek sites and conspecific attraction (Lane et al. 2001;[START_REF] Alonso | Status and recent trends of the great bustard (Otis tarda) population in the Iberian peninsula[END_REF][START_REF] Alonso | Distribution dynamics of a great bustard metapopulation throughout a decade: influence of conspecific attraction and recruitment[END_REF]. According to present knowledge, maintaining a dense network of suitable patches and maximizing local habitat conditions in order to improve reproductive output, and thus the number of potential dispersers, is the key to secure the persistence of the great bustard metapopulation in Spain and in other similarly fragmented landscapes.

  of the haplotypes H_18 and H_19 not found in Spain b The sequenced fragments correspond to positions 280-945 in the control region of the complete Ciconia ciconia mitochondrial genome (AB026818). Letters designate base substitutions. Nucleotide position numbers of the great bustard mtDNA correspond to positions in the alignment *P<0.05 **P<0.01 n.d. Not determined because of limited sample size. available in Spain. The highest-quality habitats were wide areas showing gentle terrain undulations covered by unobstructed cereal fields mosaics with a low diversity of strata, but a high herbaceous cover (intermediate values of overall greenness), which may offer adequate concealment for birds.

Fig. 1

 1 Fig.1Sample locations plotted on a map of suitable habitats for great bustard in Spain derived from Suárez-Seoane et al.(2002). The habitat suitable index ranges from 0 (unsuitable) to 1 (highly suitable). Predefined regions (ellipses) are: Albacete (1), Aragón (2), Burgos (3), Cáceres (4), Castilla (5), Córdoba (6), Huelva (7), Madrid N (8),

Fig. 3 Fig. 2

 32 Fig.3The 95% ellipses for the pooled haplotypes 16, 17 and 21 (rhombus-shaped) and the pooled remaining haplotypes (circles). The ellipses are centred on the sample means of the respective longitudes and latitudes. The standard deviations of longitude and latitude determine the major axes, and the covariance between longitude and latitude, the ellipses' orientations

Table 1

 1 Geographic distribution of mtDNA haplotypes, number of individuals (N), haplotype diversity (h), and nucleotide diversity (π) within predefined regions throughout the breeding range of the great bustard in the Iberian peninsula

		Albacete Aragón Burgos Cáceres Castilla Córdoba Huelva Madrid
		Accession	no.
	Variable sites b	58111122222445	70333744569392	356056248358
	Haplotype a

Table 3

 3 Rank correlations of genetic diversity with habitat quality and demographic variables across great bustard study sites

		Genetic diversity	Patch size	Population size	Population density	Habitat quality
	Genetic diversity	-	0.436	0.624*	0.756**	0.620*
	Patch size	-0.211	-	0.961**	0.519	0.829**
	Population size	0.191	0.890**	-	0.682 *	0.909**
	Population density	0.591	-0.342	-0.078	-	0.773**
	Correlation coefficients (Spearman's r) above and partial rank correlation coefficients controlling for HSI below the diagonal. Populations in
	Albacete, Aragón, Burgos, and Toledo West were excluded because of their small sample size (<10 individuals)
	*P<0.05
	** P<0.01

Eur J Wildl Res (2011) 57:411-419

Acknowledgments We thank Javier Alonso, Carlos Palacín, Marina Magaña and Beatriz Martín for collaborating during fieldwork. Additional help was provided by Manuel Morales and Enrique Martín. We also thank Anke Schmidt and Dietmar Lieckfeldt for technical assistance, and P.E. Osborne for his participation in building the habitat suitability model. We also greatly appreciate the valuable comments from the editor and two anonymous reviewers. The study was financed by the Dirección General de investigación (projects PB94-0068, PB97-1252, BOS2002-01543 and CGL2008-02567).


Jozef Noga
email: jozef.noga@savba.sk

Miroslav Urban
email: urban@fns.uniba.sk

Towards the Full Coupled Cluster Model with Triples: A Retrospective after 25 Years

Keywords: coupled cluster theory, triple excitations, non-iterative correction, many body perturbation theory, one-particle relaxation coupled cluster theory, triple excitations, non-iterative correction, many body perturbation theory, one-particle relaxation

The CCSD+T(CCSD) approach was proposed 25 years ago together with a suggestion of a coupled cluster hierarchy including triply excited configurations, known as CCSDT-n. This retrospective recapitulates the essence of those approaches in the contemporary context. Finally, the relation to other developments in this field is briefly outlined with emphasis on work perfomed at the Quantum Theory Project.

More than fourty years ago, using diagrammatic techniques Čížek [1][2][3] worked out the exponential ansatz of Coester and Kümmel [4,5] in a manageable way for many-fermion systems. Thus, he fully opened the door into the world of quantum chemistry for what is nowadays known as the coupled cluster (CC) theory. Restricting ourselves to a single Slater determinant reference |Φ , the exact wave function Ψ is expressed as

|Ψ = exp ( T )|Φ . (1)
Here, in general, T = T1 + T2 + T3 + ... is a global excitation operator. Since, due to the Thouless theorem [6], exp( T1 )|Φ represents again a single Slater determinant, minimization of the energy functional

E( T † 1 , T1 ) = Φ|e T † 1 Ĥe T1 |Φ Φ|e T † 1 e T1 |Φ (3)
is not only equivalent to the Hartree-Fock (HF) method, but, as shown recently, it also provides a very practical and efficient alternative to solve the HF problem [START_REF] Šimunek | AIP Conf. Proc. in press[END_REF][START_REF] Noga | [END_REF]. Consequently, with a HF reference a CC approach restricted to mere singles (CCS) brings zero energy change and the most important contribution to the post HF correction (correlation energy) comes from double (pair) substitutions (i. e. T = T2 ), to which the operator of Eq. ( 2) was restricted in the original work of Čížek. Accordingly, he named the method Coupled Pair Many Electron Theory (CPMET) [1,2].

Undoubtedly, throughout the years the power of the CC approach has been convincingly proved in calculations of highly accurate energies and properties of molecules [9][START_REF] Urban | Methods in Computational Chemistry[END_REF][START_REF] Barlett | [END_REF][12][START_REF] Helgaker | Molecular Electronic Structure Theory[END_REF][START_REF] Crawford | [END_REF][15][START_REF] Barlett | Many-Body Methods In Quantum Chemistry[END_REF][START_REF]Recent Progress in Coupled Cluster Methods Theory and Applications[END_REF]. Although it is not explicitly seen in the original work, in order to obtain the wave function, one is solving a wave equation

H|Φ = E|Φ (4)
with a single reference determinant |Φ and similarity transformed Hamiltonian

H = e -T Ĥe T = Ĥe T C . (5)
The CC equations, either with doubles only [2,[START_REF] Bartlett | [END_REF][19][20] (CCD) or singles and doubles (CCSD) [21] involve steps whose computational complexities formally scale as ∝N 6 where N is the number of molecular orbitals. Including triples in full CCSDT [22] means a jump to ∝N 8 scaling, including full quadruple excitations in CCSDTQ causes another jump to ∝N 10 [23, 24]. More generally, each excitation level of the T operator increases the computational complexity by N 2 , starting with CCS that scales as ∝N 4 , i. e. similarly to the HF method.

A particularly important step was accomplished in QTP by introducing the CCSD with full consideration of single and double excitation operators in the exponential ansatz of Eq. (1) [21], i. e. using T = T1 + T2 . Bearing in mind the Thouless theorem, it is clear that this ansatz introduces the important orbital relaxation effects. Nevertheless, soon after introducing the CC method in "large scale" codes it became clear that in order to obtain higly accurate results inclusion of triple excitations into the wave function ansatz, more precisely to the excitation operator expressed by Eq. ( 2), is inevitable. Contribution due to triple excitations appears already at the second order in the wave function expansion in terms of the many body perturbation theory (MBPT), similarly as for the single excitations when a HF reference is considered. Consequently, the aforementioned excitations contribute to fourth order in the energy expansion. In fact, the first report of the inclusion of triples and also singles into the wave function ansatz derived from Eq. (1) goes back to the work of Paldus, Čížek and Shavitt [25] who introduced an Extended CPMET (E-CPMET) in which in addition to CCD (CPMET) the T1 and T3 containing terms in the expansion of H were restricted to ( Ĥ T1 ) C and ( Ĥ T3 ) C . The performace of this method, including the importance of demanding triples was shown for the BH 3 molecule using minimal basis. The cost and complexity of E-CPMET was essentially the same as if would have the full CCSDT. This fact, combined with fairly limited available computer power prevented the authors from demonstrating the ability of E-CPMET on larger molecular systems. Nevertheless, their work showed that there was a need for development of efficient approaches and algorithms allowing full consideration of triples or eventually even higher excitation operators.

More practical implementations of triples within the framework of MBPT intended for "larger" molecules were reported from a few laboratories at the beginning of eighties [26][27][28][29][30]. At the same time Purvis and Bartlett [21] showed that the fourth order contribution from triples combined with the energy from their newly developed CCSD approach provided a reasonable approximation to full configuration interaction (FCI). Soon after, Lee, Kucharski and Bartlett introduced the simplest iterative extension of CCSD by triples within the framework of the CC theory [31] that is at present denoted as CCSDT-1a. They pointed out that unlike CI, connected triple excitations were the most important new terms to add to CCSD, which already accounts for the disconnected triple and quadruple contributions.

Being involved in similar projects in Bratislava, at this stage the present authors started an active collaboration with Quantum Theory Project (QTP), namely with the group of Rodney Bartlett, where both authors spent some time as research associates in the 1980s of the last century. It was the time when quantum chemists and general users of quantum chemistry methods started to believe that the coupled cluster approach could be useful. Of course, this was also a time when the practical applicability of CC (similarly as other highly correlated ab-initio wave function methods) was restricted to rather small systems due to much more limited computer resources than today.

Independently, in Bratislava we had been developing MBPT codes since the late 1970s [29,32], eventually including an efficient algorithm for triples that implicitly used molecular symmetry [28]. Combining our codes with the program developed at QTP [21] enabled a prompt implementation of the methods described below and helped much in debugging both of the codes.

This retrospective summarizes essential results of the joint effort dating back to 1984 and 1986. There were two directions proposed in papers resulting from this period. First, we suggested a hierarchy of CC methods involving triple excitations known as CCSDT-n methods [33,34] eventually leading to the full treatment within FCCSDT [22]. Secondly, a substantial improvement over the aforementioned MBPT correction for triples to CCSD was offered by a modified a posteriori noniterative correction introduced as CCSD+T(CCSD) [33]. Although some of the proposed approximations are currently not so frequently used in practical applications, some of them remain topical, or served as essential foundations for even better approximations.

The CCSDT-n Hierarchy

The standard way of solving Eq. ( 4) is provided by its subsequent projection onto the reference -giving the energy-and onto excitation subspaces that are created by the action of operator (2) on the reference determinant resulting in a set of non-linear equations that determine the amplitudes of T . These sets of coupled equations are then solved iteratively. The driving force behind the development of the CCSDT-n hierarchy was to keep the slowest step in a single iteration at the computational complexity corresponding to the 4th order MBPT, i. e. ∝N 7 . Little was then known about the higher order corrections and their possible mutual cancellations. Our philosophy was to provide a hierarchy that would be "exactly" formulated with T = T1 + T2 + T3 when the Eq. ( 4) was projected onto single (Φ 1 ) and double (Φ 2 ) substitutions, i. e.

Φ µ | Ĥe T1+ T2+ T3 C |Φ = 0, µ = 1, 2 (6)
while (at least from the formal point of view) gradually improving the exponential expansion when projecting onto triply excited states (Φ 3 ). Hence, in addition to Eq. ( 6) the following equations:

Φ 3 | Ĥ T2 C |Φ = 0, (7)
Φ 3 | Ĥe T2 C |Φ = 0, (8)
Φ 3 | Ĥe T1+ T2 C |Φ = 0, (9)
define CCSDT-1b, CCSDT-2 and CCSDT-3, respectively [33,34]. If we consider the amplitudes in terms of the perturbation theory in a generalized sense [35,36] then Eq. ( 7) involves second order, Eq. ( 8) third order, and Eq. ( 9) fourth order contributions in the wave function expansion. One has to add, however, that third and higher orders are not included in full, i. e. the energy is only correct to fourth order in MBPT. In fact, prior to suggesting this hierarchy, an approach originally called CCSDT-1 had been implemented at QTP by Lee Kucharski and Bartlett [31] which differs from CCSDT-1b in projection onto double excitations, when in CCSDT-1 the disconnected clusters involving T3 are disregarded in Eq. ( 6) and which was also correct to fourth order in the MBPT energy expansion. The latter method was renamed to CCSDT-1a in Ref. [33]. Implementation of the CCSDT-2 and CCSDT-3 has been reported in Ref. [34] together with results for BH, HF, NH 3 , CO, C 2 H 2 at equilibrium geometries, and H 2 O for equilibrium bond distance (R e ) together with symmetrically stretched bonds with R = 1.5R e and R = 2.0R e . Except for CO and C 2 H 2 the FCI benchmarks were available. Somewhat surprisingly, for equilibrium geometries the methods that are "higher" in the CCSDT-n hierarchy provided results which deviated from the FCI ones more than CCSDT-1. Nevertheless, with substantially stretched geometry (2R e ), CCSDT-3 energy was closest to FCI. When we extended the hierarchy with CCSDT-4 -including ( Ĥ T3 ) C into Eq. ( 9)-and to full CCSDT [22] led us to conclude that, generally, linear terms in the exponential expansion lower the energy while nonlinear terms act in the opposite direction. A short note regarding the correction added to our original paper a year later is appropriate at this point. Despite careful checking, it happened that a seventh order contribution of ( Ĥ T 3 1 ) C was accounted twice in CCSDT. This was identified after a discussion and numerical comparisons with Gustavo Scuseria whom one of us (JN) met at the Sanibel Symposium in 1988 and who was just finishing implementation of CCSDT at the Center for Computational Quantum Chemistry, University of Georgia [37]. After a quarter of a century, it is now clear that the main impact of Ref. [33] to the field of practical quantum chemistry was not provided by the suggested CCSDT-n hierarchy, but instead by a rather simple approximation to the energy contribution stemming from triply excited configurations that was simply added 'a posteriori' to the CCSD energy. This approximation was originally named "CCSD+T(CCSD)", later a shorthand notation "CCSD[T]" was suggested [38].

What was the essence of T(CCSD)? Even though in CCSDT-1 -CCSDT-3 the computational scaling of the most demanding step is the same as in MBPT(4), i. e. ∝N 7 , such a step has to be repeated in the iterative procedure of solving the nonlinear Eqs. ( 6)-( 9) twice in each iteration. This is, of course, not very favorable, and any plausible approximation that would improve upon such unfavorable scaling was more than welcome. Hence, shortly after introducing the CCSD approach it was suggested to "complete" the energy in order to be correct to the fourth order of MBPT simply by adding the triple-excitation contribution from MBPT(4) to the CCSD energy [30]. The aforementioned contribution can be expressed as

E [4] T = -Φ|( T † 3 F T3 ) C |Φ , ( 10
)
where F is the Fock operator. In the spinorbital basis, F can be always made diagonal . Hence, having T3 , this particular contribution involves ∝N T3 algebraic operations, where N T3 is the number of independent triple-excitation amplitudes. Having

Ĥ = F + Ŵ (11)
T3 is obtained from

-Φ 3 |( F T3 ) C |Φ = Φ 3 |( Ŵ T (1) 2 ) C |Φ . (12)
On the r.h.s the T

2 amplitudes come from the first order wave function within MBPT. Let us remind that if F is a diagonal operator, T3 is obtained in a single step simply by dividing the r.h.s. matrix elements by the pertinent energy denominator.

The idea behind CCSD[T] consists in substituting the first order double excitation amplitudes in Eq. ( 12) by converged double excitation amplitudes from CCSD. Formally, and computationally Eq. ( 12) remains the same. However, through the T2 amplitudes arising from CCSD, the inherently hidden higher order information is also transfered to triples. Indeed, this idea had been partially prompted by an incidental observation. During the development of CCSDT-1b, we used a checking output number from Eq. (10) using T2 of the current iteration in Eq. ( 12). In our iteration scheme the first iteration provided exactly E [4] T when we started from the first order MBPT wave function. Surprisingly, when we started from the T2 -amplitudes after converged CCSD equations, we noticed that if we added the initial "check" number to the CCSD energy, we approximately obtained the CCSDT-1 energy (in this context, "1a" or "1b" is immaterial). Hence, the idea of CCSD[T] was born. In the reasoning supporting this approach we have argued in terms of a symmetrical energy functional [2] related to the CC ansatz of Eq. ( 1)

E( T † , T ) = Φ|e T † Ĥe T |Φ Φ|e T † e T |Φ = Φ| e T † Ĥe T C |Φ . (13)
Čížek has shown that for exact T Φ| Ĥe

T C |Φ = Φ| H|Φ = Φ| e T † Ĥe T C |Φ . (14)
With T = T1 + T2 + T3 and after expanding the exponential in a Taylor series Eq. ( 13) can be rewritten as:

E( T † , T ) = Φ|(1 + T † 1 + T † 2 )( Ĥe T1+ T2 ) C |Φ (15)
+ Φ| T † 3 ( F T3 + Ŵ T2 ) C |Φ (16)
+ Φ|( T † 2 Ŵ T3 ) C |Φ (17)
+ Φ|( T † 2 F T3 ) C |Φ (18)
+ Φ|( T † 1 Ŵ T3 ) C |Φ (19)
+ ...other fifth and higher order terms

The first term (15) corresponds to the asymmetric functional of Arponen's type [39] from which the CCSD equations can be derived by making it stationary with respect to T † 1 and T † 2 . Hence, if we start from CCSD amplitudes this term is exactly equal to the CCSD energy. In T(CCSD) we replace in Eq. ( 12) the first order amplitudes by those from the converged CCSD. Due to the validity of Eq. ( 12) the contribution ( 16) is zero and the term (17) becomes equivalent with Eq. [START_REF] Urban | Methods in Computational Chemistry[END_REF]. What remains (for closed shell systems or the UHF reference) are fifth and higher order terms.

Independently, and in parallel to our work on T(CCSD), Raghavachari [40] used a non-iterative correction due to single and triple excitations added to the CCD energy. The method was named CCD+ST(CCD) and the a posteriori energy corrections were calculated exactly as in our case, using fourth order formulae with converged T2 amplitudes, which has been often an improvement over the CCD+ST(MP4) used earlier at QTP [9,21]. The better performance of T(CCSD) originates mainly from the much better recovering of the orbital relaxation (vide infra) that results from the interplay between the T1 and T2 amplitudes.

Towards a more accurate approximation

The CCSD[T] approximation was the forerunner of the currently widely used CCSD(T) approximation [41] -often denoted as a "golden standard" of ab initio quantum chemistry. In the latter the contribution of ( 19) [E (5) ST ] is taken into account in addition to the T(CCSD) correction with the closed shell and the UHF references. One of the arguments for using this particular T1 -containing fifth order contribution was the fact that this term participates in the final CCSDT-n energies twice, i. e. as in 5th order MBPT, Retrospective: Towards full CCSDT when T † 1 and T † 3 are implicitly restricted to second order in the MBPT wave function expansion. In CCSD[T] only half of the total contribution is inherently included via CCSD. Thus, the CCSD(T) energy includes all the fifth order terms that are present in CCSDT-1. In fact, there are many fifth order terms [42] that are comparable and often larger [43,44]. In our opinion, the most plausible reason for including of E (5) ST together with T(CCSD) was given by Stanton somewhat later [45]. Instead of focusing on the missing contributions in terms of MBPT he concentrated on a perturbation approach based on the CCSD reference. This automatically accounts for the asymmetric nature of the (approximated) coupled cluster theory and the aforementioned contributions represent the leading-order corrections.

Φ|T † 1 Ĥ T3 |Φ + Φ|T † 3 Ĥ T1 |Φ , (20)
The original T(CCSD) approximation was based on arguments using a symmetrical energy functional, while the amplitudes are obtained using a non-Hermitian similarity transformed Hamiltonian Eq. ( 5), i. e. the symmetry in T(CCSD) is forced. For critical cases, when this approximation starts to be inappropriate, significantly improved results are obtained using the left vector amplitudes, instead of assuming that the amplitudes of T † are equivalent to amplitudes of T [46][47][48].

As follows from the Thouless theorem mentioned in Sec. 1, the action of exp( T1 ) on |Φ in the coupled cluster treatment takes care of the orbital relaxation, i.e. transforms the the original reference to a different Slater determinant that is closer to an optimal reference within the given CC ansatz, while the "closeness" depends on the specific evaluation of T1 . Hence, as mentioned already in the context of CCSD [21] the exp( T1 ) could be included in the reference. Then, however, the |Φ has to be determined self consistently. Consequently, integral transformation is needed in each iteration step. This has been done in the so-called Brueckner doubles theory [49]. Since the T1 =0, with this reference the non-iterative correction for triples [BD(T)] is calculated exactly as in ST(CCD) or T(CCSD), though the actual 5th order content (as compared to the HF reference) differs somewhat from CCSD(T). Comparison of results using CC hierarchy up to triples using both HF and Brueckner orbitals has been given somewhat later as a QTP contribution [50].

Starting from open shell ROHF reference

In CC calculations, the use of an ROHF reference for open-shell molecules was first reported in 1988 by Rittby and Bartlett [51]. For the latter reference or more general reference wave functions created from non-canonical orbitals the CCSD[T] (and accordingly also CCSD(T)) approximation has to be modified. First, if the Brillouin theorem is not obeyed, the occupied-virtual block of F is not zero and hence the contribution of ( 18) is non-zero, as well. If the pertinent block of the Fock matrix is treated as a perturbation then the latter term contributes in the 4th order of the MBPT energy expansion and it has to be included in order to stay within arguments behind T(CCSD). This term is, however, easy to calculate and, computationally, it only requires a small fraction of the cost needed for the evaluation of T3 via Eq. ( 12). More problems are caused by the off-diagonal elements of the occupiedoccupied ( Foo ) and virtual-virtual ( Fvv ) blocks of F . With these terms the accurate solution of Eq. ( 12) requires an iterative procedure involving ∝N 7 steps in each iteration, which is undesired. Alternatively, one can abandon the frame of ROHF orbitals and separately diagonalize the Foo and Fvv blocks. Afterwards, one works with semicanonical orbitals leading to different orbitals for different spins. Careful analysis of the problem of triples within iterative and non-iterative CC procedures was presented by Watts, Gauss and Bartlett [52] who prefer the semicanonical approach. If one decides to stay within the spin-restricted scheme, the T3 in Eq. (12) has to be obtained iteratively in order to arrive at results that are invariant to rotations within the blocks of occupied and/or virtual orbitals [52]. In this case the whole Foo and Fvv blocks are treated as the unperturbed Hamiltionian. For the efficient treatment of triples which takes the advantage of having the same spatial orbitals for different spins in the restricted part of the ROHF reference one can eventually give up the orbital invariance and restrict the procedure to its first iteration fully neglecting the off-diagonal part of F with two options [53]: i) using the orbital energies instead of the full diagonal elements in contributions calculated according to Eqs. ( 10) and (12), which would correspond to "traditional" perturbational viewpoint or ii) using the complete diagonal part of F that is in open shell ROHF different from orbital energies. Although such approaches are not invariant with respect to orbital rotations, it has been numerically demonstrated that the differences from invariant values were small, in particular for CCSD(T) d [53].

One has to realize that even when using the ROHF reference which is a proper wave function of the spin operators, the resulting CCSD wave function and the triples contribution in CCSD(T) is not. Approximate spin adaptation of the double excitation amplitudes used for the triples was suggested by Neogrady and one of us (MU) [53]. For a rigorous treatment of the open-shell molecules and the spin adaptation see e.g. Refs. [54,55].

Prompting alternative ansätze

As mentioned, the T(CCSD) approximation was based on arguments related to the expansion of the expectation value Eq. ( 13) in terms of MBPT. Accordingly, we had also investigated the one-particle properties as expectation values of the pertinent operators [56], showing that the error was relatively small, as soon as singles were included. The paper was in press as I (JN) revisited QTP in 1988 for a short time. One night, we started to discuss the details with Rod Bartlett spending about two hours at the blackboard drawing tens of diagrams, while our discussion went far beyond the initial issue. Not far from the early morning we ended understanding how and to which order the terms mutually cancel. To my surprise, the next day Rod came in with the handwritten first draft of the manuscript on the expectation value CC [XCC(n)] hierarchy [35]. Soon after, this was followed by a fruitful collaboration with Staszek Kucharski resulting in a hierarchy based on the Unitary CC method [57] and methods involving 5th order energy terms, including the fifth order MBPT [43], the first implementation of quadruples within the CC methods, namely within the XCC(5) ansatz [58], or non-iterative corrections [44] based on the generalized PT scheme [36]. Methods including quadruple and higher excitations emerged from QTP subsequently [24,[59][60][61].

Conclusions

Presently, CC methods have developed as a most powerfull method for accurate calculations of molecular properties. The hierarchical structure of the CC theory, the possibility of controlling the applicability of a single determinant reference, e. g. through checking the magnitudes of excitation amplitudes, enables obtaining "good results for a good reason". We believe that the work done 25 years ago at QTP and reminded in this retrospective has significantly contributed to the former statement. Besides, the suggested T(CCSD) approach that had been a forerunner of the pragmatic and widely used CCSD(T), we tried to suggest a systematic hierarchy that eventually led to full inclusion of triples within the coupled cluster treatment. In the mean time, there are certainly new hierarchies suggested. Among them the CCn hierarchy [62,63] is based on the argument the orbital relaxation and hence particularly suited for molecular properties. Another interesting concept is the hierarchy based on the method of moments in coupled cluster theory [64][65][66] that builds a common roof above the variety of CC approximations. Along with the progress related to molecules well represented by a single determinant reference we are witnessing a progression of CC methods applicable for more general molecules, particularly those which show some extent of quasidegeneracy. Either in the treatment via higher excitations [43,46,[67][68][69] or in the treatment via multireference approaches [12,15,[70][71][72][73]. In several of these approaches the ideas put forward in our original paper on approximations to iterative as well to noniterative triples were further extended.

Improved theoretical formulations transformed into sophisticated computer codes with high level of parallelism [74][75][76][77], employing localized orbitals [78][79][80][81], or the possibility of reducing the virtual space by methods like OVOS (also initiated at QTP) [82,[START_REF] Neogrády | Recent Progress in Coupled Cluster Methods Theory and Applications[END_REF] contribute to making gradually larger molecules tractable by the CC methods.

3 .

 3 Non-iterative correction for triples.

Acknowledgements

Since the first scientific stay of one of us (MU) at QTP in 1984 we visited QTP repeatedly in years 1984 -1993 for a longer or shorter periods of time. We appreciate very much the unforgettable creative and friendly social and scientific atmosphere at QTP. We particularly wish to thank Rod Bartlett for his wonderful hospitality and for our fruitful collaboration throughout many years. Currently, the work of the authors in this field has been supported by the Grant Agency of the Ministry of Education of the Slovak Republic and Slovak Academy of Sciences (VEGA projects No. 2/0079/09 and 1/0520/10) and also benefits from the Centers of Excellence program of the Slovak Academy of Sciences (COMCHEM, Contract no. II/1/2007).