text:
# XAFS spectroscopy. I. Extracting the fine structure from the absorption spectra

## I Introduction

X-ray-absorption fine-structure (XAFS), $`\chi `$, is determined by :

$$\chi (E)=[\mu (E)\mu _0(E)]/[\mu _0(E)\mu _b(E)],$$
(1)
where $`\mu `$ is the measured absorption, $`\mu _0`$ is the “atomic” absorption due to electrons of considered atomic level, $`\mu _b`$ is the absorption of other processes. Since the electronic state of an embedded atom is, in general, different from its state in gaseous phase, $`\mu _0`$ is not the same as for isolated atom and cannot be found experimentally. Therefore a demand arises for an artificial construction of $`\mu _0`$.

Usually, $`\mu _b`$ is approximated by a Victoreen polynomial $`P=aE^3+bE^4`$  or by a more general polynomial $`P`$, coefficients of which are found by the least squares method from $`\mu (E)=P(E)`$ at energies lower than the edge.

Further, energy dependence is transformed to the photoelectron wave number dependence: $`k=\sqrt{2m_e(EE_0)}/\mathrm{}`$, where $`E_0`$ is the energy of the corresponding absorption edge. Usually, to the $`E_0`$ the energy at half the step is assigned or the energy of inflection point of $`\mu (E)`$. In most practical works the deviation of $`E_0`$ from true value, $`\mathrm{\Delta }E_0`$, is one of the fitting parameters.

The most difficult procedure in extracting of XAFS from the measured absorption is the construction of $`\mu _0`$ since one cannot definitely distinguish the environmental-born part of absorption from the atomic-like one. All methods for determination of the post-edge background are based on the assumption of its smoothness, and the only criterion for its validity is the absence of low-frequncy structure in $`\chi (k)k^w`$, i. e. the small absolute value of the Fourier transform (FT) $`\rho (r)`$ at low $`r`$. The review of existing post-edge background methods and the propositions of some new is the main purpose of the article.

Special attention must be paid to the estimation of noise and uncertainties in XAFS data. Experimental noise is shown to be essentially smaller than the errors of the background approximation, and it is the latter that determines the variances of structural parameters in subsequent fitting. The corresponding section of the present article is closely related with the next article devoted to the determining the errors of structural parameters .

All described in the article methods for background removal, its error estimations, and XAFS-function corrections are realized in the freeware program viper  which allows one to vary several parameters by hand and watch the results simultaneously.

## II Methods of $`\mu _0`$ construction

### A Smoothing spline

Owing to fast algorithm and easy program realization, the approximation of $`\mu _0`$ by the smoothing spline has become widespread. Let $`N+1`$ experimental values of $`\mu _i`$ are defined on the mesh $`E_i`$. The smoothing spline $`\mu _0`$ minimizes the functional

$$J(\mu _0,\mu )=_{E_{\mathrm{min}}}^{E_{\mathrm{max}}}[\mu _0^{\prime \prime }(E)]^2𝑑E+\frac{1}{\alpha }\underset{i=0}{\overset{N}{}}(\mu _{0i}\mu _i)^2.$$
(2)
The smoothing parameter (or regularizer) $`\alpha `$ is the measure of compromise between smoothness of $`\mu _0`$ and its deviation from $`\mu `$. At $`\alpha =0`$ the smoothing spline exactly coincides with $`\mu `$, at $`\alpha \mathrm{}`$ it degenerates to $`\mu _0=const`$. Optimal regularizer should lead to $`\mu _0`$ containing only low-frequency oscillations and, hence, to $`\chi `$ containing only structural oscillations. The formulation of a new criterion for optimal $`\alpha `$ we shall consider below.

First, we address another problem, the well-known spline instability with respect to the small variations of input parameters: number of nodes, nodal values of the processed function, and limits on integral. In our case the spline is most sensitive to $`E_{\mathrm{min}}`$ due to fast growth of $`\mu `$ in the edge. To raise the stability the method was put forward in viper program which lies in the use of a prior information specifying the shape of $`\mu _0(E)`$ dependence. It is known in advance that the absorption edge without so-called white line constitutes nearly smooth step; the white line, if presents, is added to the step. Denote this prior function as $`p(E)`$. Now we will tend the second derivative of the sought $`\mu _0(E)`$ not to zero (at the specified deviation of $`\mu _0`$ from $`\mu `$) but to the second derivative of $`p(E)`$. The sought $`\mu _0(E)`$ is now minimizes the functional

$$J^{}(\mu _0,\mu )=_{E_{\mathrm{min}}}^{E_{\mathrm{max}}}[\mu _0^{\prime \prime }(E)p^{\prime \prime }(E)]^2𝑑E+\frac{1}{\alpha }\underset{i=0}{\overset{N+1}{}}[\mu _{0i}\mu _i]^2.$$
(3)
As seen, in fact there is no need to know $`p(E)`$ itself, its second derivative is sufficient. The explicit presence of $`p(E)`$ in the following formulas should be taken as a consequence of the technical trick applied: at first $`p(E)`$ is subtracted from the data, then it is added to the found spline.

Represent the second derivatives in finite-difference approximation, introduce $`\stackrel{~}{\mu }_{0i}=\mu _{0i}p_i`$, and denote $`\mathrm{\Delta }_i=E_{i+1}E_i`$:

$$J^{}(\mu _0,\mu )=\underset{i=1}{\overset{N}{}}[\stackrel{~}{\mu }_{0i1}\mathrm{\Delta }_{i1}^1\stackrel{~}{\mu }_{0i}(\mathrm{\Delta }_{i1}^1+\mathrm{\Delta }_i^1)+\stackrel{~}{\mu }_{0i+1}\mathrm{\Delta }_i^1]^2+\frac{1}{\alpha }\underset{i=0}{\overset{N+1}{}}[\stackrel{~}{\mu }_{0i}(\mu _ip_i)]^2=J(\stackrel{~}{\mu }_0,\mu _ip_i).$$
(4)
Thus, the problem is reduced to the preceding one in which instead of initial data $`\mu _i`$ the difference $`\mu _ip_i`$ is appeared. The sought $`\mu _0`$ is found from the smooth $`\stackrel{~}{\mu }_0`$ as $`\mu _{0i}=\stackrel{~}{\mu }_{0i}+p_i`$. In Fig. 1 is shown an example of the atomic-like absorption approximation by the smoothing spline with and without the use of prior function***Here and hereafter for examples is used the spectrum at Bi $`L_3`$ absorption edge in Ba0.6K0.4BiO3 at 50 K recorded in transmission mode at D-21 line (XAS-13) of DCI (LURE,Orsay, France) at positron beam energy 1.85 GeV and the average current $`250`$ mA. Energy step — 1 eV, counting time — 1 s. Energy resolution of the double-crystal Si  monochromator (detuned to reject 50% of the incident signal in order to minimise harmonic contamination) with a 0.4 mm slit was about 2–3 eV at 13 keV.. Energy $`E_0`$ was determined at half the step height. Here, we constructed $`p(E)`$ in the following manner. Found the average value $`\overline{\mu }`$ of $`\mu (E)`$ in region $`20E70`$ eV above the absorptance maximum. Moving from the beginning of spectrum, assign $`p=\mu `$ until $`\mu >\overline{\mu }`$, further $`p=\overline{\mu }`$. Then $`p(E)`$ was smoothed 5 times on 3 points. To perform the Fourier transform, $`\chi (k)k^2`$ was brought into the uniform scale with $`\delta k=0.03`$ Å-1 and multiplied by a Kaiser-Bessel window with parameter $`A=1.5`$. As seen, the use of $`p(E)`$ has led to disappearance of the spurious peak on the absolute value of FT at $`r0.5`$ Å.

So far we have considered the atomic-like absorption $`\mu _0`$ to be a smooth function with no peculiarities. However, in some spectra $`\mu _0`$ itself has a fine structure  originating from resonance scattering within absorbing atom or from multi-electron transitions. If in these cases, based on theoretical calculations, experimental information, or empirical considerations, one can nearly indicate the location of peculiarities, their width and weight relatively to the step height, then one would readily construct the prior function $`p(E)`$ and find the correct $`\mu _0`$. Instead of constant value above absorption edge, the prior function would have corresponding valleys and/or peaks.

Let us now define the criterion for determination of smoothing parameter. An attempt to solve the problem was made in Ref. , where the requirement was proposed: $`H_RH_N0.05H_M`$, where $`H_R`$ is the average value of the weighted Fourier transform magnitude between 0 and 0.25 Å, $`H_M`$ is the maximum value in the transform magnitude between 1 and 5 Å, $`H_N`$ is the average value of the transform magnitude between 9 and 10 Å attributed to the noise. Obviously, that this criterion cannot pretend to the generality since depends on the weighting (op. cit., $`k^3`$) and the relative contribution of noise and the first coordination shell into spectra.

In the program viper we have proposed another approach to the problem based on the consideration of heights of FT peaks as functions of regularizer $`\alpha `$ (see Fig. 2). On increasing $`\alpha `$ from zero, $`\mu _0`$ starts to deviate from the experimental absorption $`\mu `$, $`\rho (r)`$ is growing and then saturates, the peaks at larger $`r`$ being saturated earlier. Clearly, that $`\alpha `$ should be determined by the first peak height since it is the last to saturate. Define the start of saturation on the minimum of second derivative of the first peak squared with respect to $`\mathrm{ln}\alpha `$. Declare the value of $`\alpha `$ in the minimum to be optimal. It is seen that the increase of $`\alpha `$ from the optimal leads to unwanted rapid growth of $`\rho `$ at low $`r`$.

In the example in Fig. 1 the regularizer is optimized following our new criterion.

Unfortunately, the method of smoothing spline does not include any approach to the estimations of uncertainties in the $`\mu _0`$ obtained, in contrast to the following two methods.

### B Interpolation spline drawn through the varied knots

