Elisa Rebolini
email: erebolini@kjemi.uio.no

Julien Toulouse
email: julien.toulouse@upmc.fr

Andrew M Teale
email: andrew.teale@nottingham.ac.uk

Trygve Helgaker
email: t.u.helgaker@kjemi.uio.no

Andreas Savin
email: savin@lct.jussieu.fr

Calculating excitation energies by extrapolation along adiabatic connections

In this paper, an alternative method to range-separated linear-response time-dependent densityfunctional theory and perturbation theory is proposed to improve the estimation of the energies of a physical system from the energies of a partially interacting system. Starting from the analysis of the Taylor expansion of the energies of the partially interacting system around the physical system, we use an extrapolation scheme to improve the estimation of the energies of the physical system at an intermediate point of the range-separated or linear adiabatic connection where either the electron-electron interaction is scaled or only the long-range part of the Coulomb interaction is included. The extrapolation scheme is first applied to the range-separated energies of the helium and beryllium atoms and of the hydrogen molecule at its equilibrium and stretched geometries. It improves significantly the convergence rate of the energies toward their exact limit with respect to the range-separation parameter. The range-separated extrapolation scheme is compared with a similar approach for the linear adiabatic connection, highlighting the relative strengths and weaknesses of each approach.

I. INTRODUCTION

The calculation of excitation energies in densityfunctional theory (DFT) is nowadays mostly done by means of linear response theory in the timedependent framework.

Linear-response timedependent DFT (TDDFT) [START_REF] Casida | Recent Adv. Density Funct. Methods[END_REF] exhibits an excellent price-performance ratio and is, within the usual adiabatic semi-local approximations, very successful at describing excitations to low-lying valence states. However, these approximations introduce several limitations, especially for the treatment of static correlation [START_REF] Gritsenko | [END_REF], Rydberg and charge-transfer excitations [3,4], and double or multiple excitations [5].

Time-dependent theory is, however, not mandatory for calculating excitation energies, as stated by the Hohenberg-Kohn theorem [6]. Indeed, several time-independent DFT approaches for calculating excitation energies exist and are currently being developed. These include ensemble DFT [7][8][9][10][11][12], ∆SCF [13][14][15][16][17] and related methods [18][19][20][21], or perturbation theory [22][23][24][25] along the standard adiabatic connection using the noninteracting Kohn-Sham (KS) Hamiltonian as the zero-order Hamiltonian.

Range-separated DFT constitutes an alternative to standard KS DFT [6,26] where the physical electronic Hamiltonian is replaced not by an effective non-interacting Hamiltonian but by a partially interacting Hamiltonian that incorporates the long-range part only of the electronelectron interaction [START_REF] Savin | Recent Dev. Appl. Mod. Density Funct. Theory[END_REF][START_REF] Yang | [END_REF][29][30][31]. This partially interacting Hamiltonian corresponds to an intermediate point along a range-separated adiabatic connection, which links the KS Hamiltonian to the physical Hamiltonian by progressively switching on the long-range part of the two-electron interaction, whilst simultaneously modifying the oneelectron potential so as to maintain a constant ground-state density.

In range-separated time-dependent DFT, the excitation energies of the long-range interacting Hamiltonian act as starting approximations to the excitation energies of the physical system and are corrected using a short-range density-functional kernel, in the same manner as the KS excitation energies act as starting approximations in linearresponse TDDFT. Several such range-separated linear-response schemes have been developed, in which the short-range part is described by an approximate adiabatic semi-local density-functional kernel and the long-range linear-response part is treated at the Hartree-Fock [32][33][34][35], multiconfiguration self-consistent field (MCSCF) [34,35], second-order polarization-propagator approximation (SOPPA) [35], or density-matrix functional theory (DMFT) [36] level.

Within the time-independent framework, a standard method for improving upon the excitation energies of the partially interacting Hamiltonian would be to use perturbation theory. However, given that perturbation theory in its standard Rayleigh-Schrödinger based formulation does not keep the ground-state density constant at each order in the perturbation, it has not led to a systematic improvement [37].

In this work, we propose a time-independent alternative method for correcting the excitation energies of the partially-interacting system, based on extrapolation along the range-separated adiabatic connection. Given that the long-range part of the interaction is included in the partially interacting system, its excitation energies constitute better approximations to the energies of the physical system than do the excitation energies of the KS system. The analysis of the Taylor expansion of the energies in the range-separation parameter µ about the physical system (µ → +∞) presented in Ref. 38 shows that the energies of the partially interacting system converge towards their physical limits as µ -2 . Using this information, it is possible to develop a scheme for extrapolating the energies of the physical system from the energies of the partially interacting system by following the ideas of Refs. 39, 40. This extrapolation scheme involves low-order derivatives of the energies with respect to µ and constitutes an alternative to perturbation theory and to range-separated TDDFT [32,41,42].

The extrapolation scheme is also applied to the linear adiabatic connection, where the interaction is scaled by a parameter λ going from 0 to 1, and where the analysis of the excitation energies around λ = 1 provides the required information to improve the estimation of the energies of the physical system from an intermediate point of the connection.

The expression for the energies of the partially interacting system and for their extrapolations are given in Section II for the range-separated and linear adiabatic connections. The extrapolation is subsequently applied to the range-separated energies of the helium and beryllium atoms and of the hydrogen molecule at its equilibrium and stretched geometries; for helium, we also use the linear adiabatic connection. The computational details are given in Section III and the results are discussed in Section IV.

II. ENERGY EXTRAPOLATION ALONG THE RANGE-SEPARATED AND LINEAR ADIABATIC CONNECTIONS

A. Range-separated adiabatic connection

Range-separated DFT uses a partially interacting system, where the long-range part of the Coulomb interaction is included instead of the more traditional non-interacting KS system-see, for example, Ref. 31. In terms of the long-range (lr) electron-electron interaction operator Ŵ lr,µ ee =

1 2 w lr,µ ee (r 12 )n 2 (r 1 , r 2 )dr 1 dr 2 , (1)
where n2 (r 1 , r 2 ) is the pair-density operator and w lr,µ ee (r 12 ) is the error-function interaction

w lr,µ ee (r 12 ) = erf(µr 12 ) r 12 , (2)
the Hamiltonian of the partially interacting system is given by Ĥlr,µ = T + Vne + Ŵ lr,µ ee + V sr,µ Hxc .

The parameter µ controls the range of the separation, with 1/µ acting as a smooth cut-off radius. This Hamiltonian also contains the shortrange Hartree-exchange-correlation potential operator V sr,µ Hxc , whose role is to ensure that the ground-state density of the partially-interacting system

n 0 (r) = Ψ µ 0 |n(r)|Ψ µ 0 , (4)
is equal to the ground-state density of the physical system for all µ. Here Ψ µ 0 is the ground-state wave function of the partially interacting Hamiltonian and n(r) is the density operator. The remaining terms in the Hamiltonian of Eq. ( 3) are the usual kinetic-energy operator T and nuclear-electron interaction operator Vne = v ne (r)n(r)dr.

The eigenvectors and eigenvalues of Ĥlr,µ are the ground-and excited-state wave functions |Ψ µ k and energies E µ k of the partially interacting system

Ĥlr,µ |Ψ µ k = E µ k |Ψ µ k . (5)
These excited-state wave functions and energies provide natural first approximations to the excited-state wave functions and energies of the physical system. For µ = 0, they reduce to the single-determinant eigenstates and associated energies of the non-interacting KS Hamiltonian,

ĤKS |Φ KS k = E KS k |Φ KS k , (6)
while, for µ → ∞, they reduce to the excitedstate wave functions and energies of the physical Hamiltonian