The method was put forward in Ref. . $`N`$ knots are equally spaced in $`k`$ space, through them an interpolation spline is drawn. The ordinates of the knots are varied to minimize $`\rho `$ or $`|\rho \rho _{\mathrm{st}}|`$ in the chosen low-$`r`$ region $`0rr_0`$, where $`\rho _{\mathrm{st}}`$ is the absolute value of the FT of a “standard” $`\chi _{\mathrm{st}}(k)k^w`$, calculated or experimental. The number of knots must not exceed the value $`N_{\mathrm{max}}=2r_0\mathrm{\Delta }k/\pi +1`$,  where $`\mathrm{\Delta }k`$ is the $`k`$ range of useful data. In the Ref.  was asserted that one need to know the “standard” $`\chi _{\mathrm{st}}(k)k^w`$ merely approximately since it used only to get an estimate of the leakage from the first shell to the region minimized. The strange thing is that having omitted the question on the accuracy of found knots (as we show below, rather poor), the authors of the cited work made a fine comparison between several theoretical models for $`\chi (k)`$ calculations.

In Fig. 1 is shown an example of the method application. Ordinates of the 13 knots ($`N_{\mathrm{max}}=13.2`$) were varied to minimize the difference $`\rho \rho _{\mathrm{st}}`$ at $`0r1.05`$ Å. The function $`\chi _{\mathrm{st}}(k)`$ was calculated using feff6 program  (as was pointed above, a crude estimate is sufficient, so details omitted). In the minimized region the $`\rho (r)`$ is somewhat better than that obtained by the previous method. However, at $`k>15`$ Å-1 one can distinguish the obviously wrong behavior of $`\chi (k)k^2`$, and the first peak on $`\rho (r)`$ becomes quite distorted.

Consider now the problem of the accuracy of knot positions $`y_j`$, $`j=1,\mathrm{},N`$ in fitting $`\rho (r)`$ to $`\rho _{\mathrm{st}}(r)`$. As a figure of merit, the $`\chi ^2`$-statistics appears:

$$\chi ^2=\frac{N_{\mathrm{max}}}{M}\underset{m=1}{\overset{M}{}}\frac{[\rho (r_m)\rho _{\mathrm{st}}(r_m)]^2}{\sigma _m^2},$$
(5)
where $`\sigma _m`$ are the errors of $`\rho (r_m)`$. It can be shown (detailed analysis see in the next article ) that under the assumption of uncorrelated knot positions, the mean-square deviation of $`y_j`$ from the obtained through the fit optimal value $`\widehat{y}_j`$ equals $`\delta (y_j)=(\frac{1}{2}^2\chi ^2/y_j^2)^{1/2}`$, where the partial derivatives are calculated in the fitting procedure at the minimum. $`\sigma _m`$ are assumed to be constant and equal to the root-mean-square average of $`\rho (r)`$ between 15 and 25 Å, where solely the noise is present. The errors $`\epsilon _j=\delta (y_j)/[\mu _0(E_j)\mu _b(E_j)]`$ found under such assumptions are shown in Fig. 4 as open circles with the solid line. Notice, that the ussumption that the knot positions are not correlated gives quite optimistic $`\epsilon _j`$. Actually, several first knot positions appear to be highly correlated; the proper taking into account of the correlations (here we do not present these calculations) raises $`\epsilon _j`$ at the least as twice. But even these underestimated $`\epsilon _j`$ are appreciably larger than those given by the following method.

### C Bayesian smooth curve

Ideologically similar to the smoothing spline method is the method of bayesian smoothing (see Appendix on p. Bayesian smoothing and deconvolution) proposed in the program viper. This method also finds the regularized function $`\mu _0`$, the regularizer $`\alpha `$ is the measure of compromise between smoothness of $`\mu _0`$ and its deviation from $`\mu `$. In comparison with smoothing spline method, this method has some advantages. (i) Various prior information on $`\mu _0`$ can be considered. (ii) In this method the posterior distributions of all $`\mu _{0j}`$ are sought for. From those distributions one can find not only average values but also any desirable momenta, which appears to be an additional difficulty for other methods. (iii) In the framework of the method it is possible also to deconvolute $`\mu `$ with the monochromator rocking curve. The weakness of the method is its low speed (comparing with method II A, not with II B!). On a modern PC the curve drawn through $`N500`$ points is smoothed for a few minutes.

In Fig. 5 the bayesian smoothing was done on the mesh of 536 experimental points above $`E_0`$, without and with the prior function (its construction is described in Sec. II A). Besides, in the last case another information was used: the atomic-like absorption must coincide with the total absorption (minus pre-edge background) at energies $`E<E_0`$. Therefore, we demanded from the bayesian curve to pass through a point nearest (at left) to $`E_0`$. The values $`\overline{\mu }_{0j}`$ and $`\delta ^2(\mu _{0j})`$ were found by formulas (64) and (67). Since the smoothed values do not lie within the limits of $`\pm \delta (\mu _{0j})`$ from $`\mu _j`$, we did not look for the most probable smoothness (see. Appendix), instead we considered the regularizer to be known and equal to the optimal one found in the method II A. The introduction of the prior information has significantly diminished the errors of $`\chi (k)`$ extraction (see dotted curves in Fig. 4) which were defined as $`\epsilon _j=\delta (\mu _{0j})/(\overline{\mu }_{0j}\mu _{bj})`$. This is quite natural: any decrease of our ignorance about $`\mu _0`$ should narrow the posterior distribution of $`\mu _{0j}`$ for all $`j`$. Of course, this concerns the experimental information as well: errors $`\epsilon _j`$ are the less the more measured points $`N`$ the spectrum has. Comparing Fig. 1( ) and Fig. 5(c), it is seen practically perfect coincidence of the results of bayesian smoothing and smoothing spline. From this one can assume the equality of the errors which both methods give.

Could we take into account possible systematic errors in the framework of the method? Yes, if we have the information on their nature and are able to translate it into the mathematics language; such a translation might be rather non-trivial. In any case, now we have the tool to extract from the prior and experimental information not only the sought values but their errors as well.

### D Other methods

Consider briefly the methods for $`\mu _0`$ construction not included into the viper program.

A rich variety of computer programs for XAS spectra processing is collected on the International XAFS Society Web-site . The vast majority of them use as an approximation for the atomic-like absorption a smoothing spline or more general piecewise-polynomial representation. For example, in the method of Ref. , the construction of $`\mu _0`$ is divided into several stages: $`\mu _0`$ is approximated by a low-degree polynomial, obtained $`\chi (k)`$ is multiplied by $`k^w`$, additional $`\mu _0^{}`$ is drawn again as a low-degree polynomial and subtracted, a smoothing spline then approximates one more additional $`\mu _0^{\prime \prime }`$. The sum of all $`\mu _0`$’s gives the total atomic-like background. The necessity of the preliminary stages was not discussed op. cit., however, clearly it was caused by the instability of spline with respect to the small variations of input parameters. And the point is not that the preliminary stages make the process stable, but that for each specific spectrum, auxiliary parameters (degrees of polynomials) could provide an acceptable construction of the atomic-like background. Above (in Sec. II A) we proposed the way to rise the stability of spline making the preliminary stages to be redundant.

In Ref.  an iterative approach to “atomic background” removal was developed. First a spline is used to obtain a rough estimate of the background; this alone is enough to have a reliable $`\chi `$ at $`k>56`$ Å-1. Over that range the $`\chi `$ obtained is fitted to the theoretical $`\chi _{\mathrm{th}}`$ in $`r`$-space. The resulting fit parameters are used then to generate $`\chi _{\mathrm{th}}(k)`$ that extends down to low $`k`$. This function is transformed back into $`e`$-space and $`\mu _0`$ is obtained as $`\mu _0=\mu /(\chi _{\mathrm{th}}+1)`$ that need be a little smoothed or fitted by an additional spline. Since the logic of reasoning was inversed: not “find $`\mu _0`$ to find $`\chi `$,” but “find $`\chi `$ to find $`\mu _0`$,” the method is suited for the quest of peculiarities on $`\mu _0`$ curve, not for structural XAFS-researches. Besides, the range of accuracy of the model appears to be unknown in principle: all, that is not described by the model, is included in $`\mu _0`$; the errors of the background approximation are also undefined.

In the old work  for the determination of the background absorption $`\mu _0`$ was considered the damping of the XAFS amplitude resulting from the measurements with low resolutions (with a large slit width). The superimposition of two spectra measured with different energy resolutions gives the intersection points, the part of which belong to the $`\mu _0`$. Then through the obtained nodal points a smoothing spline is drawn. As the authors of Ref.  noted, the measurements of the spectra with worsened resolution are not necessary; the spectra could be damped by the convolution with a “rocking curve,” approximated by a Gaussian function. Of course, the method is correctly works only with a small variation of the Gaussian curve width since for the large width not only the XAFS amplitude is damped but the very edge is washed out. Because of this only the extended part of a spectrum could be reliably determined.

The damping of the XAFS amplitude can be due to other reasons. For instance, as was pointed in Ref. , the nodal points may be obtained from the variable-temperature study. This idea was realized in Ref.  and is more sound since the atomic-like background is really independent of temperature and with temperature the XAFS amplitude is changed, not the shape of the edge. But for all that it is important that the phase difference between XAFS of different temperatures was negligible, which is true only for low wave numbers. Unfortunately, the method is suitable only for some particular cases (to say nothing of need for measured temperature series of spectra). Op. cit. it was demonstrated for the x-ray-absorption data for the $`L_3`$ edge of solid Pb. In those spectra the first crossing of $`\mu `$ and $`\mu _0`$ occurs already at $`15`$ eV above edge. In our sample spectra the first crossing occurs only at $`30`$ eV, which allows one to find at most 2–3 points and the first of them being situated at $`k2.5`$ Å-1.

An interesting approach to the problem of $`\mu _0`$ determination was reported in Ref. . It is based on the simple identity that relates the FT of some function with the FT of its $`n`$-th derivative:

$$\mathrm{FT}[f^{(n)}(k)]=(2ir)^n\mathrm{FT}[f(k)],$$
(6)
where the conjugate variables are $`k`$ and $`2r`$. Since the atomic-like background is smooth enough, the higher derivatives $`\mu ^{(n)}(k)`$ ($`n2`$) are oscillatory near zero. Performing the FT of $`\mu ^{(n)}(k)k^w`$ and using Eq. (6), one obtains the FT of unnormalized would-be $`\chi (k)k^w`$ (see Fig. 6). Op. cit. the low-$`r`$ part (which in our example is $`0r1.1`$) was cut off, and then the back FT was done. As a result, one has the unnormalized $`\chi (k)k^w`$ and, having subtracted it from the $`\mu (k)`$, the atomic-like background on which some peculiarities due to multi-electron excitations can be distinguished. Like the method of Ref. , this method is suited for the quest of peculiarities on the $`\mu _0`$ curve, not for structural XAFS-researches because of evident distortion of the first peak on the FT by the contribution from the atomic-like background. To illustrate this assertion, in Fig. 6 we show the FT of the second derivative of the $`\mu _0(k)`$ that was found by the present method. As seen, this contribution is not as small.

If the electronic states of an absorbing atom in gaseous phase and in the compound of interest may be considered as equivalent, $`\mu _0`$ can be set equal to the measured absorption in gas, as was done in Ref  for solid, liquid, and gaseous Kr. Some differences in energy positions and relative weights of double-electron excitation channels were taken into account by a model using simple empirical functions which were transferred then to the spectra of liquid and solid Kr. Notice that the proposed in the present paper prior function for the methods of smoothing spline and bayesian smoothing can include additional items corresponding to the multi-electron contributions.

## III Errors in $`\mu _0`$ constructing, noise, and choice of limits $`k_{\mathrm{min}}`$ $`k_{\mathrm{max}}`$

For what we need to know the errors of XAFS-function extraction? First, without knowing of these values one cannot in principle aim at their minimization. Second, they are used in the definition of $`\chi ^2`$-statistics in the fitting problems; their underestimation is a source of unjustified optimistic errors of fitting parameters. Third, along with analysis of the noise, the errors of $`\mu _0`$ construction allow us to choose the limits of reliable EXAFS signal, $`k_{\mathrm{min}}`$ and $`k_{\mathrm{max}}`$.

Unfortunately, the issue of quality of XAFS extraction from the measured absorption has not been addressed properly. We see several reasons for that. On the one hand, not having a correctly developed approach to the estimation of the errors of final results (interatomic distances, Debye-Waller factors etc. found via fitting), the errors of EXAFS extraction are useful. On the other hand, only a few methods include approaches to their estimations.

Easily one can compare the errors of different methods (see Fig. 4) and then choose the most reliable one. The problem of plausible limitations on the absolute value of the errors is more difficult. Define “signal” as the envelope of $`\chi (k)`$ (solid line with gray filling in Fig. 4). It is quite reasonable to demand that the errors of $`\mu _0`$ construction were less than XAFS signal. For the method of the interpolation spline drawn through the varied knots to meet this requirement leads to the restriction on the photoelectron wave numbers: $`2k14`$ Å-1. For the bayesian curve $`a`$ this range is $`0k14`$ Å-1, for the bayesian curves $`b`$ and $`c`$ this range is wider: $`0k16`$ Å-1.

Another factor that limitates the spectrum length is the presence of noise. To determine the noise is a straightforward task for $`r`$-space, where XAFS signals at high $`r`$ have clearly noise character. By Parseval’s identity the noise in $`r`$-space is related with the noise in $`k`$-space :

$`{\displaystyle _{k_{\mathrm{min}}}^{k_{\mathrm{max}}}}|n_kk^w|^2𝑑k=2{\displaystyle _0^{\pi /2dk}}|n_r|^2𝑑r.`$ (7)
Substitute the mean value over the range $`15<r<25`$ Å of the FT magnitude squared for $`|n_r|^2`$. Then

$`n_k^2=|n_r^2|{\displaystyle \frac{\pi }{dk}}{\displaystyle \frac{2w+1}{k_{\mathrm{max}}^{2w+1}k_{\mathrm{min}}^{2w+1}}}.`$ (8)
As seen from the formula, $`n_k`$ depends on $`dk`$, the size of evenly-spaced $`k`$-grid. Although above we already have used the Fourier transform, the question of choice of $`dk`$ was not raised yet. The algorithm of fast FT needs the transformed function to be set on a uniform grid. Having chosen a small $`dk`$, we artificially obtain the large number of “experimental” values. Naturally, this trick would not give more information than we have, and the errors $`n_k`$ must be large at the small $`dk`$. In our example the choice of $`dk`$ (0.03 Å-1) was based on the equality of numbers of experimental points and the nodes of the grid. The signal-to-noise ratio obtained is greater than unity for all the spectrum (see Fig. 4). There was no doubt in that: the signal is visually distinguished even for the very extended end of the spectrum (see Fig. 1(b) and Fig. 5(b)).

The noise can be estimated based on the bayesian considerations . Let after measurements we have the values of counts from the solid-state or gas-filled detectors and let there is a positive real number $`\lambda `$ such that the probability that a single count occurs in the time interval $`dt`$ is

$`P(1|\lambda )=\lambda dt.`$ (9)
It can be shown  that merely from this assumption follows that the counts obey the Poisson distribution law:

$`P(N|\lambda ,T)={\displaystyle \frac{(\lambda T)^N\mathrm{exp}(\lambda T)}{N!}},`$ (10)
where $`T`$ is the sampling time. The problem is to find the intensity $`\lambda `$ and its variance. Using Bayes theorem and introducing prior probabilities $`P(N)=1/N`$ and $`P(\lambda )=1/\lambda `$ , one obtains:

$`P(\lambda |N,T)={\displaystyle \frac{P(N|\lambda ,T)P(\lambda )}{P(N)}}={\displaystyle \frac{T(\lambda T)^{N1}\mathrm{exp}(\lambda T)}{(N1)!}},`$ (11)
that is after measurement the variate $`2T\lambda `$ follows the $`\chi ^2`$-distribution with $`2N`$ degrees of freedom. It is easy to find that $`\overline{\lambda }=N/T`$, $`\overline{\lambda ^2}=N(N+1)/T^2`$, and the variance of intensity is $`\delta \lambda =\sqrt{N}/T`$.

Denote counts from detectors measuring $`i_0`$ and $`i_1`$ as $`I_0`$ and $`I_1`$. By definition the variate $`\xi =\frac{i_0/2I_0}{i_1/2I_1}`$ follows Fisher’s $`F`$-distribution with $`(2I_0,2I_1)`$ degrees of freedom. Its expected value and variance are known: $`\overline{\xi }=I_1/(I_11)`$, $`\delta ^2\xi =I_1^2(I_0+I_11)/((I_11)^2(I_12)I_0)`$, from where we find for the absorption in the fluorescence mode ($`\mu x=i_0/i_1`$):

$`\overline{i_0/i_1}={\displaystyle \frac{I_0}{I_11}},\delta ^2(i_0/i_1)={\displaystyle \frac{I_0(I_0+I_11)}{(I_11)^2(I_12)}}.`$ (12)
Further, the variate $`\eta =\frac{1}{2}\mathrm{ln}\xi `$ follows $`z`$-distribution (Fisher’s distribution of variance ratio) with $`(2I_0,2I_1)`$ degrees of freedom. Its expected value and variance are known: $`\overline{\eta }=0`$, $`\delta ^2\eta =\frac{1}{4}(I_1+I_1)/(I_0I_1)`$, from where we find for the absorption in the transmission mode ($`\mu x=\mathrm{ln}(i_0/i_1)`$):

$`\overline{\mathrm{ln}(i_0/i_1)}=\mathrm{ln}{\displaystyle \frac{I_0}{I_1}},\delta ^2\mathrm{ln}(i_0/i_1)={\displaystyle \frac{1}{I_0}}+{\displaystyle \frac{1}{I_1}}.`$ (13)
The noise of XAFS-function is

$`n_P={\displaystyle \frac{\delta \mu }{\mu _0\mu _b}}=\left({\displaystyle \frac{1}{I_0}}+{\displaystyle \frac{1}{I_1}}\right)^{1/2}{\displaystyle \frac{1}{\mu _0\mu _b}}.`$ (14)
In our example this noise at $`k15`$ Å-1 becomes greater than signal (see Fig. 4). What is the reason for such significant difference between the really present noise $`n_k`$ and its statistical estimate $`n_P`$? Of course, the reason is in the false premise (9). In practice this condition is realized as: $`P(c|\lambda )=\lambda dt`$. For example, the photocurrent in an ion-chamber depends on gas pressure, potential applied etc.; these dependencies are contained in $`c`$. In other words, the amplification path works in such a way that one photon gives birth to $`c`$ counts. There is no difficulty in writing the posterior distribution for the generalized premise:

$`P(\lambda |N,T)={\displaystyle \frac{T(\lambda T)^{\frac{N}{c}1}\mathrm{exp}(\lambda T)}{\mathrm{\Gamma }(N/c)}},`$ (15)
with $`\overline{\lambda }=N/(cT)`$ and $`\delta \lambda =\sqrt{N/c}/T`$. Thus, having unknown $`c`$ (and implicitly assigning $`c=1`$), we got wrong variances for $`i_0/i_1`$ and $`\mathrm{ln}(i_0/i_1)`$. Unfortunately, in the most of real experiments the association between the probability of a single count event and the radiation intensity (via $`c`$) is unknown. In spite of this, the Poisson counting statistics is traditionally used for a long time. For example, in Ref.  signal-to-noise ratios are evaluated (assuming $`c=1`$) for the different detection schemes.

Practically all programs for XAFS spectra processing  to estimate the noise use the Fourier analysis. But then it is the noise that they use as uncertainties $`\epsilon _i`$ of $`\chi (k)`$ determination in definition of $`\chi ^2`$-statistics:

$`\chi ^2={\displaystyle \frac{N_{\mathrm{max}}}{M}}{\displaystyle \underset{i=1}{\overset{M}{}}}{\displaystyle \frac{[(\chi _{\mathrm{exp}})_i(\chi _{\mathrm{mod}})_i]^2}{\epsilon _i^2}}.`$ (16)
It would be more correct to consider as $`\epsilon _i`$ the larger from the two: the noise and the errors of the construction of $`\mu _0`$. In our case (and as a rule) the latter are essentially greater (especially in the method II B) than the noise. In the following paper  we shall show how the understated $`\epsilon _i`$ lead to optimistic errors of structural parameters.

## IV XAFS-function correction

Because of one reason or another the experimental XAFS might be distorted. Consider some of them.

(i) Let the counts ($`I`$) from detectors are associated with the intensities ($`i`$) as $`i_0=\varkappa _0I_0`$ and $`i_1=\varkappa _1I_1`$. Then the absorption (in the transmission mode) equals:

$`\mu x=\mathrm{ln}(i_0/i_1)=\mathrm{ln}(I_0/I_1)+\mathrm{ln}(\varkappa _0/\varkappa _1).`$ (17)
The second term is a slightly varied function of energy and can be taken into account in independent experiments. Such a distortion appears, for instance, if the absorptance of the gas in ion-chamber detectors depends on energy.

(ii) If some part of incident radiation is not attenuated in the sample as much as expected (due to the pinholes in the sample, harmonics in the incoming beam etc.), that is $`i_0=\varkappa _0I_0+b`$, then the real absorption is connected with the measured $`I_0`$ and $`I_1`$ in a complicated way. In Ref.  the possible decrease of XAFS amplitude shown to be essential even at low $`b/(\varkappa _0I_0)`$ but thick samples. At known ratio $`b/(\varkappa _0I_0)`$, the correcting factor can be easily obtained.

(iii) In the fluorescence mode, due to absorptance of the fluorescent signal in the sample itself XAFS spectra strongly depend on the detection geometry. In Ref.  the correcting functions are found explicitly.

(iv) The problem of glitches is widespread in the XAFS analysis. The glitches are due to multiple Bragg reflection being satisfied simultaneously and for each given monochromator are manifested in the strictly determined spectral positions. In most cases the glitches seen on curves $`I_0(E)`$ and $`I_1(E)`$, vanish on $`I_0/I_1`$ ratio. If not, one can easily get rid of them. For instance, the glitch area, usually extremely thin, is smoothed or, with fixed ends, replaced by a straight-line segment. The main thing in the correct analysis of glitches is their detection.

To detect a glitch on curves $`\mu `$ or $`\chi k^w`$ is practically impossible. For this, one needs the primary data $`I_0(E)`$ and $`I_1(E)`$, not $`\mathrm{ln}(I_0/I_1)`$ nor $`I_0/I_1`$. Out of glitches the intensity of incident radiation smoothly, ignoring the noise, depends on energy (see Fig. 7). The idea of detection of glitches via critical level for the derivative $`|d\mathrm{ln}I_0/dE|_c`$ is self-evident. For the presented in Fig. 7 $`I_0(E)`$ curve the absolute value of the derivative in the glitch areas is greater than the critical value chosen to be equal to $`1.7710^3`$-1. Having extracted the XAFS, one can see that the first (paired) glitch $`a`$ is not manifested on $`\chi (k)k^2`$, the last two ($`c`$ $`d`$) are obscured by the noise, solely glitch $`b`$ is clearly pronounced. Now, being in the firm belief that this is not a part of the XAFS, one can eliminate the glitch with ease. Here, we fixed its ends on $`\mu (E)`$, replaced it by the straight-line segment, and constructed $`\chi (k)k^2`$ again.

## V Conclusion

In this paper we have considered all stages of XAFS function extraction from the measured absorption. We focused our attention on the most important stage, construction of the atomic-like absorption $`\mu _0`$.

For the wide-spread method of approximation of $`\mu _0`$ by a smoothing spline we have proposed the way to raise the stability by including the prior information about absorption edge shape (“nearly step” or “nearly step with a white line”). Besides we have propose a new reliable criterion for determination of the smoothing regularizer.

A new method for approximation of $`\mu _0`$ is proposed, the method of bayesian smoothing. It can include various prior information, which raises the accuracy of XAFS determination. Following this method one finds the distributions of $`\mu _0`$ in each experimental point, from which one can find not only average values but also any desirable momenta, which appears to be an additional difficulty for other methods. This method was shown to give more accurate atomic-like background than that obtained by the method of Ref .

Particular attention has been given to the analysis of noise. We have discussed the difficulties of its estimates on the basis of statistical approach. More reliable is the determination of noise from the Fourier transform. We have shown that the experimental noise is essentially less than the errors of $`\mu _0`$ construction, and the use of values of noise in the $`\chi ^2`$-statistics definition appears to be erroneous since leads to the unjustified optimistic errors of structural parameters inferred in fitting procedures. For detailed consideration of the accuracy of fitting parameters see the following paper .

## Bayesian smoothing and deconvolution

### 1 Posterior distribution for smoothed data

Consider general linear problem of data smoothing with the use of statistical methods (for introduction see review by Turchin et al.  and the articles from Web-site bayes.wustl.edu). Let data $`𝐝`$ are defined on the mesh $`x_1,\mathrm{},x_N`$ and consist of the true values $`𝐭`$ and the additive noise $`𝐧`$:

$$d_i=t_i+n_i,i=1,\mathrm{},N.$$
(18)
The problem of smoothing is to find the best estimates of $`𝐭`$. For an arbitrary node $`j`$, find the probability density function for $`t_j`$ given the data $`𝐝`$:

$$P(t_j|𝐝)=\mathrm{}𝑑t_{ij}\mathrm{}P(𝐭|𝐝),$$
(19)
where $`P(𝐭|𝐝)`$ is the joint probability density function for all values $`𝐭`$, and the integration is done over all $`t_{ij}`$. According to Bayes theorem,

$$P(𝐭|𝐝)=\frac{P(𝐝|𝐭)P(𝐭)}{P(𝐝)},$$
(20)
$`P(𝐭)`$ being the joint prior probability for all $`t_i`$, $`P(𝐝)`$ is a normalization constant. Assuming that the values $`n_i`$ are independent in different nodes and normally distributed with zero expected values, the probability $`P(𝐝|𝐭)`$, so-called likelihood function, is given by

$$P(𝐝|𝐭,\sigma )=(2\pi \sigma ^2)^{N/2}\mathrm{exp}\left(\frac{1}{2\sigma ^2}\underset{k=1}{\overset{N}{}}(d_kt_k)^2\right),$$
(21)
where the standard deviation of the noise, $`\sigma `$, appears as a known value. Later, we apply the rules of probability theory to remove $`\sigma `$ from the problem.

Now define prior probability $`P(𝐭)`$. Let we know in advance that the function $`t(x)`$ is smooth enough. To specify this information, introduce the norm of the second derivative and indicate its expected approximate value:

$$\mathrm{\Omega }(t(x))=\left(\frac{d^2t}{dx^2}\right)^2𝑑x\omega .$$
(22)
Denote $`\mathrm{\Delta }_i=x_{i+1}x_i`$, $`i=1,\mathrm{},N1`$ and represent the second derivative in the finit-difference form:

$$\mathrm{\Omega }(t(x))\mathrm{\Omega }(𝐭)=\underset{i=2}{\overset{N1}{}}[t_{i1}\mathrm{\Delta }_{i1}^1t_i(\mathrm{\Delta }_{i1}^1+\mathrm{\Delta }_i^1)+t_{i+1}\mathrm{\Delta }_i^1]^2\underset{k,l=1}{\overset{N}{}}\mathrm{\Omega }_{kl}t_kt_l.$$
(23)
$`\mathrm{\Omega }_{kl}`$ is a five-diagonal symmetric matrix with the following non-zero elements:

$`\mathrm{\Omega }_{11}`$ $`=`$ $`\mathrm{\Delta }_1^1\mathrm{\Delta }_2^2,\mathrm{\Omega }_{22}=\mathrm{\Delta }_2^1(\mathrm{\Delta }_1^1+\mathrm{\Delta }_2^1)^2+\mathrm{\Delta }_2^2\mathrm{\Delta }_3^1,\mathrm{\Omega }_{12}=(\mathrm{\Delta }_1\mathrm{\Delta }_2)^1(\mathrm{\Delta }_1^1+\mathrm{\Delta }_2^1),`$ (25)
$`\mathrm{}\mathrm{}\mathrm{}\mathrm{}`$
$`\mathrm{\Omega }_{ii}`$ $`=`$ $`\mathrm{\Delta }_i^1(\mathrm{\Delta }_i^1+\mathrm{\Delta }_{i1}^1)^2+\mathrm{\Delta }_i^2\mathrm{\Delta }_{i+1}^1+\mathrm{\Delta }_{i1}^3,`$ (26)
$`\mathrm{\Omega }_{i1,i}`$ $`=`$ $`\mathrm{\Delta }_{i1}^2(\mathrm{\Delta }_{i1}^1+\mathrm{\Delta }_{i2}^1)(\mathrm{\Delta }_{i1}\mathrm{\Delta }_i)^1(\mathrm{\Delta }_{i1}^1+\mathrm{\Delta }_i^1),`$ (27)
$`\mathrm{\Omega }_{i2,i}`$ $`=`$ $`\mathrm{\Delta }_{i2}^1\mathrm{\Delta }_{i1}^2,`$ (29)
$`\mathrm{}\mathrm{}\mathrm{}\mathrm{}`$
$`\mathrm{\Omega }_{NN}`$ $`=`$ $`\mathrm{\Delta }_{N1}^3,\mathrm{\Omega }_{N1,N1}=\mathrm{\Delta }_{N1}^1(\mathrm{\Delta }_{N1}^1+\mathrm{\Delta }_{N2}^1)^2+\mathrm{\Delta }_{N2}^3,\mathrm{\Omega }_{N1,N}=\mathrm{\Delta }_{N1}^2(\mathrm{\Delta }_{N1}^1+\mathrm{\Delta }_{N2}^1).`$ (30)
In order to introduce the minimum information in addition to that contained in (23), from all normalized to unity functions $`P(𝐭)`$ which satisfy the condition (23) we choose a single one that contains minimum information about t that is minimizes the functional