Ĥ|Ψ k = E k |Ψ k . (7)
In Ref. 38, it was shown that the asymptotic expansion of the total energy of state k around the physical system is

E µ k = E k + 1 µ 2 E (-2) k + 1 µ 3 E (-3) k + O 1 µ 4 , (8)
where

E (-2) k and E (-3) k
are the corrections entering at the second and third powers of 1/µ, respectively. Following the scheme proposed in Refs. 39, 40, it is possible to estimate the energy of the physical system E k from the energy of the partially interacting system E k and its first-and second-order derivatives with respect to µ.

¿From the Taylor expansion of the energies when µ → ∞, the first-order derivatives of the energies with respect to µ behave as

∂E µ k ∂µ = - 2 µ 3 E (-2) k + O 1 µ 4 , (9)
around the real system. Inserting this into Eq. ( 8), the exact energies E k can be written as a function of the energies along the adiabatic connection and of their first-order derivative as

E k = E µ k + µ 2 ∂E µ k ∂µ + O 1 µ 3 . ( 10
)
This scheme gives extrapolated energies

E EE,µ k = E µ k + µ 2 ∂E µ k ∂µ , (11)
that are correct up to and including the second power of 1/µ relative to the energies of the physical system. The correction given by the extrapolation scheme vanishes at µ = 0 by construction, but should improve the description of the energies as soon as the interaction is switched on. One should note that the absence of a correction at µ = 0 is only due to the choice of 1/µ k as the basis for the expansion. Other basis functions such as µ 2 /(a + µ 5 ) would lead to a correction at µ = 0 but are not considered in this work.

A more elaborate scheme can be developed by using also the correction E (-3) k and the secondorder derivative. In this case, the first-and secondorder derivatives are given by

∂E µ k ∂µ = - 2 µ 3 E (-2) k - 3 µ 4 E (-3) k + O 1 µ 5 , (12) ∂ 2 E µ k ∂µ 2 = 6 µ 4 E (-2) k + 12 µ 5 E (-3) k + O 1 µ 6 , (13)
and, after eliminating E

(-2) k and E (-3) k
, the extrapolated energies become

E EE2,µ k = E µ k + µ ∂E µ k ∂µ + µ 2 6 ∂ 2 E µ k ∂µ 2 . (14)
Higher-order derivatives should further reduce errors. Additionally, several points along the adiabatic connection could be used to perform the extrapolation to increase the accuracy of the extrapolated energies. However, only first-and secondorder corrections at a single point of the adiabatic connection are considered hereinafter.

B. Linear adiabatic connection

If the linear adiabatic connection is performed, then the partially interacting Hamiltonian is defined as Ĥλ = T + λ Ŵee + V λ where V λ is adjusted to keep the ground-state density constant. This potential can be expressed in terms of the connecting parameter λ as

V λ = Vne + (1 -λ) VHx + Vc -V λ c , (15)
where V λ c enters at second order in λ and is equal to Vc at λ = 1. The energies of the partially interacting system can then be expanded around the physical system as

E λ k = E k + (λ -1)E (1) k + (λ -1) 2 E (2) k + O(λ -1) 3 , (16) where E
(1)

k and E (2)
k are the contributions entering at the first and second power of (λ -1), respectively. As in the range-separated case, by differentiation with respect to λ, it is then possible to extrapolate the energies of the physical system at first order by considering only the correction

E (1) k as E EE,λ k = E λ k + (1 -λ) ∂E λ k ∂λ . (17)
When λ = 0, this extrapolation is equivalent to the first-order correction of Görling-Levy perturbation theory [23,25].

A second-order correction can be obtained by using also the correction

E (2)
k . The first-and second-order derivatives are

∂E λ k ∂λ = E (1) k + 2(λ 1 )E (2) k + O(λ -1) 2 , (18)
∂ 2 E λ k ∂λ 2 = 2E (2) k + O(λ -1), (19)
and the extrapolated energies become

E EE2,λ k = E λ k + (1 -λ) ∂E λ k ∂λ + 1 2 (1 -λ) 2 ∂ 2 E λ k ∂λ 2 . ( 20
)
III. COMPUTATIONAL DETAILS Calculations were performed for the He and Be atoms and the H 2 molecule with a development version of the DALTON program [START_REF] Angeli | DALTON, a molecular electronic structure program[END_REF], see Refs. 44-46. Following the procedure of Ref. 38, a full CI (FCI) calculation was first carried out to get the exact ground-state density within the basis set considered: uncontracted t-aug-cc-pV5Z for He, uncontracted d-aug-cc-pVDZ for Be, and uncontracted d-aug-cc-pVTZ for H 2 . A Lieb optimization of the short-range potential v sr,µ (r) was then performed to reproduce the FCI density with the long-range electron-electron interaction w lr,µ ee (r 12 ). Finally, an FCI calculation was carried out with the partially-interacting Hamiltonian constructed from w lr,µ ee (r 12 ) and v sr,µ (r) to obtain the zerothorder energies and wave functions.

Starting from the analytical form of the fit given in the supplementary material of Ref. 38, it is then straightforward to calculate the analytical derivatives of the energies with respect to µ. In the linear case, a cubic fit of the energies was performed. The extrapolated energies were calculated using Eqs. ( 11), ( 14), (17), and (20).

All the unextrapolated curves shown hereinafter correspond to the curves of Ref. 38.

IV. RESULTS AND DISCUSSION

A. Range-separated adiabatic connection of the helium atom

Ground-state energy

The results of the first-and second-order extrapolation schemes on the ground-state total energy of the helium atom are shown in Figure 1 (top). By construction, the extrapolation has no effect at µ = 0 and the ground-state energy of the KS system is, therefore, unaffected by the extrapolation. However, for µ > 0, the extrapolated energies show a systematic improvement with respect to the unextrapolated ground-state energy, that is, the ground-state energy of the partially interacting Hamiltonian without any correction.

Without extrapolation, a range-separation parameter of about 6 bohr -1 is needed to give an error smaller than 10 mhartree relative to the energy of the physical system. With the first-and second-order corrections added, the same accuracy is achieved with a range-separation parameter of only 2.8 and 1.5 bohr -1 , respectively. µ (bohr -1 ) Excitation energies (hartree) µ (bohr -1 ) Excitation energies (hartree) k = E µ k -E µ 0 calculated without extrapolation (full lines), with first-order extrapolation (dashed) and second-order extrapolation (dot-dashed) as a function of µ. The dotted horizontal lines are the physical energies. The colored regions represent errors of ±10 and ±1 mhartree for the ground-state and excitation energies, respectively.

1 1 S → 2 1 S Extrap. (1st) Extrap. (2nd) 1 1 S → 2
1 1 S → 2 1 P Extrap. (1st) Extrap. (2nd) 1 1 S → 2 3 P Extrap. (1st) Extrap. (2nd)

Rydberg excitation energies

Figure 1 also shows the effects of extrapolation on the lowest Rydberg S (middle) and P (bottom) excitation energies of helium. The convergence of the excitation energies towards their physical limit is overall improved by the first-and secondorder corrections with respect to the unextrapolated curves. In fact, the range-separation parameter required to achieve an accuracy of 1 mhartree is divided by approximately a factor 2 or a factor 3 by the first-and second-order schemes, respectively. For the excitation energies considered here, a range-separation value of 2 and 1 bohr -1 suffices to reduce the error to less than 1 mhartree with the first-and second-order schemes, respectively.