$$I[P(𝐭)]=P(𝐭)\mathrm{ln}P(𝐭)𝑑𝐭+\beta \left[1P(𝐭)𝑑𝐭\right]+\gamma \left[\omega \mathrm{\Omega }(𝐭)P(𝐭)𝑑𝐭\right],$$
(31)
where $`\beta `$ and $`\gamma `$ are the Larrange multipliers. In minimizing $`I[P(𝐭)]`$, one obtains the equation set

$`\mathrm{ln}P(𝐭)+1\beta \gamma \mathrm{\Omega }(𝐭)`$ $`=`$ $`0`$ (32)
$`{\displaystyle P(𝐭)𝑑𝐭}`$ $`=`$ $`1`$ (33)
$`{\displaystyle \mathrm{\Omega }(𝐭)P(𝐭)𝑑𝐭}`$ $`=`$ $`\omega ,`$ (34)
that has a solution:

$$P(𝐭)=(\lambda _1\mathrm{}\lambda _N)^{1/2}\left(\frac{2\pi \sigma ^2}{\alpha }\right)^{N/2}\mathrm{exp}\left(\frac{\alpha }{2\sigma ^2}\mathrm{\Omega }(𝐭)\right),$$
(35)
where $`\alpha /2\sigma ^2=\gamma =N/2\omega `$, and $`\lambda _1,\mathrm{},\lambda _N`$ are the eigenvalues of the matrix $`\mathrm{\Omega }_{kl}`$. The regularizer $`\alpha `$ will be used to control the smoothness of $`𝐭`$. The prior distribution obtained is a “soft” one, that is does not demand from the solution to have a strictly prescribed form.

Thus, we have for the probability density function:

$`P(t_j|𝐝,\sigma ,\alpha )`$ $``$ $`{\displaystyle \mathrm{}𝑑t_{ij}\mathrm{}\sigma ^{2N}\alpha ^{N/2}\mathrm{exp}\left(\frac{\alpha }{2\sigma ^2}\underset{k,l=1}{\overset{N}{}}\mathrm{\Omega }_{kl}t_kt_l\right)\mathrm{exp}\left(\frac{1}{2\sigma ^2}\underset{k=1}{\overset{N}{}}(d_kt_k)^2\right)}`$ (36)
$`=`$ $`{\displaystyle \mathrm{}𝑑t_{ij}\mathrm{}\sigma ^{2N}\alpha ^{N/2}\mathrm{exp}\left(\frac{1}{2\sigma ^2}\left[𝐝^22\underset{k=1}{\overset{N}{}}d_kt_k+\underset{k,l=1}{\overset{N}{}}g_{kl}t_kt_l\right]\right)},`$ (37)
where

$`g_{kl}=\alpha \mathrm{\Omega }_{kl}+\delta _{kl},𝐝^2={\displaystyle \underset{k=1}{\overset{N}{}}}d_k^2.`$ (38)
Since there is no integral over $`t_j`$, separate it from the other integration variables:

$`P(t_j|𝐝,\sigma ,\alpha )`$ $``$ $`\sigma ^{2N}\alpha ^{N/2}\mathrm{exp}\left({\displaystyle \frac{1}{2\sigma ^2}}[𝐝^22d_jt_j+g_{jj}t_j^2]\right)`$ (40)
$`\times {\displaystyle }\mathrm{}dt_{ij}\mathrm{}\mathrm{exp}({\displaystyle \frac{1}{2\sigma ^2}}[\underset{k,l=1}{\overset{N}{{\displaystyle }^j}}g_{kl}t_kt_l2\underset{k=1}{\overset{N}{{\displaystyle }^j}}[d_kg_{kj}t_j]t_k]),`$
Here, the symbol $`j`$ near the summation signs denotes the absence of $`j`$-th item. Further, find the eigenvalues $`\lambda _i^{}`$ and corresponding eigenvectors $`𝐞_i`$ of the matrix $`g_{kl}`$ in which the $`j`$-th row and column are deleted, and change the variables:

$$b_i=\sqrt{\lambda _i^{}}\underset{k=1}{\overset{N}{^j}}t_ke_{ik},t_k=\underset{i=1}{\overset{N}{^j}}\frac{b_ie_{ik}}{\sqrt{\lambda _i^{}}}(i,kj).$$
(41)
Using the properties of eigenvectors:

$$\underset{k=1}{\overset{N}{^j}}g_{lk}e_{ik}=\lambda _i^{}e_{il},\underset{k=1}{\overset{N}{^j}}e_{lk}e_{ik}=\delta _{li}(l,ij),$$
(42)
one obtains:

$`P(t_j|𝐝,\sigma ,\alpha )`$ $``$ $`\sigma ^{2N}\alpha ^{N/2}\mathrm{exp}\left({\displaystyle \frac{1}{2\sigma ^2}}[(𝐝^2𝐡^2)2t_j(d_j\mathrm{𝐡𝐮})+t_j^2(g_{jj}𝐮^2)]\right)`$ (44)
$`\times {\displaystyle }\mathrm{}db_{lj}\mathrm{}\mathrm{exp}({\displaystyle \frac{1}{2\sigma ^2}}\underset{i=1}{\overset{N}{{\displaystyle }^j}}[b_ih_i+u_it_j]^2),`$
where new quantities were introduced:

$`h_i={\displaystyle \frac{1}{\sqrt{\lambda _i^{}}}}\underset{k=1}{\overset{N}{{\displaystyle }^j}}d_ke_{ik},u_i={\displaystyle \frac{1}{\sqrt{\lambda _i^{}}}}\underset{k=1}{\overset{N}{{\displaystyle }^j}}g_{kj}e_{ik},`$ (45)
$`𝐡^2=\underset{i=1}{\overset{N}{{\displaystyle }^j}}h_i^2,𝐮^2=\underset{i=1}{\overset{N}{{\displaystyle }^j}}u_i^2,\mathrm{𝐡𝐮}=\underset{i=1}{\overset{N}{{\displaystyle }^j}}h_iu_i.`$ (46)
Evaluating the $`N1`$ integrals in (44), one finally obtains the posterior probability for $`j`$-th node:

$$P(t_j|𝐝,\sigma ,\alpha )\sigma ^{(N+1)}\alpha ^{N/2}\mathrm{exp}\left(\frac{1}{2\sigma ^2}[(𝐝^2𝐡^2)2t_j(d_j\mathrm{𝐡𝐮})+t_j^2(g_{jj}𝐮^2)]\right).$$
(47)
### 2 Eliminating nuisance parameters

In most real problems $`\sigma `$ and $`\alpha `$ are not known. To eliminate $`\sigma `$ is a quite straightforward problem:

$$P(t_j|𝐝,\alpha )=𝑑\sigma P(t_j,\sigma |𝐝,\alpha )=𝑑\sigma P(\sigma )P(t_j|𝐝,\sigma ,\alpha ),$$
(48)
one needs only to know a prior probability $`P(\sigma )`$. Having no specific information about $`\sigma `$, a Jeffreys prior $`P(\sigma )=1/\sigma `$ is assigned . Then

$`P(t_j|𝐝,\alpha )`$ $``$ $`{\displaystyle _0^{\mathrm{}}}𝑑\sigma \sigma ^{(N+2)}\mathrm{exp}\left({\displaystyle \frac{1}{2\sigma ^2}}[(𝐝^2𝐡^2)2t_j(d_j\mathrm{𝐡𝐮})+t_j^2(g_{jj}𝐮^2)]\right)`$ (49)
$``$ $`[(𝐝^2𝐡^2)2t_j(d_j\mathrm{𝐡𝐮})+t_j^2(g_{jj}𝐮^2)]^{(N+1)/2}.`$ (50)
Introducing the substitution

$`w_j^2=N{\displaystyle \frac{(g_{jj}𝐮^2)^2}{(𝐝^2𝐡^2)(g_{jj}𝐮^2)(d_j\mathrm{𝐡𝐮})^2}}\left(t_j{\displaystyle \frac{d_j\mathrm{𝐡𝐮}}{g_{jj}𝐮^2}}\right)^2,`$ (51)
one obtains the Student $`t`$-distribution with $`N`$ degrees of freedom:

$`P(w_j|𝐝,\alpha )\left(1+{\displaystyle \frac{w_j^2}{N}}\right)^{(N+1)/2}`$ (52)
with zero average and the variance $`N/(N2)`$. From where one finds for $`t_j`$:

$`\overline{t}_j={\displaystyle \frac{d_j\mathrm{𝐡𝐮}}{g_{jj}𝐮^2}},\delta ^2(t_j)={\displaystyle \frac{(𝐝^2𝐡^2)(g_{jj}𝐮^2)(d_j\mathrm{𝐡𝐮})^2}{(g_{jj}𝐮^2)^2}}{\displaystyle \frac{1}{N2}}.`$ (53)
Thus, we have got rid of unknown $`\sigma `$ and found the expressions for mean values $`t_j`$ and their dispersions at known regularizer $`\alpha `$. To eliminate the latter is more difficult. The idea is not to find the smoothest solution, but the solution of the most probable smoothness. For that we will find the posterior probability:

$`P(\alpha |𝐝)={\displaystyle 𝑑𝐭𝑑\sigma P(\alpha ,\sigma ,𝐭|𝐝)}={\displaystyle 𝑑𝐭𝑑\sigma P(\alpha ,\sigma )P(𝐭|\alpha ,\sigma ,𝐝)}.`$ (54)
Assuming that $`\alpha `$ and $`\sigma `$ are independent and using Bayes theorem (20), one obtains:

$`P(\alpha |𝐝){\displaystyle 𝑑𝐭𝑑\sigma P(\alpha )P(\sigma )P(𝐭|\alpha ,\sigma )P(𝐝|𝐭,\alpha ,\sigma )}.`$ (55)
Substituting (35) for the prior probability $`P(𝐭|\alpha ,\sigma )`$, (21) for the likelihood, and a Jeffreys prior $`P(\sigma )=1/\sigma `$ and $`P(\alpha )=1/\alpha `$, one obtains the posterior distribution for the regularizer $`\alpha `$:

$`P(\alpha |𝐝)`$ $``$ $`{\displaystyle 𝑑𝐭𝑑\sigma \sigma ^{2N1}\alpha ^{N/21}\mathrm{exp}\left(\frac{\alpha }{2\sigma ^2}\underset{k,l=1}{\overset{N}{}}\mathrm{\Omega }_{kl}t_kt_l\right)\mathrm{exp}\left(\frac{1}{2\sigma ^2}\underset{k=1}{\overset{N}{}}(d_kt_k)^2\right)}`$ (56)
$`=`$ $`{\displaystyle 𝑑𝐭𝑑\sigma \sigma ^{2N1}\alpha ^{N/21}\mathrm{exp}\left(\frac{1}{2\sigma ^2}\left[𝐝^22\underset{k=1}{\overset{N}{}}d_kt_k+\underset{k,l=1}{\overset{N}{}}g_{kl}t_kt_l\right]\right)},`$ (57)
where matrix $`g_{kl}`$ was defined in (38). After its diagonalization, analogously to what was done above, finally one obtains:

$`P(\alpha |𝐝)(\lambda _1^{}\mathrm{}\lambda _N^{})^{1/2}\alpha ^{N/21}[𝐝^2𝐡^2]^{N/2},`$ (58)
where $`𝐡^2`$ is given by

$`𝐡^2={\displaystyle \underset{i=1}{\overset{N}{}}}h_i^2,h_i={\displaystyle \frac{1}{\sqrt{\lambda _i^{}}}}{\displaystyle \underset{k=1}{\overset{N}{}}}d_ke_{ik},`$ (59)
and $`\lambda _i^{}`$ and $`𝐞_i`$ are the eigenvalues and eigenvectors of $`g_{kl}`$. Having found the maximum of the posterior probability (58) or having averaged over it the expression (53), one has the sought $`𝐭`$ with the most probable smoothness. However it is necessary to point out that this procedure narrows the applicability of the bayesian smoothing down to the class of tasks where the smoothed values lie in most within the limits $`\pm \sigma `$ from the most probable. In practice, there possible other tasks where the condition (18) is treated more wider and the smoothed values exceed the bounds of noise.

### 3 Expressions for smoothed values and their variances

The formulas (53) appear useless in practice since require to find the eigenvalues and eigenvectors for the matrix of rank $`N1`$ on each node. Those formulas have merely a methodological value: the explicit expressions for posterior probabilities enable one to find the average of arbitrary function of $`t_j`$. However, $`\overline{t}_j`$ and $`\delta ^2(t_j)`$ could be found significantly easier. Using (48) and (36), represent $`\overline{t}_j`$ as:

$`P(\sigma )P(t_j|𝐝,\sigma ,\alpha )dt_j`$ (60)
$``$ $`{\displaystyle 𝑑𝐭𝑑\sigma \sigma ^{2N1}t_j\mathrm{exp}\left(\frac{1}{2\sigma ^2}\left[𝐝^22\underset{k=1}{\overset{N}{}}d_kt_k+\underset{k,l=1}{\overset{N}{}}g_{kl}t_kt_l\right]\right)}.`$ (61)
Performing the diagonalization, one obtains:

$`\overline{t}_j`$ $``$ $`{\displaystyle 𝑑𝐛𝑑\sigma \sigma ^{2N1}\mathrm{exp}\left(\frac{1}{2\sigma ^2}[𝐝^2𝐡^2]\right)\left(\underset{i=1}{\overset{N}{}}\frac{b_ie_{ij}}{\sqrt{\lambda _i^{}}}\right)\mathrm{exp}\left(\frac{1}{2\sigma ^2}\underset{i=1}{\overset{N}{}}[b_ih_i]^2\right)}`$ (62)
$``$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{h_ie_{ij}}{\sqrt{\lambda _i^{}}}}{\displaystyle 𝑑\sigma \sigma ^{N1}\mathrm{exp}\left(\frac{1}{2\sigma ^2}[𝐝^2𝐡^2]\right)},`$ (63)
from where

$`\overline{t}_j={\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{h_ie_{ij}}{\sqrt{\lambda _i^{}}}}.`$ (64)
Analogously, for the variance $`\delta (t_j)`$ one has:

$`\delta ^2(t_j)`$ $``$ $`{\displaystyle 𝑑𝐛𝑑\sigma \sigma ^{2N1}\mathrm{exp}\left(\frac{1}{2\sigma ^2}[𝐝^2𝐡^2]\right)\left(\underset{i=1}{\overset{N}{}}\frac{(b_ih_i)e_{ij}}{\sqrt{\lambda _i^{}}}\right)^2\mathrm{exp}\left(\frac{1}{2\sigma ^2}\underset{i=1}{\overset{N}{}}[b_ih_i]^2\right)}`$ (65)
$``$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{e_{ij}^2}{\lambda _i^{}}}{\displaystyle 𝑑\sigma \sigma ^{N1}\mathrm{exp}\left(\frac{1}{2\sigma ^2}[𝐝^2𝐡^2]\right)\sigma ^2}.`$ (66)
Normalizing, one finally obtains:

$`\delta ^2(t_j)`$ $`=`$ $`{\displaystyle \frac{𝑑\sigma \sigma ^{N+1}\mathrm{exp}\left([𝐝^2𝐡^2]/2\sigma ^2\right)}{𝑑\sigma \sigma ^{N1}\mathrm{exp}\left([𝐝^2𝐡^2]/2\sigma ^2\right)}}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{e_{ij}^2}{\lambda _i^{}}}`$ (67)
$`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{N}{2}1\right)}{\left([𝐝^2𝐡^2]/2\right)^{N/21}}}{\displaystyle \frac{\left([𝐝^2𝐡^2]/2\right)^{N/2}}{\mathrm{\Gamma }\left(\frac{N}{2}\right)}}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{e_{ij}^2}{\lambda _i^{}}}={\displaystyle \frac{[𝐝^2𝐡^2]}{N2}}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{e_{ij}^2}{\lambda _i^{}}}.`$ (68)
Now we got the usable formulas, which require to find the eigenvalues and eigenvectors for the matrix of rank $`N`$ just one time.

### 4 Addenda to the bayesian smoothing

(i) Let the curvature of the function $`t(x)`$ is approximately known in advance. To specify this information, introduce the norm of the difference between $`d^2t/dx^2`$ and approximately known second derivative $`d^2f/dx^2`$:

$$\mathrm{\Omega }(t(x))=\left(\frac{d^2t}{dx^2}\frac{d^2f}{dx^2}\right)^2𝑑x\omega .$$
(69)
Notice, that there is no need to know $`f(x)`$ itself, its second derivative is sufficient. The explicit presence of $`f(x)`$ in the following formulas should be taken as a consequence of the technical trick applied: at first $`f(x)`$ is subtracted from the data, then it is added to the found solution.

Everywhere in formulas (2367) make the substitutions:

$$\stackrel{~}{t}_i=t_if_i,\stackrel{~}{d}_i=d_if_i,i=1,\mathrm{},N.$$
(70)
Performing the described above procedure for smoothing, one finds $`\stackrel{~}{t}_i`$, from which by inverse transformation the sought vector is given by $`𝐭=\stackrel{~}{𝐭}+𝐟`$.

(ii) In some tasks the value on the starting (zero) node is known without measurement. This sort of prior information represents a “hard” one, that is it restricted the class of possible solutions; in the given case the solution must pass through the known zero node. The quadratic form $`\mathrm{\Omega }(𝐭)`$ (or $`\mathrm{\Omega }(\stackrel{~}{𝐭})`$ in the case of approximately known second derivative) in the expression for the prior probability has changed:

$$\mathrm{\Omega }(𝐭)=\underset{i=1}{\overset{N1}{}}[t_{i1}\mathrm{\Delta }_{i1}^1t_i(\mathrm{\Delta }_{i1}^1+\mathrm{\Delta }_i^1)+t_{i+1}\mathrm{\Delta }_i^1]^2\underset{k,l=1}{\overset{N}{}}\mathrm{\Omega }_{kl}t_kt_l+\mathrm{\Omega }_{00}t_0^2+2\mathrm{\Omega }_{01}t_0t_1+2\mathrm{\Omega }_{02}t_0t_2,$$
(71)
the first few matrix elements of $`\mathrm{\Omega }_{kl}`$ now are:

$`\mathrm{\Omega }_{00}`$ $`=`$ $`\mathrm{\Delta }_0^2\mathrm{\Delta }_1^1,\mathrm{\Omega }_{01}=(\mathrm{\Delta }_0\mathrm{\Delta }_1)^1(\mathrm{\Delta }_0^1+\mathrm{\Delta }_1^1),\mathrm{\Omega }_{02}=\mathrm{\Delta }_0^1\mathrm{\Delta }_1^2,`$ (72)
$`\mathrm{\Omega }_{11}`$ $`=`$ $`\mathrm{\Delta }_1^1(\mathrm{\Delta }_0^1+\mathrm{\Delta }_1^1)^2+\mathrm{\Delta }_1^1\mathrm{\Delta }_2^2,\mathrm{\Omega }_{12}=\mathrm{\Delta }_1^2(\mathrm{\Delta }_1^1+\mathrm{\Delta }_0^1)(\mathrm{\Delta }_1\mathrm{\Delta }_2)^1(\mathrm{\Delta }_1^1+\mathrm{\Delta }_2^1).`$ (73)
If $`t_0=0`$ (or $`\stackrel{~}{t}_0=0`$), none further changes to the formulas of smoothing (2367) are needed; at $`t_00`$ the changes are evident: instead of the scalar product $`\mathrm{𝐝𝐭}`$ in (36) will be $`(𝐝\widehat{𝐝})𝐭`$, where $`\widehat{d}_1=\alpha t_0\mathrm{\Omega }_{01}`$, $`\widehat{d}_2=\alpha t_0\mathrm{\Omega }_{02}`$, all remaining $`\widehat{d}_i=0`$; to the $`𝐝^2`$ the term $`\alpha t_0^2\mathrm{\Omega }_{00}`$ will be added.

(iii) Making some changes in the considered above problem of smoothing allows one to solve the problem of deconvolution. If the experimental value $`d_j`$ on some node $`j`$ is determined not only by $`t_j`$ but also by the values of some neighboring nodes, then instead of (18) we have:

$`d_i={\displaystyle \underset{j=1}{\overset{N}{}}}r_{ij}t_j+n_i,i=1,\mathrm{},N,`$ (74)
where $`r_{ij}`$ is the grid representation of the impulse response function. Instead of expression (21), for the likelihood now we have:

$$P(𝐝|𝐭,\sigma )=(2\pi \sigma ^2)^{N/2}\mathrm{exp}\left(\frac{1}{2\sigma ^2}\underset{k=1}{\overset{N}{}}\left[d_k\underset{i=1}{\overset{N}{}}r_{ik}t_i\right]^2\right),$$
(75)
and instead of (36), the posterior probability for $`t_j`$ is now expressed as:

$`P(t_j|𝐝,\sigma ,\alpha ){\displaystyle \mathrm{}𝑑t_{ij}\mathrm{}\sigma ^{2N}\alpha ^{N/2}\mathrm{exp}\left(\frac{1}{2\sigma ^2}\left[𝐝^22\underset{k=1}{\overset{N}{}}D_kt_k+\underset{k,l=1}{\overset{N}{}}G_{kl}t_kt_l\right]\right)},`$ (76)
where

$`G_{kl}=\alpha \mathrm{\Omega }_{kl}+{\displaystyle \underset{i=1}{\overset{N}{}}}r_{ik}r_{il},D_k={\displaystyle \underset{i=1}{\overset{N}{}}}r_{ik}d_i.`$ (77)
Further steps for finding of $`𝐭`$ are analogous to the described above.


format: text

text:
# Ground state cooling, quantum state engineering and study of decoherence of ions in Paul traps

## 1 Quantum information processing with trapped ions

Trapped and laser cooled ions in Paul traps are used for an implementation of quantum information processing. Internal electronic states of individual ions serve to hold the quantum information (qubits) and an excitation of common vibrational modes provides the coupling between qubits, which is necessary for quantum logic operations between qubits, more specifically, for the realization of gate operations between two ions. The Cirac-Zoller proposal  requires that initially the ions are cooled to the ground state of motion and that the whole system can be coherently manipulated and controlled. The time scale of decoherence and coupling to the environment is required to be much smaller than the time scale of coherent manipulation.

The paper is organized as follows: In section II, we outline the techniques which are used for storing and detecting single ions in Paul traps. Section III decribes high resolution laser spectroscopy on the S$`{}_{1/2}{}^{}`$ D5/2 transition, which is the basis of coherent manipulation of the internal electronic and external vibrational quantum state of a trapped ion. Ground state cooling for one and two ions is reported in section IV. For a single ion, we reached 99.9 $`\%`$ of motional ground state occupation. Most recently, to scale up our results from one ion towards strings of ions, we have also shown ground state cooling for two ions. After a single ion is has been cooled, it is prepared for coherent manipulation (section V) of the ion’s vibrational state, including the generation of Fock states. We measured the coherence time of this process as well as the time scale of motional heating. Our experimental finding of 1 ms coherence time, and a heating rate of only one phonon in 190 ms, shows good conditions for quantum information processing. Finally, in section VI, we investigate the maximum speed of gate operations and find that with our system approximately 30 - 50 gate operations are possible within the coherence time.

## 2 Paul traps for single Calcium ions

Ions are stored in Paul traps  under ultra high vacuum ($`10^{\mathrm{10..11}}`$ mbar) conditions. In one of the experiments we use a spherical quadrupole Paul trap with a 1.4 mm ring diameter (Fig. 1.a) , where we observe motional frequencies ($`\omega _x`$, $`\omega _y`$, $`\omega _z`$)/(2$`\pi `$) of up to (2.16, 2.07, 4.51) MHz along the respective trap axes. The axial direction of the trap is denoted by $`z`$, the degeneracy of radial directions $`x`$ and $`y`$ is lifted by small asymmetries of the setup. Strings of ions are held in a linear trap (Fig. 1.b) , where we typically reach motional frequencies of ($`\omega _x`$, $`\omega _y`$, $`\omega _z`$)/(2$`\pi `$)= (2.0, 2.0, 0.7) MHz, e.g. three ions crystallize at an axial distance of 7$`\mu `$m.

40Ca+ ions have a single valence electron and no hyperfine structure. All relevant transitions are accessible by solid state or diode lasers (see Fig. 1.c) . In our experiment, we apply Doppler cooling to the ion and detect the internal state on the S1/2 to P1/2 transition at 397 nm, excited with a frequency-doubled Ti:Sapphire laser. This transition has a natural linewidth of 20 MHz and is not closed since the ion may decay to the metastable D3/2 level. A diode laser at 866 nm serves to repump the ion via the P1/2 state thus closing the cooling cycle. As the upper internal level for quantum state engineering and sub-Doppler cooling, we employ the metastable D5/2 level with a natural lifetime of approx. 1 s. The S$`{}_{1/2}{}^{}`$ D5/2 quadrupole transition at 729 nm is excited with a Ti:Sapphire laser. We can detect whether a transition to D5/2 occurred by applying the beams at 397 nm and 866 nm and monitoring the fluorescence of the ion (electron shelving technique). The internal state of the ion is discriminated with an efficiency close to 100$`\%`$ . Another diode laser at 854 nm is used to repump the ion from the D5/2 level to the electronic ground state via the P3/2 level. We observe the ions’ fluorescence on a photomultiplier and an intensified CCD camera.

## 3 Spectroscopy

For coherent spectroscopic investigation and state engineering on the narrow S$`{}_{1/2}{}^{}`$ D5/2 transition at 729 nm we use a pulsed technique which consists of five consecutive steps.

1. Laser light at 397 nm, 866 nm, and 854 nm is used to pump the ion to the S1/2 ground state. At the Doppler limit, $`E=\mathrm{}\mathrm{\Gamma }/2`$ , the thermal vibrational state corresponds to a mean vibrational quantum number $`n`$ 10 for $`\omega =(2\pi )\mathrm{\hspace{0.25em}1}`$ MHz.
2. A constant magnetic field of a few Gauss splits the 10 Zeeman components of the S$`{}_{1/2}{}^{}`$ D5/2 transition in frequency space. The S$`{}_{1/2}{}^{}(m=1/2)`$ sub-state is prepared by optical pumping with $`\sigma ^{}`$ radiation at 397 nm.
3. Optional sideband cooling step (see also section 4): The S$`{}_{1/2}{}^{}(m=1/2)`$ D$`{}_{5/2}{}^{}(m=5/2)`$ transition is excited on one of the red sidebands. Therefore, the laser frequency is detuned red by one motional frequency ($`\omega _{laser}=\omega _{SD}\omega _{trap}`$). The laser power is chosen so that approximately 1 mW laser power is focused to a waist size of 30 $`\mu `$m (spherical trap), or 6 $`\mu `$m (linear trap). The laser at 854 nm is switched on to broaden the D5/2 level at a power level which is set for optimum cooling. Optical pumping to the S$`{}_{1/2}{}^{}(m=+1/2)`$ level is prevented by interspersing short laser pulses of $`\sigma ^{}`$-polarized light at 397 nm. The duration of those pulses is kept at a minimum to prevent unwanted heating.
4. Spectroscopy, or alternatively, state engineering step: We excite the S$`{}_{1/2}{}^{}(m=1/2)`$ D$`{}_{5/2}{}^{}(m=5/2)`$ transition at 729 nm with a single laser pulse, or a number of pulses of well controlled frequency, power, and timing. These parameters are chosen according to the desired state manipulation.
5. Final state analysis: The ion’s fluorescence is collected under excitation with laser light at 397 nm and 866 nm and thus we detect whether a transition to the shelving level D5/2 has been previously induced.

This sequence is repeated typically 100 times to measure the D5/2 state population $`P_D`$ after the state engineering step. We study the dependence of $`P_D`$ on the experimental parameters, such as the detuning $`\delta \omega `$ of the light at 729 nm with respect to the ionic transition, or the length of one of the excitation pulses in step four. The duration of a single sequence is typically 20 ms, so we can synchronize the sequence with the ac line frequency at 50 Hz to reduce ac-magnetic field fluctuations.

A detailed study of the Zeeman and vibrational structure of the S$`{}_{1/2}{}^{}`$ D5/2 transition for a single ion in the linear trap was performed by Nägerl et. al. . For the carrier transition (electronic excitation only, $`nn`$) we measure a 1 kHz linewidth. We attribute this residual broadening to laser and magnetic field fluctuations.

Recently, we excited Ramsey fringes using two consecutive pulses. See Fig. 2 for the Ramsey signal of the S$`{}_{1/2}{}^{}`$ D5/2 carrier transition for a single ion in the spherical Paul trap. Ramsey spectroscopy allows us to investigate the coherence of superposition states $`\frac{1}{\sqrt{2}}\{|S_{1/2}+|D_{5/2}\}`$ which are excited with the first pulse. The fit to the data of Fig. 2, allows us to estimate the purely transversal decay constant to be 2 kHz. We attribute this value to laser phase fluctuations and mechanical vibrations of the trap electrodes.