The 1 S and 3 S excitation energies change monotonically with increasing µ. Accordingly, extrapolation provides a systematic improvement, the sign of the derivative pulling the excitation energies towards their physical limits at both first-and second-order levels.

The 3 P excitation energy also changes monotonically with µ and the first-order extrapolation provides therefore a systematic improvement. However, the first-order extrapolated energy does not converge monotonically towards its physical limit (not visible), leading to a slight degradation of the excitation energies around µ = 1.5 bohr -1 at second order.

Finally, the 1 P excitation energy shows a nonmonotonic behavior even before extrapolation, exhibiting a "bump" for small µ that is probably a basis-set effect [38]. In fact, for such small µ, only the very long-range part of the interaction is modified, which is poorly described by Gaussian basis functions. The higher a given state is in energy, the more sensitive it becomes to this basis-set defect.

As a consequence, the 1 P excitation energy approaches its physical limit from above, its firstorder derivative changing sign around 0.7 bohr -1 . In this region, the extrapolated energies become less accurate than the unextrapolated energy. However, this behavior is observed only in a small region. As soon as the excitation energy recovers a monotonic convergence towards its physical limit (for µ larger than 0.7 bohr -1 ), the energy is improved by the extrapolation and converges faster to its physical limit.

B. Range-separated adiabatic connection for the valence excitation of the beryllium atom

The ground-state energy of the beryllium atom is shown in Figure 2 (top). Since beryllium has a core orbital, the convergence of its total energies µ (bohr -1 ) Excitation energies (hartree) k = E µ k -E µ 0 (bottom) of beryllium as a function of µ. The unextrapolated energies are shown as full lines, the first-order extrapolated energies are plotted in dashed lines and the secondorder ones in dot-dashed lines. The energies of the physical system are given as horizontal dotted lines. An error of ±50 mhartree is colored around the physical ground-state energy and an error of ±2 mhartree is colored around the physical excitation energies.

1 1 S → 1 1 P Extrap. (1st) Extrap. (2nd) 1 1 S → 1 3 P Extrap. (1st) Extrap. (2nd)
is slower than for helium as the density is more contracted and a larger range-separation parameter is needed to describe correctly the core region. However, this affects all valence states in a similar fashion (not shown here).

Extrapolation systematically improves the convergence of the ground-state energy along the adiabatic connection. First-order extrapolation reduces the error to less than 50 mhartree with µ ≈ 5 bohr -1 , an order of magnitude smaller than the error without the extrapolation correction but still large. Second-order extrapolation gives the same error reduction already with µ ≈ 3 bohr -1 .

The effect of the extrapolation on the valence excitation energies of beryllium is shown in Figure 2 (bottom). As the errors associated with the core largely cancel in the excitation energies, the unextrapolated excitation energies already converge faster than do the total energies. With first-order extrapolation, an error smaller than 2 mhartree is reached with µ ≈ 0.5 bohr -1 , to be compared with a much larger error of 4 hartree in the total energies with the same µ value. The second-order extrapolation allows one to reach the same accuracy with a range-separation parameter as small as 0.3 bohr -1 . However, once again, some bumps are observed in the extrapolated energies probably due to the limited size of the basis set. This fast convergence of the excitation energies with respect to the range-separation parameter is due to the fact that in beryllium, static correlation is important and the multi-configurational character of the wave function is quickly established when the interaction is switched on; see Ref. 47.

C. Range-separated adiabatic connection for the hydrogen molecule

Finally, we consider extrapolation of the lowest excitation energies of the hydrogen molecule along the range-separated adiabatic connection, at the equilibrium geometry and at a stretched geometry. The results of the first-and second-order extrapolations on the singlet and triplet Σ + g → Σ + u excitation energies at the equilibrium geometry are shown in Figure 3 (top). First-and second-order extrapolations provide a systematic improvement in the excitation energies, µ ≈ 2 bohr -1 for first order and µ ≈ 1 bohr -1 for second order being sufficient to reproduce the physical energies to within 1 mhartree.

Having stretched the hydrogen molecule to three times the equilibrium distance, we apply extrapolation to the singlet and triplet excitations to the 1Σ + u state and to the double excitation to the 2Σ + g state-see the bottom part of Figure 3. Again, the improvement is systematic. The triplet extrapolated energy shows a monotonic behavior with respect to µ, whereas the singlet energy shows a slight bump at 0.8 bohr -1 . However, all extrapolated excitation energies converge faster than their unextrapolated counterparts. Extrapolation works remarkably well, reducing errors to less than 5 mhartree with µ ≈ 0.6 bohr -1 , compared with 2 bohr -1 without extrapolation. In particular, extrapolation allows us to describe double and single excitation energies equally well. In this case, one should note that the second-order scheme does not improve significantly the convergence of the µ (bohr -1 ) Excitation energies (hartree) Req µ (bohr -1 ) Excitation energies (hartree) 3Req

1 1 Σ + g → 1 1 Σ + u Extrap. (1st) Extrap. (2nd) 1 1 Σ + g → 1 3 Σ + u Extrap. ( 1st
1 1 Σ + g → 2 3 Σ + g (σ + u ) 2 Extrap. (1st) Extrap. (2nd) 1 1 Σ + g → 1 1 Σ + u Extrap. (1st) Extrap. (2nd) 1 1 Σ + g → 1 3 Σ + u Extrap. (1st) Extrap. (2nd)
Figure 3: (Color online) Unextrapolated (full lines), first-order extrapolated (dashed lines) and secondorder extrapolated (dot-dashed lines) excitation energies ∆E µ k = E µ k -E µ 0 of the H2 molecule at the equilibrium internuclear distance Req (top) and three times the equilibrium distance (bottom) as a function of µ. The excitation energies of the physical system ∆E k = ∆E µ→∞ k are plotted as horizontal dotted lines. An error of ±1 mhartree is colored around the physical excitation energies at equilibrium and an error of ±5 mhartree at stretched geometry. The total energies of the helium atom along the linear adiabatic connection are plotted in Figure 4 (top). When λ = 0, no interaction is included, so the KS energies are recovered as for µ = 0. When λ = 1, the full interaction is present and the energies of the physical system are recovered, which corresponds to the limit µ → ∞. The two limiting cases are, therefore, identical for the two adiabatic connections but the way they are connected differs.

1 1 Σ + g → 2 3 Σ + g (σ + u ) 2 and 1 1 Σ + g → 1 1 Σ + u excitation
The evolution of the total energies with respect to λ is almost linear. Although this behavior is easier to predict and should provide an efficient framework for extrapolations, the value of λ required to have an error of 10 mhartree is very close to 1; in the range-separation case, an intermediate value of µ is sufficient. In general, we note that, in the linear case, the calculation of the wave function is by and large equally expensive at all points along the connection since the electron-electron cusp must always be described (except at λ = 0); in the range-separated case, by contrast, the calculation becomes less expensive with decreasing µ since the description of the cusp is then avoided.

Excitation energies and comparison with the range-separated case

The lowest (Rydberg) excitation energies of helium along the linear adiabatic connection are also plotted in Figure 4 (middle and bottom). As for the total energies, the end points are the same as in the range-separated case but the behavior of the energies along the connection is more linear. As in the range-separated case, the 1 P excitation energy does not evolve monotonically with λ, probably because of basis-set limitations.

When the first-and second-order extrapolation corrections are added, a systematic improvement is observed for the triplet excitation energies. With the second-order correction, the physical energies are already reproduced to 1 mhartree at λ = 0. The singlet excitation energies are less affected by the correction but are still overall improved, the amount of interaction required to reproduce the physical limit within an accuracy of 1 mhartree dropping to 50%. Moreover, unlike in the range-separated case, the KS excitation energy also benefits from this correction, which no longer vanishes in this limit. Indeed, at λ = 0, the excitation energy matches the results obtained in Ref. 25 Figure 5: Error in the 3 S excitation energy along the range-separated and linear adiabatic connections for the helium atom as functions of µ and λ. The unextrapolated energies are given in full lines and the extrapolated ones in dashed lines. An error smaller than 1 mhartree around the physical limit is given by the colored region.

perturbation theory.

To compare the range-separated and linear cases, the effects of first-order extrapolation on the lowest excitation energy of helium are shown in Figure 5 in both cases. Without extrapolation, the scaling parameter λ must be greater than 0.95 to reproduce energies within 1 mhartree. By contrast, a small change in µ near zero gives a large change in the energies; thus, a range-separation parameter of 2 bohr -1 is sufficient to ensure the same accuracy. Clearly, the range-separated connection includes the most significant region of the interaction first, whereas the linear connection treats all ranges equally, independently of their importance for the excitation energies.

For the KS system, it is obviously better to use the correction obtained from the linear connection as the corresponding range-separated correction vanishes. For a partially interacting system, the comparison is more difficult.

V. CONCLUSION

In this work, we have exploited the asymptotic behavior of the energies of a partially interacting system along the range-separated adiabatic connection to design an energy correction that allows us to extrapolate to the physical energies of the system from its partially interacting energies. The simplest possible extrapolations were obtained by using either only the first-order derivative of the energies with respect to the range-separation parameter at a given point or by using the first-and second-order derivatives.

This extrapolation scheme was tested at the FCI level of theory on the helium and beryllium atoms and on the hydrogen molecule (at equilibrium and at stretched geometry), where it significantly improves the convergence of energies and excitation energies towards their physical limits. Moreover, the improvements are systematic, except at µ = 0 (where the correction is zero by construction) and in a few cases where the partially interacting energies present a bump for small µ. In all cases, with respect to the unextrapolated case, the extrapolation schemes reduce the smallest value of the range-separation parameter required to reproduce the physical energies of the system with a given accuracy by approximately a factor of 2 with the first-order scheme and by a factor of 3 with the second-order scheme. This is of particular relevance for truncated wave functions, as the smaller the range-separation parameter is, the fewer Slater determinants are needed to describe the wave functions with an equivalent accuracy.

Finally, the extrapolation scheme was applied along the linear adiabatic connection, where it also improves significantly the description of the excitation energies along the connection.

All results discussed here were obtained without the use of approximate functionals. The proposed extrapolation scheme should now be tested in a more pragmatic case, where the potential is not obtained by Lieb optimization but from different approximations such as the (semi)local approximations or more interestingly approximations where the long-range behavior of the potential is correct, such as the optimized-effective potential (OEP) [48,49]. The effects of the inclusion of higher-order derivatives and of multiple points on this extrapolation should also be explored.

Figure 1 :

 1 Figure 1: (Color online) Helium ground-state energy (top) E µ 0 , and first S (middle) and P (bottom) excitation energies ∆E µk = E µ k -E µ 0 calculated without extrapolation (full lines), with first-order extrapolation (dashed) and second-order extrapolation (dot-dashed) as a function of µ. The dotted horizontal lines are the physical energies. The colored regions represent errors of ±10 and ±1 mhartree for the ground-state and excitation energies, respectively.

Figure 2 :

 2 Figure 2: (Color online) Ground-state energy E µ 0 (top) and excitation energies ∆E µk = E µ k -E µ 0 (bottom) of beryllium as a function of µ. The unextrapolated energies are shown as full lines, the first-order extrapolated energies are plotted in dashed lines and the secondorder ones in dot-dashed lines. The energies of the physical system are given as horizontal dotted lines. An error of ±50 mhartree is colored around the physical ground-state energy and an error of ±2 mhartree is colored around the physical excitation energies.

Figure 4 :

 4 Figure 4: (Color online) Helium ground-and excitedstate energy (top), S excitation energies (middle), and P excitation energies (bottom) calculated without extrapolation (full lines), with first-order extrapolation (dashed), and second-order extrapolation (dotdashed) along the linear connection. The dotted horizontal lines are the physical energies. The colored regions represent errors of ±1 mhartree for the excitation energies.

  energies because of their nonmonotonicity probably due to the limited basis set.

		-1				1 1 S
						2 3 S
						2 1 S
	Total energies (hartree)	-2.5 -2 -1.5				1 3 P 1 1 P
		-3
		0	0.2	0.4	0.6	0.8	1
				λ
		0.77
	Excitation energies (hartree)	0.73 0.74 0.75 0.76				1 1 S → 2 1 S Extrap. (1st) Extrap. (2nd) 1 1 S → 2 3 S Extrap. (1st) Extrap. (2nd)
		0.72
		0	0.2	0.4	0.6	0.8	1
				λ
		0.79
						1 1 S → 1 1 P
						Extrap. (1st)
	Excitation energies (hartree)	0.77 0.775 0.78 0.785				Extrap. (2nd) 1 1 S → 1 3 P Extrap. (1st) Extrap. (2nd)
		0.765
		0	0.2	0.4	0.6	0.8	1
				λ

  , using first-order Görling-Levy

				λ
		0	0.2	0.4	0.6	0.8	1
	Error in the excitation energies (hartree)	0 0.005 0.01 0.015 0.02			E µ -E E λ -E µ Extrap. (1st) λ Extrap. (1st)
		0	2	4	6	8	10
				µ (bohr -1 )

VI. ACKNOWLEDGMENTS

This work was supported by the Norwegian Research Council through the CoE Centre for Theoretical and Computational Chemistry (CTCC) Grants No. 179568/V30 and No. 171185/V30 and through the European Research Council under the European Union Seventh Framework Program through the Advanced Grant ABACUS, ERC Grant Agreement No. 267683. A. M. T. is grateful for support from the Royal Society University Research Fellowship scheme.


Jérôme Fort
email: fort.jerome@gmail.com

Yves Cherel

Ann M A Harding

Jorg Welcker

Dariusz Jakubas

Harald Steen

Nina J Karnovsky

David Grémillet

Geographic and seasonal variability in the isotopic niche of little auks Running head: Isotopic niche of little auks

Keywords: Alcid, Annual cycle, Copepod, Diet, North Atlantic, Pelagic ecosystem, Seabird, Stable isotopes

The non-breeding season of seabirds is extremely challenging to study, because it is often spent offshore under harsh environmental conditions. We used stable isotope analysis to investigate little auk (Alle alle) feeding ecology throughout the annual cycle. The geographic distribution of little auks in the Arctic covers a wide range of oceanographic conditions. We sampled birds from five different colonies located in the most important breeding areas (Greenland and Spitsbergen) to examine how individuals breeding in contrasting marine environments differ in their trophic niche throughout the year. We found differences in summer δ 15 N values amongst the colonies, suggesting different target species despite overall low δ 15 N values in blood which indicates a diet primarily composed of copepods. A rise of δ 15 N values between summer and autumn indicated that adults changed their trophic status to feed at a higher trophic level. During autumn, a large overlap in feathers δ 13 C values between colonies suggests a common moulting area off North-East Greenland.

During winter, the isotopic signatures show that the trophic status of Greenland and Spitsbergen-birds differed, with birds from Greenland feeding at low trophic levels, probably mostly on copepods, while birds from Spitsbergen maintained a higher trophic level. These findings highlight contrasted seasonal and regional diet in little auk populations, and reveal possible population overlaps during their autumn moult. We found substantial trophic variability in little auks, which may indicate unsuspected capabilities to adapt to current, drastic environmental change in the North Atlantic.

INTRODUCTION

In polar and temperate regions, the non-breeding season has a profound impact on seabird ecology. During this period, they have to cope with harsh and highly variable environmental conditions, which can directly affect their body condition, energy balance and survival (Grosbois & Thompson 2005, Rolland et al. 2009) and ultimately shape their population dynamics (Barbraud & Weimerskirch 2003, Daunt et al. 2006). However, the underlying mechanisms and the different factors involved often remain a 'black box' and further studies are therefore needed to explore seabird ecology throughout the non-breeding period. Between breeding seasons, seabirds are typically far offshore and are consequently inaccessible for conventional scientific studies (Gaston 2004). In recent years, technical developments have started to alleviate this problem. For example, miniaturized electronic tracking tags allow investigations of winter movements and migration routes (e.g. Bost et al. 2009, Egevang et al. 2010), fatty acids analyses can provide information on specific dietary components (e.g. Williams et al. 2008), and modelling has improved our knowledge of seabird winter energetics (e.g. Fort et al. 2009). These studies are nonetheless often restricted temporally or spatially, thereby ignoring potential environmental variability which may affect strategies across individuals and populations (Grémillet & Boulinier 2009). Stable isotope analysis is a powerful technique in the study of animal feeding ecology under variable conditions (West et al. 2006). This approach indeed permits the identification of animal isotopic niches (i.e. the δspace used by an organism and defined by the two isotopic axes δ 15 N and δ 13 C; Newsome et al. 2007) where the nitrogen isotopic ratio reflects trophic position/diet of the predator, while the carbon ratio reflects its spatial foraging distribution (Kelly 2000, Hobson et al. 1994).

These isotopic axes therefore help defining the ecological niche of the animal (Newsome et al. 2007, Cherel 2008). In seabirds, stable isotope analysis has been used to track both migration and diet (e.g. Hobson 1999, Cherel et al. 2007). Indeed, since different body tissues incorporate the isotopic signatures of resources at different rates, the analysis of different tissues allows tracking changes in feeding ecology over different time periods and at different time scales, including the non-breeding period (Hobson 1993, Cherel et al. 2008).

The little auk (Alle alle) is a small diving seabird that feeds almost exclusively on zooplankton in Arctic and North Atlantic waters, with a breeding distribution covering a wide range of ocean current regimes and contrasted water masses (Gaston & Jones 1998, Stempniewicz 2001). With recent estimations of more than 80 million individuals (Gaston & Jones 1998, Kampp et al. 2000, Isaksen & Gavrilo 2000, Egevang et al. 2003), this species is the most abundant seabird of the North Atlantic and one of the most abundant seabirds in the world (Stempniewicz 2001). It therefore plays a crucial role within Arctic ecosystems, notably in terms of energy and organic matter transfer. For example, breeding little auks foraging in the North Water Polynya off North-West Greenland account for 92-96% of the carbon flux to seabirds, and consume up to 24% of the copepod standing stock in this region (Karnovsky & Hunt 2002). During the non-breeding season, knowledge of their spatial distribution and diet are more elusive. Vertical migration of their main prey (Calanus spp.; Harding et al. 2008) during winter, which tend to descend to depths unreachable to little auks (>100m; Falk-Petersen et al. 2009), suggests a strong seasonal shift in little auk diet. Only one study (Karnovsky et al. 2008) has investigated the non-breeding diet of little auks, and concludes that these birds show dietary shifts during the different seasons, especially during autumn. However, this study was limited to one particular site, the North Water Polynya, where random birds were sampled during different seasons. This procedure does not allow the comparative analysis of little auk diet across the annual cycle for birds of known origin. Such investigations are nonetheless essential to assess dietary flexibility in little auk populations across a large portion of the North Atlantic, and their capacity to respond / adapt to contrasted and fluctuating environments and resource availability.

Therefore, using isotopic signatures of different body tissues from birds breeding at different geographical sites, the present study had two objectives: (1) to assess trophic level variability in little auks throughout the annual cycle to investigate how birds respond to a change of their feeding environment, and (2) to compare the isotopic niches across little auk populations at different spatial scales in order to determine how environmental conditions experienced by each population may affect little auk feeding ecology during the different seasons.

MATERIALS AND METHODS

Sample collection.

This study was conducted during the 2007 breeding season at four colonies on Spitsbergen and one colony at East Greenland. A total of 291 breeding adults and 20 chicks of little auks were sampled, and 28 chick meals collected (Table 1, Fig. 1).

Hereafter, these five colonies are referred to as EG for the East Greenland colony, and as S1 to S4 following a north-south gradient for the Spitsbergen colonies (Fig. 1).

Blood samples (~0.2 ml) were collected from the brachial vein. Seventy percent ethanol was then added to the whole blood, which was kept frozen at -20°C until isotopic analysis. To check if results were not sex-biased, an additional small amount of blood was taken at EG, S2 and S4 for subsequent molecular sexing as detailed in Fridolfsson and Ellegren (1999).

Little auks have two distinct moults per year: one complete moult in autumn (September-October) that involves the replacement of the complete body plumage, and a partial moult in winter (March) when only feathers from the neck and head are replaced (Gaston & Jones 1998, Stempniewicz 2001). Therefore, two batches of cover feathers were plucked on each breeding adult in summer: one from their body (back or belly, reflecting the autumn period), and one from the head (cheek, neck or throat, reflecting the winter period).

These 2 batches are hereafter called 'body feathers' and 'head feathers'. All feathers were randomly plucked to avoid a bias due to an eventual patchy moult pattern. For chicks, only newly-grown cover feathers were plucked during the late chick-rearing period, when down had totally disappeared. All collected feathers were kept at ambient temperature in sealed plastic bags until analysis.

Chick meals were collected from adult birds, which carry food to the nest in a sublingual (gular) pouch (Stempniewicz 2001). Adult birds were caught in the colony using mist nets or noose carpets. Each food load was gently scooped out of the gular pouch, and immediately preserved in 70% ethanol. Caught adults were released unharmed after 5-10 min of handling. In the laboratory, a random sample of each chick meal was taken for dietary analyses. Each diet sample was divided into broad prey classes (copepods, amphipods and others). Prey items were then numbered and identified after Keast & Lawrence (1990), Kwasniewski et al. (2003), and Walkusz (personal communication). To estimate the composition by mass of the chick meals, we dried items of each taxon and calculated individual dry mass. We then reconstituted the proportion by dry mass of each taxon in the different chick meal samples.

Twenty additional birds legally shot at sea by Greenlandic hunters for commercial sale were bought. Birds were shot off Nuuk (64°10′N 51°45′W; West Greenland, hereafter symbolized as WG) in winter (January 2007) and immediately kept frozen until dissected in the laboratory. Blood samples were collected from the cardiac clot, 'body feathers' randomly plucked from belly, back or throat, and birds were visually sexed from gonads. Blood and body feathers were stored as for breeding birds. No head feathers were sampled from these birds because both body and head feathers in winter plumage originate from the same moult (autumn moult). Stomach contents were also removed, but the advanced digestion stage of the prey precluded identifying them.

Sample preparation and isotopic analysis. Prior to isotopic analysis, blood samples were dried for 48 h at 60°C and homogenized. Feathers were rinsed in a 2:1 chloroform:methanol solution, rinsed 2x in a methanol solution, dried for 48 h at 60°C and homogenized with scissors. Food samples were dried for 48 h at 60°C and ground to powder. Lipids and carbonates were then removed from food samples by rinsing them successively in cyclohexane and 1N HCl. All isotopic analyses were performed by Mylnefield Research Services Stable Isotope Laboratory (SCRI, Dundee, Scotland). Analyses were performed on 1 mg subsamples of dried and homogenized materials loading into tin cups and combusting at 1000°C in an ANCA-GSL elemental analyzer (Sercon, Crewe, UK). Resultant CO 2 and N 2 gases were then analyzed for carbon-13 and nitrogen-15 isotope abundance in continuousflow mode using a SerCon 20:20 isotope ratio mass spectrometer (SerCon, Crewe, UK).

Measured data were scale-calibrated using two international reference materials (IAEA-600 and IAEA-CH6) as well as one in-house standard (leucine). Stable isotope abundances were expressed in δ notation as the deviation from standards in parts per thousand (‰) according to the following equation: δX = [(R sample /R standard )-1] x 1000, where X is 13 C or 15 N and R is the corresponding ratio 13 C/ 12 C or 15 N/ 14 N. The R standard values were based on the Vienna-PeeDee Belemnite (VPDB) for 13 C and atmospheric N 2 (AIR) for 15 N. Data analysis. Using allometric equations between body mass and carbon half-life in avian blood (Carleton & del Rio 2005), stable isotopes half-lives in little auk blood was estimated to 12-15 days. Given this result, and following little auk moulting sequences (see above), blood collected from breeding birds reflects summer isotopic signature (from early to late July, following sampling dates at each colony), body feathers the autumn signature (September-October), and head feathers the winter signature (March) (Bearhop et al. 2002). Blood collected on hunted birds also reflects the winter isotopic niche (early January) of little auks.

In order to compare blood and feather isotopic signatures and, consequently, to investigate the bird isotopic niche at different periods of their annual cycle, we needed to take into account the tissue-dependent metabolic routing and enrichment factors for little auks (Quillfeldt et al. 2008a, Podlesak & McWilliams 2006). We thus compared stable isotope ratios obtained from feathers and blood samples simultaneously collected on chicks, and calculated correction factors to compare δ 13 C and δ 15 N values from these two tissues in adult birds. To avoid effects of growth on δ 15 N blood values (Sears et al. 2009) and therefore to allow comparison of chicks and adults, chicks were sampled during the late chick-rearing period, just before fledging, when growth was reduced (see Harding et al. 2009b). We found differences between chick feathers and blood of -0.25‰ for nitrogen and 0.65‰ for carbon, and used these values as correction factors. Therefore, we corrected adult body and head feathers isotopic data by subtracting the correction factors before comparing them with blood and zooplankton (uncorrected) isotopic results.

Statistics were computed using Statistica 6.0. We used multivariate analysis of variance (MANOVA) with Wilk's lambda statistics to compare simultaneously δ 15 N and δ 13 C values between colonies, as well as analysis of variance (ANOVA) followed by post-hoc range tests (Tukey multiple comparison test for unequal sample sizes) to compare independently δ 15 N and δ 13 C values between sexes, seasons and colonies. Presented values are means ± SD, and statistical significance was assumed at p < 0.05.

RESULTS

For all colonies on which birds were sexed and for all tissues, no difference between δ 15 N and δ 13 C values of adult males and females were observed (ANOVA, δ 15 N: F 5,162 = 1.27, p = 0.28; δ 13 C: F 5,162 = 0.52, p = 0.76). Therefore, results for the two sexes were pooled in the following analyses.

Chick meals and adult δ 15 N values

Chick meals collected at EG and S4 were mostly (63% and 78%) composed of copepods (almost exclusively Calanus hyperboreus and C. glacialis, respectively; Fig. 2).

However, chick meals from EG also included the ice-associated amphipod Apherusa glacialis (32%). The values of δ 15 N measured on chick meals was similar for both colonies (Mann-Whitney test, U = 72.0, p = 0.28, n = 12 and n = 16 for S4 and EG, respectively) (Table 1, Fig. 3). Moreover, chicks and adults from EG had similar δ 15 N values measured in blood samples (t-test, t = 1.56, df = 38, p = 0.13) (Table 1).

Adult breeding birds at EG showed significant variation in their δ 15 N isotopic ratio during the different seasons (ANOVA, F 2,55 = 43.99, p < 0.001). Post hoc Tukey multiple comparison tests indicated an increase of this ratio between summer and autumn (p < 0.001), followed by a decrease between autumn and winter (p < 0.001). δ 15 N isotopic values were similar between summer and winter (p = 0.59) (Fig. 3). Birds breeding in Spitsbergen also present a significant variation in δ 15 N values between seasons (F 2,230 = 541.97, p < 0.001).

Like the Greenland birds, this ratio increased between summer and autumn (p < 0.001).

However, it remained high in winter with a similar mean value to that in autumn (p = 0.100). This winter ratio was significantly higher than in summer (p < 0.001) (Fig. 3).

Seasonal and colony variations

During summer, the different colonies were segregated by their overall isotopic signatures (MANOVA, Wilk's lambda, F 8,174 = 95.50, p < 0.001) and in an univariate analysis, by their δ 13 C (ANOVA , F 4,88 = 282.20, p < 0.001) andδ 15 N values (F 4,88 = 51.15, p < 0.001) measured for blood samples (Fig. 4A). Post hoc Tukey multiple comparison tests indicated that among Spitsbergen birds, the 2 southern colonies (S3 and S4) had higher δ 15 N values than the 2 northern colonies (S1 and S2). Also S3 had higher δ 13 C values than S2 and S4, whereas S4 had lower δ 13 C values than S1. Moreover, δ 13 C and δ 15 N values measured in EG were respectively lower and higher than in all Spitsbergen colonies (all p < 0.001) During autumn, birds from the different colonies were also segregated by their overall isotopic signatures (F 10,224 = 7.31, p < 0.001), as well as by their δ 13 C (F 5,113 = 10.30, p < 0.001) and δ 15 N values (F 5,113 = 4.00, p = 0.002). However, multiple comparison tests indicated that only birds from S3 differed from birds from the other sites, because of higher δ 13 C values (all p < 0.01). Furthermore, only birds from S4 had significantly higher δ 15 N values than those from S3, EG and WG (p = 0.02, p = 0.003 and p = 0.03, respectively) (Fig. 4B).

During winter, birds from all colonies were again segregated by their overall isotopic signatures (F 10,226 = 20.10, p < 0.001) and in an univariate analysis by their δ 13 C (F 5,114 = 12.75, p < 0.001) and their δ 15 N values (F 5,114 = 30.14, p < 0.001). However, multiple comparison tests indicated that δ 13 C values measured on birds from S3 and WG were similar (p = 1.00) but differed from all other colonies (all p < 0.01), whereas δ 15 N values of Greenland birds (EG and WG) were similar (p = 1.00) but differed from all Spitsbergen colonies (all p < 0.001) (Fig. 4C). DISCUSSION Detailed knowledge of seabirds feeding ecology through the non-breeding season is essential to understand the determinants of winter mortality, and how individuals respond to environmental constraints and variability, notably in terms of prey availability. Karnovsky et al. (2008) studied little auk diet and found large seasonal changes throughout the year.

However, their study was restricted to one site (the North Water Polynya off North-West Greenland) and did not take into account intra-individual and spatial variations. Additional investigations were therefore needed to define the trophic niche of this essential seabird component of Arctic marine food webs across the North Atlantic.

The use of stable isotope analysis is particularly relevant in this context, but results and interpretations in terms of distribution and trophic status have to be formulated carefully, notably by integrating information from baseline isoscapes (geographical distribution of isotope values, here δ 15 N and δ 13 C). Indeed, while a difference between measured values might reflect contrasted habitats occupied or prey consumed, it might also result from a variation (in space and/or time) in isotopic values at the base of the food web. A recent overview of isoscapes at large geographical scale (Graham et al. 2010) enables such an integration, and therefore provide guidance for future isotopic research. Little auk distribution is limited to the North Atlantic and the Arctic (above ca. 50ºN; Gaston and Jones 1998), where δ 15 N values are homogeneous (Graham et al. 2010). This homogeneity is confirmed by similar δ 15 N signatures measured on zooplankton (chick meals) from East Greenland and Spitsbergen (Fig. 3), and on adult diet at Hornsund (S4) in a previous year (Harding et al. 2008). This therefore suggests minor spatio-temporal changes in δ 15 N values at the base of the food web within the distribution range of little auks, and allows comparing nitrogen signatures assumed to reflect their relative trophic position on the basis of homogeneous δ 15 N baseline levels. In great contrast, δ 13 C values are spatially structured throughout this region, showing contrasted isotopic habitats. Nevertheless, this geographical gradient is not continuous, and some isotopic habitats represent vast areas where very distant regions (e.g. on both sides of Greenland) share similar baseline values (Graham et al. 2010). Interpretation of δ 13 C signatures due to indicate the movements and distribution of marine predators such as little auks must therefore integrate such essential isoscape information.

Our study shows strong inter-seasonal variation in δ 15 N and δ 13 C isotopic signatures in birds from a given colony, thus confirming the results of Karnovsky and colleagues (2008). Moreover, we observed that little auks from the different study colonies (East Greenland and Spitsbergen) occupy variable isotopic niches throughout the different seasons of their annual cycle, with some important variations between bird populations at different spatial scales.

Interestingly, isotopic niches measured are similar for males and females during all seasons. A previous investigation also found no sexual differences in the isotopic niche of little auks during the summer breeding season (Harding et al. 2008). The present investigation further extends these findings to the autumn and winter. Like guillemots (Uria spp.), little auks switch from bi-parental to uni-paternal care prior to chick fledging (Harding et al. 2004).

Females usually leave the colony earlier, while males accompany chicks during fledging and the first period at sea (Stempniewicz 2001). Given these behavioural differences, we expected autumn moulting grounds and migration routes to be different for the two sexes, therefore leading to contrasted δ 15 N and δ 13 C isotopic values in autumn. Contrary to this expectation, our results show that both sexes have an identical isotopic niche during this season. While this result indicates that males and females are likely to moult in the overall same geographic area, it cannot confirm at smaller spatial scale whether both sexes moult together or within a close area. Further specific studies, using methods complementary to stable isotopes, will have to be performed to investigate sex-specific strategies and potential segregation more accurately.

Spatial variability in the isotopic niche of little auks

The coasts of Spitsbergen and Greenland, on both sides of the Greenland Sea, are characterized by water masses from two highly contrasted sources (Buch 2000; Fig. 1). The West coast of Spitsbergen is dominated by the northward flowing Norwegian Atlantic current which is an extension of the relatively warm North Atlantic current. In contrast, the East Greenland coast is characterized by cooler, less saline water of the East Greenland current flowing southward from the central Arctic basin (Buch 2000). The two distinct summer isotopic habitats (δ 13 C values) of the East Greenland and Spitsbergen birds (Fig. 1) reflect these two different oceanographic environments, and confirm recent isoscape meta-analyses (Graham et al. 2010).

The summer diet of little auks has been described in detail in previous studies (e.g. Karnovsky et al. 2003, Jakubas et al. 2007, Karnovsky et al. 2008, Harding et al. 2008). Here, only chick meals were collected, but the similar δ 15 N isotopic values measured both in chick and adult blood samples strongly suggest that adult and chick little auks feed on similar prey.

Therefore, results obtained in this study (isotopes and prey identification) confirm that all birds from the different colonies and populations mostly feed at low trophic levels and on Calanus copepods in summer. However, several differences are apparent amongst colonies at large and small spatial scales (Fig. 2 and Fig. 4A). For example, prey identification shows that the main copepod species consumed by East Greenland (EG; C. hyperboreus) and Spitsbergen (S4; C. glacialis) birds differ, a result that mirrors prey availability at both sites (Karnovsky et al. 2003, Falk-Petersen et al. 2009) and confirms previous investigations (Harding et al. 2008, 2009a, Kwasniewski et al. in press). Moreover, the presence of sea-ice near the East Greenland colony allowed breeding birds to feed on the ice-associated amphipod Apherusa glacialis (Fig. 2).

In Spitsbergen, the two northern colonies are also segregated from the two southern colonies by their δ 15 N values (Fig. 4A). However, the difference observed is only marginally significant, a fact that might be explained in two different ways. Firstly, and most likely, northern birds may feed at a slightly lower trophic level than southern birds, thereby reflecting different targeted prey. It is established that during summer, adult little auks perform long foraging trips, allowing them to feed in areas over 200 km from their breeding sites (Welcker et al. 2009). Birds from S1 and S2 are therefore able to reach distant feeding grounds in the Arctic Ocean, where multi-year ice prevails as well as associated amphipods Apherusa glacialis (Søreide et al. 2008), which is the prey with the lowest nitrogen isotopic signature consumed by little auks (Harding et al. 2009a, Tamelander et al. 2006). First investigations of little auk feeding grounds, using GPS-tracking, confirmed that some breeding adults from Magdalenefjorden (S1) reach the sea ice edge at 130 km from their colony (Jakubas et al. ). Consequently, it is likely that adults from the two northern colonies feed on a higher proportion of Apherusa glacialis (and therefore at a lower trophic level) than those from the two southern colonies which cannot reach northern, cooler waters. On the other hand, we cannot exclude a slight and regional baseline isotopic variation between North and South Spitsbergen, which may also explain the small difference between both groups, with birds from all Spitsbergen colonies feeding at similar trophic levels. To clarify this point, further, smaller-scale investigations are now required to define δ 15 N baseline values across the marine habitats of little auks from Spitsbergen.

Overall, the foraging distance of breeding little auks is constrained by the need to frequently feed the chick at the colony (Welcker et al. 2009). This constraint drives the behaviour of the birds, with all individuals from a given colony feeding in a similar habitat, and on the same range of prey available in their immediate environment. This phenomenon is reflected by a low variance in δ 15 N and δ 13 C isotopic values during summer (Fig. 4A).

Conversely, during the non-breeding season, birds are free to occupy any area of their winter range, and to feed on any prey item. This is confirmed by our study, which clearly shows that the isotopic niche occupied by little auks from the same colony is wider during the nonbreeding season, as found in other seabird species elsewhere (Cherel et al. 2007(Cherel et al. , 2008)). This strongly suggests that non-breeding little auks disperse across contrasted water masses of the North Atlantic, where they feed on a wider range of prey.

Temporal variability in the isotopic niche of little auks

Between summer and autumn, little auks from all studied colonies showed a drastic change in their trophic status (reflected by δ 15 N isotopic ratio). Indeed, their δ 15 N increased by approximately 3‰ between both seasons. This change of more than one trophic level (Hobson & Welch 1992) strongly suggests a dietary shift between summer and autumn. This change is consistent with that observed in little auks from the North Water Polynya (Karnovsky et al. 2008), and with our current knowledge of copepod ecology (Falk-Petersen et al. 2009). Indeed, soon after summer, copepods (including Calanus spp. species) caught by little auks while breeding (Fig. 2) are known to perform a vertical migration to depths of several hundred meters to undergo diapause, thereby becoming inaccessible to the birds, which can only dive to 40 meters (Harding et al. 2009a). These copepods are thought to migrate back to surface waters a few months later to match the phytoplankton spring bloom (Falk-Petersen et al. 2009, Henson et al. 2009). During autumn, little auks are thus expected to feed on alternate prey. For instance, larger, energy-richer amphipods feeding at higher trophic levels might be consumed by little auks in early autumn (Stempniewicz 2001).

During autumn, birds from all colonies displayed similar δ 13 C values (except birds from S3; Fig. 4B). As mentioned above, and relying on this similarity only, it is difficult to interpret these results in terms of common/different moulting places. However, existing information from long-term ring recoveries (e. The pattern of changes of bird trophic levels between autumn and winter is different for the East Greenland and Spitsbergen colonies. Birds breeding in East Greenland display a second dietary shift in winter since they feed at a lower trophic level than during autumn.

Interestingly, their δ 15 N winter values, similar to those of summer, are in agreement with birds preying upon copepods during this period (Sato et al. 2002, Karnovsky et al. 2008). Little auks from East Greenland mainly winter off Newfoundland (Stempniewicz 2001, Mosbech et al. unpublished data). In this region, phytoplankton spring bloom starts earlier than in the rest of the North Atlantic (i.e. from February;Henson et al. 2009) and it is likely that copepods, which migrate back to surface waters following this bloom, are available to diving little auks during their winter moult (Fort et al. 2010). In contrast, birds from all Spitsbergen colonies adopt a different feeding strategy during winter. Their δ 15 N values show that little auks maintain a high trophic level during the winter that is similar to that in the autumn.

During the winter moult, Greenland and Spitsbergen populations share similar isotopic habitat (reflected by similar δ 13 C isotopic ratios). However, δ 15 N values show that birds feed at different trophic levels (see 'Results'), suggesting contrasted targeted prey, and therefore distinct moulting grounds, where prey availability is different. Interestingly, during all the non-breeding season, birds from S3 clearly adopted a different strategy and moult in different water masses than birds from all other colonies (Fig. 4). Mechanisms underlying this difference remain unknown and merit further investigation, notably by combining geolocation and isotope analysis (e.g. Bost et al. 2009).

The isotopic niche of birds wintering off South-West Greenland

During the autumn, birds collected off South-West Greenland (WG) had identical isotopic niches as those from other colonies (with the exception of S3). This suggests that birds were moulting in a similar isotopic habitat, feeding at the same trophic level and most likely on the same prey. However, during winter, they occupy a different isotopic niche than birds from Spitsbergen and a different habitat (different δ 13 C values) than EG birds, even if feeding at similar trophic levels (similar δ 15 N values). Therefore, it seems that wintering birds caught off Nuuk (WG) represent a different population than those sampled during summer, though our study does not permit to determine their origin. Interestingly, this result also indicates that birds from both East Greenland and Spitsbergen do not winter off southwest Greenland, but rather in different isotopic/feeding habitats, thereby confirming previous investigations based upon ring recoveries and biotelemetry. However, it is important to note that the analysis performed on winter bird blood samples reflects isotopic signatures from January, a period during which birds are still on their wintering grounds. In contrast, analyses performed on breeding bird head feather samples reflect isotopic signatures from the March moult, during which birds already start their spring migration (Stempniewicz 2001).

Consequently, and even if both tissues represent the winter trophic status and isotopic niche, they are separated by a period of about 1-2 months during which birds probably moved, potentially affecting their δ 13 C isotopic ratios. Therefore, these isotopic results are not totally comparable and have to be interpreted carefully, supporting the fact that further studies are necessary to confirm bird location during the different seasons.

Conclusions

Consideration of temporal and spatial variation in seabirds feeding ecology is essential to a better understanding of how they adapt to environmental stochasticity, and of how varying constraints affect their survival. Few studies have investigated seabird diet across seasons using stable isotopes on different tissues (e.g. Cherel et al. 2008, Quillfeldt et al. 2008b), while others have compared different populations during a specific period (e.g. Cherel et al. 2006). However, hardly any study considered both variables simultaneously in a seabird species. Presenting two distinct moults per year, both occurring during a short time window, and with a geographic distribution covering a wide range of current regimes and water characteristics (Stempniewicz 2001), little auks are excellent candidates to study spatiotemporal variation in feeding ecology by way of stable isotopic analyses. Using these particularities, our study is therefore one of the first to track trophic niche variation in a seabird across the annual cycle while comparing strategies adopted by birds from distant colonies (but see Bearhop et al. 2000).

We showed that during their non-breeding season, individual little auks display an important seasonal change in trophic status, and highlighted significant variations between colonies and populations at both small and large spatial scales. Our study also provides one of the first descriptions of the winter feeding ecology of little auks (e.g. Fort et al. 2010) and gives important insights into their moulting and wintering grounds, which need be confirmed with complementary methods such as geolocation (see Bost et al. 2009).

Little auks from the North Atlantic are able to modulate their feeding ecology depending on the temporal availability of their prey and also to exploit a wide range of prey, ranging from calanoid copepods to fish larvae (present study, Karnovsky et al. 2008). Such information is essential to understand how this species responds to a constraining and fluctuating environment. Interestingly, we also hypothesize that such a spatio-temporal dietary flexibility might enable little auks to feed on zooplankton species which currently do not occur in the North Atlantic, but which may extend their northern distribution as a result of global warming (Beaugrand & Reid 2003) and may therefore invade North Atlantic waters and little auk foraging areas. The capacity of little auks to survive upon these new resources may determine the fate of their populations in a warming Arctic.

• AMAP ( 1998

  unpublished data). Moreover, the prey Apherusa glacialis was present in chick meals collected at S1 in2007 and 2008 (Kwasniewski et al. in press) while it has not found at S4 (this study,Kwasniewski et al. in press

  g. Stempniewicz 2001) and recent biotelemetry studies (Mosbech et al. unpublished data) suggested that: (1) Spitsbergen birds move quickly after the breeding season to ice-filled areas of the West Greenland Sea, and start to moult there (Stempniewicz 2001), while (2) East Greenland birds move northerly after summer to moult in the Greenland Sea (Mosbech et al. unpublished data). Combined results therefore indicate that during autumn, birds from East Greenland and Spitsbergen moult in a similar isotopic habitat/water mass, most likely in the Greenland Sea.

Fig. 3 .Fig. 4 .

 34 Fig. 1. Map of the studied little auk breeding colonies which are labelled as follows in the

Acknowledgements. We thank all assistants for their hard work in the field, as well as M. Munck and NANU Travel for their logistical support. We are grateful to C. Egevang for his help in obtaining winter birds, to B. Moe and E. Noreen for genetic sexing, and to W. Walkusz for his advice on zooplankton identification. We also thank D. Ainley, S. Bearhop and S. Votier for their comments on an earlier version of the manuscript. Funding for this study was provided by the French Polar Institute Paul Emile Victor (Grant 388 to D.G.), the National Science Foundation (grant 0612504 to A.H. and N.K. and Grant 0301469 to N.K.) and the Polish Ministry of Higher Education and Science (grant IPY/25/2007 to D.J.). This study was conducted under permits of the Ethics Committee of the French Polar Institute (MP/12/24/05/05) and with the permission of the Greenland Home Rule Government, Ministry of Environment and Nature (Danish Polar Center Scientific Expedition Permit 512-240). J.F. was supported by a grant from the University Louis Pasteur (Strasbourg, France).