In the last part of this section we will focus on spectroscopy of ion crystals: Ions in a string are strongly coupled by the Coulomb interaction. Small displacements from their equilibrium positions are described in terms of normal modes of the entire ion crystal vibrating at distinct frequencies . As an example, consider two ions confined in a Paul trap. The first normal mode corresponds to an oscillation of the entire crystal of ions moving back and forth as if they were rigidly joined. This oscillation is referred to as the center-of-mass mode (COM) . The second normal mode, the so-called breathing mode, corresponds to an oscillation where the ions move in opposite directions. The frequency of this mode is at $`\sqrt{3}\omega _{COM}`$.

In the experiment two ions were trapped in the spherical trap, and the excitation spectrum exhibits additional sideband resonances. Apart from the center-of-mass vibration frequencies, which remain at $`\omega _x`$, $`\omega _y`$, $`\omega _z`$, one observes additional sidebands which are identified in Fig. 3. The understanding of the full motional sideband structure presents the basis of further experiments, such as ground state cooling and driving coherent dynamics.

## 4 Cooling to the vibrational ground state

Preparation of the motional ground state is accomplished by a two-stage cooling process. First, the ion is pre-cooled on the S1/2 \- P1/2 dipole transition. In the second stage, a resolved-sideband cooling scheme, similar to the one used in Ref. , is applied on the narrow S1/2 \- D5/2 quadrupole transition: The laser frequency is detuned from the line center by the trap frequency, $`\omega _{laser}=\omega _{SD}\omega _{trap}`$ (first ”red sideband” excitation), thus removing one phonon with each electronic excitation process. The cooling cycle is closed by a spontaneous decay to the ground state which conserves the phonon number in the Lamb-Dicke regime. When the vibrational ground state $`|n=0`$ is reached, the ion decouples from the laser excitation. Due to the weak coupling between light and atom on a quadrupole transition one would expect long cooling times. However, the cooling rate is greatly enhanced by (i) strongly saturating the transition and (ii) shortening the lifetime of the excited state via coupling to a dipole-allowed transition.

After sideband cooling, the ground state occupation is determined by probing the sideband absorption immediately after the cooling pulse. Fig. 4 shows the excitation probability $`P_D`$ for frequencies $`\omega `$ centered around the red and blue $`\omega _z`$ sidebands. For the quantitative determination of the vibrational ground state occupation probability $`p(n=0)`$ after sideband cooling, we compare $`P_D`$ for excitation at $`\delta \omega =\omega `$ and $`\delta \omega =+\omega `$, i.e. on the red and on the blue sideband. For a thermal phonon probability distribution $`p(n)`$, the ratio of excitation on the red and blue sidebands is given by $`P_{red}/P_{blue}=n/(1+n)`$. A Lorentzian fit of the sideband heights yields a ground state occupation of $`p_0`$ = 99.9$`\%`$ for the axial mode when $`\omega _z=(2\pi )\mathrm{\hspace{0.25em}4.51}`$ MHz . By cooling the radial mode with $`\omega _y=(2\pi )\mathrm{\hspace{0.25em}2}`$ MHz, we transfer 98$`\%`$ of the population to the motional ground state. Ground state cooling is also possible at lower trap frequencies, however slightly less efficient. At trap frequencies of $`\omega _z=(2\pi )\mathrm{\hspace{0.25em}2}`$ MHz, $`\omega _y=(2\pi )\mathrm{\hspace{0.25em}0.92}`$ MHz we achieve $`n_z=0.95`$ and $`n_y=0.85`$, respectively. The x-direction is left uncooled because it is nearly perpendicular to the cooling beam. We also succeeded in simultaneously cooling all three vibrational modes by using a second cooling beam and alternating the tuning of the cooling beams between the different red sidebands repeatedly . For two ions in the spherical trap we achieved ground state cooling on the center-of-mass mode in the y and z directions, and on the breathing-mode with more than 96$`\%`$ ground state probability.

In the linear trap, our goal is to implement small scale quantum processing, thus, ions have to be kept apart so that they can be individually addressed. Thus, the laser beam at 729 nm is tightly focused into the linear Paul trap. We have demonstrated individual manipulation of single ions in strings . However, the addressing resolution is restricted by the quality of the optical system to 6 $`\mu `$m. If we demand individual addressing, this resolution restricts our ion - ion distances to greater than 6 $`\mu `$m, which corresponds to axial trap frequencies below 700 kHz (with up to three ions). This low trap frequency aggravates the difficulty of ground state cooling, since the corresponding Doppler cooling limit is close to 20 phonons. Still, we succeed in cooling the axial mode at 700 kHz and reach $`p_0`$= 90$`\%`$ ground state. We also achieved 95$`\%`$ in the ground state for cooling on the ”rocking mode” at $`\omega _{rock}=\sqrt{\omega _{radial}^2\omega _{axial}^2}`$, see Fig. 5. More advanced techniques, such as pre-cooling on the second sideband, prior to the usual first sideband cooling, will probably further increase the ground state occupation.

## 5 Coherent manipulation

Starting from the vibrational ground state, arbitrary quantum states can be created. For a demonstration of coherent state engineering, and in order to investigate decoherence, in step 4 we excite Rabi oscillations with the ion initially prepared in Fock states of its motion. Radiation at 729 nm is applied on the blue sideband transition $`|`$S$`,n_z|`$D$`,n_z+1`$ for a given interaction time $`\tau `$ and the excitation probability $`P_D`$ is measured as a function of $`\tau `$ . The Rabi flopping behaviour allows us to analyse the purity of the initial state and its decoherence . Fig. 5a shows $`P_D(\tau )`$ for the $`|n=0`$ state prepared by sideband cooling. Rabi oscillations at $`\mathrm{\Omega }_{01}=(2\pi )\mathrm{\hspace{0.25em}21}`$ kHz are observed with high contrast indicating that coherence is maintained for times above 1 ms. For the preparation of the Fock state $`|n=1`$, we start from $`|`$S$`,n=0`$, apply a $`\pi `$-pulse on the blue sideband and an optical pumping pulse at 854 nm to transfer the population from $`|`$D$`,n=1`$ to $`|`$S$`,n=1`$. As shown in Fig. 5b for the $`|n=1`$ initial state, we also observe Rabi oscillations with a high contrast, now at $`\mathrm{\Omega }_{12}=(2\pi )\mathrm{\hspace{0.25em}30}`$ kHz. The Rabi frequencies have been theoretically investigated and Blockley et. al. found $`\mathrm{\Omega }_{n,n+1}\sqrt{n+1}`$ in the Lamb-Dicke regime , which is fulfilled in good approximation for our trap. For the ratio of Rabi frequencies $`\mathrm{\Omega }_{12}/\mathrm{\Omega }_{01}`$ in the case of Fock states $`|n=0`$ and $`|n=1`$ we thus expect $`\sqrt{2}`$. The experimental finding agrees with $`\sqrt{2}`$ within 1$`\%`$. The Fourier transform of the flopping signals also yields directly the occupation probabilities for the contributing Fock states $`|n=0,1,2,3\mathrm{}`$  and allows us to calculate the purity of the prepared and manipulated states. For the ’vacuum’ state $`|n=0`$, we obtain $`p_0=0.89(1)`$ with impurites of $`p_1=0.09(1)`$ and $`p_{n2}0.02(1)`$. For the Fock state $`|n=1`$ the populations are $`p_0=0.03(1)`$, $`p_1=0.87(1)`$, $`p_2=0.08(2)`$, and $`p_{n3}0.02(1)`$. The measured transfer fidelity is about 0.9. Note that the Rabi flopping data here were taken with less efficient cooling (at a lower trap frequency), and the number state occupation from the Fourier analysis is consistent with the occupation which we determined by sideband measurements.

## 6 Speed limits of gate operations

As a figure of merit for an ion trap quantum processor, we calculate how many gate operations are feasible within the time scale given by decoherence. We have investigated the speed limit of quantum gate operations. The proposed methods  employ the use of laser pulses to entangle the electronic and the vibrational degrees of freedom of trapped ions. A theoretical investigation shows that the proposed methods are limited mainly by the recoil frequency of the relevant electronic transition . We have experimentally studied the basic building block of a gate operation, that is a $`\pi `$-pulse excitation on the first blue sideband. If the blue sideband is driven significantly faster than the inverse of the recoil frequency, carrier transitions are excited off-resonance and the contrast of the observed Rabi oscillations decreases. In the case of the quadrupole S$`{}_{1/2}{}^{}`$ D5/2 transition we find a time of about 16 $`\mu `$s for this limit, if we demand a gate fidelity of 99 $`\%`$. For further details on fast gates, we refer to A. Steane et. al. . We thus estimate that 30 - 50 quantum gate operations will be feasible with the current setup.

## 7 Conclusion

We have engineered the quantum states of motion $`|n=0,1`$ of a single trapped ion that are relevant for quantum computation, using laser excitation on a forbidden optical transition. We have observed more than 30 periods of Rabi oscillations on the motional sidebands of this transition, thus showing that decoherence is negligible on the time scale of a few oscillations, i.e. the time required for a quantum gate operation. We attribute the observed 1 ms decoherence time to laser and magnetic field fluctuations. Heating of the motional degrees of freedom has also been observed and was measured directly to happen at least 1 order of magnitude more slowly. This confirms that in the comparatively large traps which we use, heating seems not to be a limiting process.

Furthermore, we have been able to cool two ions to the ground state, and also to reach the ground state for single and two ions in the linear trap, under conditions which allow individual addressing of ions in a string. Multiple coherent gate operations with trapped ions seem well within experimental reach.

This work is supported by the Austrian ”Fonds zur Förderung der wissenschaftlichen Forschung” within the project SFB15, and in parts by the European Commission within the TMR networks ”Quantum Information” (ERB-FMRX-CT96-0087) and ”Quantum Structures” (ERB-FMRX-CT96-0077), and the ”Institut für Quanteninformation GmbH”.