text:
# X-ray fluorescence holography: going beyond the diffraction limit

## I Introduction

X-ray fluorescence holography (XFH) is a promising method for directly imaging local atomic structure in an element-specific way that has been developed experimentally over the past six years AIP:1986 ; Tegze:1996 ; Gog:1996 . In the first mode of measurement, an atom inside an oriented system (e.g. a single crystal with low to medium mosaicity) emits a fluorescent x-ray wave; the interference between the outgoing unperturbed reference wave component and the scattered object wave component produces a hologram AIP:1986 ; Tegze:1996 . One can term this “inside source” holography. A time-reversed mode also exists, in which the exciting x-ray wave is the reference, and scattered x-rays impinging on a fluorescent emitter are the object waves Gog:1996 . The emitter is the detector of the hologram, and thus this can be called “inner detector” holography. The holographic reconstruction in either case is achieved by propagating back the wavefield from the far field.

Beyond several recent achievements of these methods Fitzgerald:2000 , e.g. in imaging the environment around a dilute semiconductor dopant Hayashi:1998 , low-Z atoms in the presence of intermediate-Z atoms Tegze:2000 , and the local environment in a quasicrystal Marchesini:2000 , another step toward the practical application of holography would be a quantitative link between the electronic charge distribution and the holographic image. It has already been demonstrated that images with accuracies at the diffraction limit can be obtained Tegze:1999 . But even in this case, the reconstructed wavefield can be viewed as a distorted image of the scatterer density distribution; i.e. it is the convolution of the charge density distribution with a point spread function due to the measurement and inversion process.

When an object is imaged through an optical system, information is lost whenever the imaging system cannot pass all the spatial frequencies contained in the scene. Further loss is produced by aberrations in the optical system. The question that we here address is how much information can be recovered if we know the point spread function (PSF) of the effective optical system associated with x-ray fluorescence holography. We will consider here several methods for the iterative deconvolution of data so as to improve image quality.

Several procedures for iterative image improvement have been developed previously in the fields of astronomy and microscopy Stark:1987 , and we begin by introducing them briefly. Van Cittert developed a very early method in 1920 Van:1970 ; this uses a self-consistent iterative solution to the problem. It can be shown that this method is close to the gradient search of the maximum likelihood (ML). Two other approaches are based on the probability theory. One developed by Lucy Lucy:1974  and Richardson Richardson:1972  is based on Bayes theorem for conditional probability. Another one, the maximum entropy method (MEM) seeks the ‘most probable’ solution in an under-determined system of equations Press:1992 . In this article, we discuss some of these methods that seem most appropriate for applying to XFH image improvement, quantitatively assess several of them via model theoretical calculations, and finally arrive at the first deconvolved atomic image at resolution going beyond the diffraction limit.

## II Basic imaging concepts in XFH

We first consider some basic imaging concepts for XFH, in particular for an arbitrary atomic charge density distribution $`\rho \left(𝐫\right)`$. For simplicity, we will also throughout this discussion consider only “inner source” holography, although the conclusions here can easily be generalized to include “inner detector” holography. Neglecting absorption or extinction effects (a reasonable approximation if we consider only near-neighbor imaging), the hologram $`\chi `$ in the far field $`𝐤`$ is given by Faigel:1999 :

$`\chi \left(𝐤\right)`$ $`=`$ $`{\displaystyle \eta (𝐤,𝐫)\rho \left(𝐫\right)d^3𝐫},\text{ with}`$
$`\eta (𝐤,𝐫)`$ $`=`$ $`2\mathrm{R}\mathrm{e}\left[\frac{r_e}{r}e^{i(kr𝐤𝐫)}f(𝐤,𝐫)\right];`$ (1)
here $`r_e`$ is the classical electron radius, $`r_ef(𝐤,𝐫)`$ is the scattering factor per electron and the phase factor is given by the path length difference between the wave emitted from the origin and the wave scattered by the electron located in $`𝐫`$. The typical reconstruction algorithm to calculate the wavefield $`U`$ in the real space $`𝐫`$ is Barton:1988 :

$`U\left(𝐫^{}\right)`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Omega }}}{\displaystyle _\sigma }\eta ^{\{1\}}(𝐫^{},𝐤)\chi \left(𝐤\right)d^3𝐤,\text{ with}`$
$`\eta ^{\{1\}}(𝐫^{},𝐤)`$ $`=`$ $`\frac{r^{}}{r_e}\mathrm{Re}\left[e^{i\left(kr^{}𝐤𝐫^{}\right)}\right];`$ (2)
here $`\sigma `$ is the measured region in k-space, $`\mathrm{\Omega }`$ the volume of $`\sigma `$, and we have included the scaling factor $`\frac{r^{}}{r_e}`$ to obtain an estimate of the charge density. Sometimes an additional factor $`f^{\{1\}}`$ is included to compensate the ‘aberration’ or distortions produced in the hologram $`\chi `$ by $`f`$Tonner:1991 , but we will not explicitly include this here. We will for simplicity assume that the scattering factor of an electron is isotropic, neglecting the angular dependence of the Thompson scattering and near field effects, The holographic reconstruction U can thus be considered to be a distorted image of the charge density distribution $`\rho `$:

$`U\left(𝐫^{}\right)`$ $`=`$ $`{\displaystyle \mu (𝐫^{},𝐫)\rho \left(𝐫\right)d^3𝐫},\text{ with}`$
$`\mu (𝐫^{},𝐫)`$ $`=`$ $`{\displaystyle _\sigma }d^3𝐤\eta ^{\{1\}}(𝐫^{},𝐤)\eta (𝐤,𝐫);`$ (3)
with $`\mu \left(𝐫^{}\text{,}𝐫\right)`$ now being the experimental point spread function. An additional factor $`\left(1+\mathrm{sign}U\left(𝐫^{}\text{,}𝐫\right)\right)/2`$ can be used to enforce positivity. In summary we have three spaces $`𝐫`$, $`𝐤`$ and $`𝐫^{}`$ for which we have three functions $`\rho \left(𝐫\right)`$, $`\chi \left(𝐤\right)`$ and $`U(𝐫^{})`$ and the propagators from one space to the other are $`\eta (𝐤,𝐫)`$, $`\eta ^{\{1\}}(𝐤,𝐫^{})`$ and $`\mu \left(𝐫^{}\text{,}𝐫\right)`$.

We now turn to specific methods for iterative deconvolution of images. For notational brevity, we in the following sections use the symbol $``$ between two functions to denote integration over the internal variable, e.g. $`f(x,y)g(y)=𝑑yf(x,y)g(y)`$.

## III Iterative image deconvolution

The general problem of the deconvolution is to obtain the unknown object $`\rho \left(𝐫\right)`$ from a knowledge of $`U(𝐫^{})`$ and $`\mu (𝐫^{},𝐫)`$. The first method developed by Van Cittert in 1920 Van:1970  is the most intuitive one and it is based on a self-consistent solution. Every step is the difference between the reconstructed experimental hologram $`U\left(𝐫\right)`$ and the reconstructed simulated hologram $`U^{\{n\}}\left(𝐫\right)=\mu (𝐫,𝐫^{})\rho ^{\{n\}}\left(𝐫\right)`$ generated by a new $`n^{th}`$-step charge density distribution. That is,

$`\rho ^{\left\{n+1\right\}}\left(𝐫\right)`$ $`=`$ $`\rho ^{\left\{n\right\}}\left(𝐫\right)+\mathrm{\Delta }\rho ^{\left\{n\right\}}\left(𝐫\right),`$
$`\mathrm{\Delta }\rho ^{\left\{n\right\}}\left(𝐫\right)`$ $`=`$ $`\left[U\left(𝐫\right)U^{\left\{n\right\}}\left(𝐫\right)\right],`$
$`\rho ^{\left\{1\right\}}\left(𝐫\right)`$ $`=`$ $`U\left(𝐫\right).`$ (4)
Thus, $`\mu (𝐫^{},𝐫)`$ is in this approximation assumed to be a delta function $`\delta (𝐫^{},𝐫)`$ over each step in the iteration. We normally enforce positivity after each step by setting to 0 the negative values of the charge density.

A similar problem is the phase retrieval of a diffraction pattern of known intensity from an object with known support $`s\left(𝐫\right)`$. The diffraction pattern provides the amplitude of Fourier transform of the charge density $`\left|\stackrel{~}{\rho }\right|`$. If the object $`\rho \left(𝐫\right)`$ is constrained within a given support (e.g. with the support equal to 1 over the object and 0 outside of it), the product of the object $`\rho \left(𝐫\right)`$ with the support function $`s\left(𝐫\right)`$ must be still equal to the object. Therefore the complex amplitude of the diffracted wave-field is not modified by a convolution with the Fourier transform of the support function $`\stackrel{~}{s}`$:

$$\rho =s\rho \stackrel{FT}{}\stackrel{~}{\rho }=\stackrel{~}{s}\stackrel{~}{\rho }.$$
(5)
If $`\stackrel{~}{s}`$ is sufficiently large (‘over-sampling’ condition Miao:1999 ) a certain number of pixels are bounded by this self-consistent equation. The main difference with iterative deconvolution methods is that in this case, the bigger $`\stackrel{~}{s}`$ is, the larger is the number of pixels bounded by (5), and the easier it is to find the solution. If the over-sampling condition is satisfied, iterative approaches such as the Gerchberg-Saxton Gerchberg:1972  and Fienup Fienup:1982  algorithms can be used to solve the phase problem. In the Fourier domain, the Gerchberg-Saxton algorithm can be expressed as:

$$\stackrel{~}{\rho }^{\left\{n+1\right\}}=\left|\stackrel{~}{\rho }_{\text{meas}}\right|e^{i\text{phase}\left\{\stackrel{~}{\rho }^{\left\{n\right\}}\stackrel{~}{s}\right\}},$$

where the $`n^{th}`$ and $`(n+1)^{th}`$ estimates for density are given, and $`\left|\stackrel{~}{\rho }_{\text{meas}}\right|`$ is the measured amplitude of the Fourier transform of the object.

Another iterative approach is given by the Lucy-Richardson method Lucy:1974 ; Richardson:1972 . This method is based on a self-consistent iterative solution of the Bayes theorem for conditional probabilities, considering $`\mu (𝐫^{},𝐫)`$ as a probability that $`𝐫^{}`$ will fall in the interval $`(𝐫^{},𝐫^{}+d𝐫^{})`$ when it is known that $`𝐫=𝐫^{}`$. A self-consistent iterative approach can in this case be expressed as Lucy:1974 :

$`\rho ^{\left\{n+1\right\}}\left(𝐫\right)=`$ $`\rho ^{\left\{n\right\}}\left(𝐫\right)`$ (6)
$`\times \left[\mu ^T(𝐫,𝐫^{}){\displaystyle \frac{U\left(𝐫^{}\right)}{\mu (𝐫^{},𝐫^{\prime \prime })\rho ^{\left\{n\right\}}\left(𝐫^{\prime \prime }\right)}}\right]`$
with $`\mu ^T(𝐫,𝐫^{})`$ being the transpose of $`\mu `$.

## IV Maximum likelihood and entropy

In another set of methods, we want to minimize the difference between the measured hologram in the far field $`\chi \left(𝐤\right)`$ and the hologram generated by a given charge density distribution $`\rho ^{\left\{n\right\}}\left(𝐫\right)`$:

$$𝐐=_\sigma d^3𝐤\left|\chi \left(𝐤\right)\eta (𝐤,𝐫)\rho ^{\left\{n\right\}}\left(𝐫\right)\right|^2.$$
(7)
where the unknowns are $`\rho ^{\left\{n\right\}}\left(𝐫_i\right),i=1\mathrm{}N`$ over some grid spanning the object. The gradient can be computed as:

$$𝐐=2\eta ^{}(𝐫,𝐤)\left[\eta (𝐤,𝐫)\rho ^{\left\{n\right\}}\left(𝐫\right)\chi \left(𝐫\right)\right].$$
(8)
with $`\eta ^{}`$ being the transpose conjugate of $`\eta `$. The Hessian of $`𝐐`$ is in this case nothing else but the point spread function and is given by:

$`𝐇`$ $`=`$ $`{\displaystyle \frac{^2𝐐}{\rho ^{\left\{n\right\}}\left(𝐫_1\right)\rho ^{\left\{n\right\}}\left(𝐫_2\right)}}`$ (9)
$`=`$ $`2\eta ^{}(𝐫_1,𝐤)\eta (𝐤,𝐫_2).`$
The single pixel approximation Wilkins:1983 ; Cornwell:1985  consists in assuming that $`𝐇`$ is diagonal, so that the steps will be:

$$\mathrm{\Delta }\rho ^{\left\{n\right\}}\left(𝐫\right)=\frac{1}{𝐇_{𝐫,𝐫}}𝐐_𝐫.$$
(10)
Under this approximation, the gradient search is equivalent to the Van Cittert method, since $`H_{𝐫,𝐫}\left(\frac{r_e}{r}\right)^2`$ (see proof in Appendix A) and $`\eta ^{\{1\}}=\frac{r^2}{r_e^2}\eta ^{}`$. However, for single energy holograms, at the inversion-symmetric position to the real atoms, there will appear a twin or ghost image. Therefore one should consider at least the symmetric term $`𝐇_{𝐫,𝐫}`$:

$`\overline{𝐇}_{𝐫,𝐫^{}}`$ $`=`$ $`𝐇_{𝐫,𝐫^{}}\left(\delta _{𝐫^{},𝐫}+\delta _{𝐫^{},𝐫}\right)`$ (11)
$`\mathrm{\Delta }\rho ^{\left\{n\right\}}\left(𝐫\right)`$ $`=`$ $`\left(\begin{array}{cccccccccccccccccccc}\overline{𝐇}_{𝐫,𝐫}& \overline{𝐇}_{𝐫,𝐫}& & & & & & & & & & & & & & & & & & \\ \overline{𝐇}_{𝐫,𝐫}& \overline{𝐇}_{𝐫,𝐫}& & & & & & & & & & & & & & & & & & \end{array}\right)^1\left(\begin{array}{cccccccccccccccccccc}𝐐_𝐫& & & & & & & & & & & & & & & & & & & \\ 𝐐_𝐫& & & & & & & & & & & & & & & & & & & \end{array}\right).`$ (16)
An additional consideration here however, is that at every energy, for a given position, and in case of a point scattering, the twin image cancels the real image, so that a divergence arises in Eq. (16). This can be avoided by introducing an anisotropy in the scattering factor between the forward and backward directions. It is known Szoke:1993  that, if we use a small Gaussian charge distribution for every point $`𝐫_i`$ instead of point scatterers, the scattering becomes stronger in the forward direction, and iterative algorithms converge more quickly. Conjugate gradient methods do not require the calculation of the Hessian, but partial knowledge of it is helpful. Extra information can also be used to speed up the process, and for example, we know that the scattering charge distribution is real and positive. This can be enforced by taking the real part at every step and by setting to $`\rho _n=0`$ when $`\rho _n\left(𝐫\right)<0`$. A further piece of information is the total number of scattering electrons $`𝐅_0`$, and we can use Lagrange multipliers to perform a constrained best fit via:

$`𝐐+\beta \left(𝐅𝐅_0\right),\text{ with}`$
$`𝐅={\displaystyle \rho \left(𝐫_i\right)}.`$ (17)
Here the minimization is performed adjusting the multiplier at every step:

$`\mathrm{\Delta }\rho _𝐫`$ $`=`$ $`\overline{𝐇}_{𝐫,𝐫^{}}^1\left(𝐐_𝐫^{}+\beta 𝐅_𝐫^{}\right),`$
$`\lambda `$ $`=`$ $`{\displaystyle \frac{\left[𝐅_0𝐅+\left(𝐅\overline{𝐇}^1𝐐\right)\right]}{𝐅\overline{𝐇}^1𝐅}}`$ (18)
The Lagrange multiplier method keeps the number of electron fixed, however after we apply each step, we enforce positivity in the charge density by setting to 0 the negative values of the charge density, therefore the actual number of electrons will tend to grow at every iteration. In the actual implementation of the algorithm, we artificially reduce the step size by a factor of 1/5 to avoid cutting too much in the negative density regions.

The deconvolution of spectra with an incomplete Fourier image is also often treated with the maximum entropy method Press:1992 . When the number of unknowns $`\rho \left(𝐫_i\right)`$ is superior to the number of equations, it is necessary to include further a-posteriori information. As more than one solution can satisfy the set of equations, we seek the most probable solution, i.e. the one with the maximum entropy. In practice we want to minimize the functional:

$`𝐒`$ $`+\lambda ^{}𝐐,`$ (19)
$`𝐒`$ $`={\displaystyle _i}\rho \left(𝐫_i\right)\mathrm{ln}\rho \left(𝐫_i\right)`$
where $`\lambda ^{}`$ is a Lagrange multiplier. The procedure is similar to the previous one, with two main differences: (i) MEM enforces positivity, so that a pure gradient minimization could lead to unphysical values that need to be chopped, (ii) $`𝐒`$ is highly nonlinear, and simple gradient search is not efficient. An interesting property of $`𝐒`$ is that its Hessian is diagonal; therefore the single pixel approximation can be used Cornwell:1985 .

## V Numerical simulations

We have initially tried all the above mentioned deconvolution methods on a test system consisting of a charge distribution described by a 3D matrix of $`50^3`$ points with 0.2 Å spacing. This matrix was set to 0 except for two Gaussian charge distributions centered in the x-z plane representing two atoms separated by 3.5 Å as shown in Fig. (1a). Figures (1b-d) now show reconstructed images via several different methods obtained from a hologram (Fig. 2a) defined at an energy of 17 keV with the points distributed over a solid angle with azimuth varying from $`0^0`$ to $`360^0`$ and polar angle from $`0^0`$ (parallel to the z-axis) to $`80^0`$. The reconstruction in Fig. (1f) has been obtained using both holograms at 17 keV and 18 keV (Fig. 2b). The standard direct Barton method in Fig. (1b) shows a weak twin image, and some distortion of the imaged charge distributions, with the difference between Fig. (1a) and Fig. (1b) essentially being the PSF. The number of electrons estimated from the Barton method is one order of magnitude bigger than the original one. This could be ascribed to the fact that many of these electrons are not contributing to the hologram, as they cancel each other. In trials with the Van Cittert and ML (Fig. 1c) methods using the double pixel approximation, we have found that the convergence is very good: less then 10 iterations are necessary for convergence. This is due to the fact that the holographic reconstruction is a relatively good image of the charge density. However the number of electrons used to fit the data was one order of magnitude bigger than that of the original image. Another relevant feature of this method is the fact that the image becomes more structured. This is due to the divergence in the inverse of the Hessian using single energy holograms when $`𝐇_{𝐫,𝐫}=𝐇_{𝐫,𝐫}`$.

Tests using the Lucy-Richardson method have shown bad convergence; further studies of this method need to be done, but one of the reasons could be some divergence caused by the denominator in equation (6).

At this point we turned to the maximum entropy method. As we had already developed the ML algorithm with single and double pixel approximations, we could easily apply the algorithm developed by Cornwell and Evans Cornwell:1985  so as to employ the maximum entropy criterion. This algorithm tries to maximize the entropy $`𝐒`$ while minimizing the likelihood $`𝐐`$ with a fixed number of electrons $`𝐅`$ using Lagrange multipliers. Several tests have shown that this algorithm (at least the one implemented by us) was not able to minimize $`𝐐`$, $`𝐒`$ and $`𝐅`$ at the same time; the result is shown in Fig. (1d). However, we did find that the maximum likelihood and the minimum number of electrons produces a good quality reconstruction, with somewhat improved resolution, and a reasonably good $`𝐅`$Marchesini:2001 ; these results are shown in Fig. (1e) and they can be compared with those in (1b). The number of electrons is only 20% bigger than the original one.

The single energy iterative reconstruction shows how the divergence in the Hessian due to the twin image problem produces several artifacts, which are removed when using slightly different energies (Fig. 1f). The reason why the number of electrons does not correspond to the original $`𝐅_0`$ even when applying the constraint is that after the steps are applied, we impose positivity by setting to 0 the negative values of the density.

Finally, we applied the constrained ML method to the same case, using a double-energy data set. This result is shown in Fig. (1f) and represents the most significant enhancement in image quality relative to the standard method.

To explore the constrained ML algorithm further, we carried out simulations based on a much larger and more realistic cluster of atoms representing CdTe; the charge density distribution in the neighborhood of the origin is shown in Fig. (3a,b). We tried to keep a minimum amount of information to show the robustness of the method: We simulated the hologram produced by a CdTe cluster of 1000 atoms at only two energies of 8 and 9 keV. At these two energies the holographic reconstruction using the Barton method shows a poor image quality as shown in Fig. (3c,d). We also introduced a randomly distributed noise level of $`10^3`$ on every data point (corresponding to 106 counts in each point), with the points distributed over a solid angle with azimuth varying from 0 to 360 degrees and polar angle from 0 (perpendicular to the  surface normal direction) to 80 degrees, with a 1-degree step size in each direction. Since we are interested in imaging the first neighbors, we finally carried out all imaging with a smoothed hologram based on 2 degree smoothing in every direction. Although the cluster size is several nanometers, we will only try to reconstruct up to 8 Å setting to 0 anything outside this range.

The holograms generated by this charge after smoothing are shown in Fig. (4a) at 8 keV and Fig. (4b) at 9 keV. We again choose as a first reference step and set of images the standard holographic reconstruction in Fig. (3a,b). Comparison between the standard reconstruction and the deconvolved images are shown in Fig. (3c,f).

The number of electrons within the reconstructed region is $`𝐅_0=1.8\times \mathrm{\hspace{0.33em}10}^4`$ in the original picture, $`𝐅=7\times 10^4`$ in the standard reconstruction, and $`𝐅=2.9\times \mathrm{\hspace{0.33em}10}^4`$ in the deconvolved image. Thus, the deconvolved image is a factor of 1.5 larger than the original image in number of electrons, but a factor of 2 better than the standard reconstruction. This result is worse than what we observed for the two Gaussian charges, and can be ascribed to the difficulty to fit point like structures (atoms) at the diffraction limit scale instead of smooth distributions. The peak heights of the charge density in the three images varies somewhat differently, the maximum number of electrons in one voxel in the original image is 50, 0.38 in the standard reconstruction, and 14 in the deconvolved image, as shown more quantitatively in the colorbars of Fig. (3). The reason for this is that in the original image, every electron of one atom is concentrated in a single voxel, while the standard reconstruction has many electrons spread over many more voxels. The deconvolved image has a total of fewer electrons than the standard image, but much higher peaks, due to the enhanced resolution relative to the standard reconstruction. Comparing image quality between the standard (Fig. 3c,d) and deconvolved (Fig. 3e,f) approaches, we see a definite enhancement of resolution in the deconvolved images, even though some artifacts remain that are of comparable strength to the actual charge-density images. The holograms generated by the cluster in Fig. (4a,b) are well described in the low frequency range by the fitted holograms in Fig. (4c) (8 keV) and Fig. (4d). Since in the iterations we only fitted the charge density within the first 8 Å setting the distribution to 0 outside this range, only the low frequency components where fitted, as the residual difference shows in Fig. (4e-f).

In Fig. (5), we present some results dealing with the rapidity and type of convergence observed. The maximum likelihood is obtained within the first 10-20 iterations (Fig. 5-top). The total number of electrons is slowly reduced in the following $`500`$ iterations (Fig. 5-center). The Lagrange multiplier $`\lambda `$ follows the behavior of the constraint: as $`𝐅`$ changes from the desired level $`𝐅_0`$, $`\lambda `$ varies in order to bring $`𝐅`$ back to the desired level.

## VI Conclusions and future possibilities

In conclusion, we have demonstrated that iterative deconvolution can be applied to XFH imaging, thus improving the quality of the images beyond the diffraction limit, and that several levels of approximation can be used. From our test cases, the constrained maximum likelihood method appears to be the most promising. In these first theoretical simulations related to image improvement, we have approximated the scattering factor as isotropic, thus neglecting the angular dependence of Thomson scattering and near-field effects. As a computationally relevant issue, we have also approximated the Hessian to be a sparse matrix where only the two diagonals are non-zero. All of these approximations could be relaxed in future applications, including terms close to the diagonals, with of course additional complexity then being incurred in the deconvolution algorithm and the length of time it takes. For our calculations, the storage and calculation of the full Hessian of $`80^6`$ points required considerable memory, but computers capable of such calculations should be available as desktops in a few years. The gradient was also calculated neglecting the angular dependence of the scattering factor. Using these approximations, every iteration involved the calculation of a hologram from a matrix of $`80^3`$ voxels, the reconstruction and required about 5 minutes on a Pentium II computer, for a total of one night of calculation. Our image deconvolution has been obtained by enforcing a real and positive charge density distribution and limited number of electrons. However, looking ahead to further image improvement, once the image is deconvolved, one could make use of ‘atomicity’ in the next step by placing atoms in the location of the peaks, and then performing an optimization using the full angular dependence of the scattering factor, while simultaneously adjusting the atomic number and position of these atoms. This study thus suggests several fruitful directions for further exploration and exploitation of image deconvolution in x-ray fluorescence holography that should considerably enhance its power for structural studies.

###### Acknowledgements.

We would like to acknowledge A. Szöke for fruitful discussions. This work was supported by the Laboratory Directed Research and Development Program of Lawrence Berkeley National Laboratory under the Department of Energy Contract No. DE-AC03-76SF00098. NATO Collaborative Linkage Grant for travel support.

## Appendix A Appendix–Calculation of the Hessian and the gradient

The Hessian does not need to be known in great detail. It can be pre-calculated before starting the iterative algorithm neglecting the angular dependence of Thompson scattering factor $`\frac{f_T}{r_e}=\frac{1+\mathrm{cos}^2\mathrm{\Theta }_{𝐤,𝐫}}{2}`$ and using the following approximation $`\underset{\sigma \left(k\right)}{}\mathrm{cos}^2\left(kr𝐤𝐫\right)k^2d^2\widehat{𝐤}\frac{1}{2}\mathrm{\Omega }_k`$ and $`\underset{\sigma \left(k\right)}{}\mathrm{sin}^2\left(𝐤𝐫\right)k^2d^2\widehat{𝐤}\frac{1}{2}\mathrm{\Omega }_k`$, with $`\mathrm{\Omega }_k=_{\sigma \left(k\right)}k^2d^2\widehat{𝐤}`$, $`\widehat{𝐤}`$ is a unit vector, and $`\sigma (k)`$ is the region where the hologram is measured at a given energy. This yields:

$`𝐇_{𝐫,𝐫}`$ $`=`$ $`{\displaystyle _k}\mathrm{\Delta }k{\displaystyle _{\sigma (k)}}\frac{r^2}{r_e^2}\mathrm{cos}^2\left(kr𝐤𝐫\right)k^2d^2\widehat{𝐤},`$
$``$ $`\frac{r^2}{r_e^2}\frac{1}{2}_k\mathrm{\Omega }_k\mathrm{\Delta }k;`$
$`𝐇_{𝐫,𝐫}`$ $`=`$ $`{\displaystyle _k}\mathrm{\Delta }k{\displaystyle _{\sigma (k)}}\frac{r}{r_e}\mathrm{cos}\left(kr𝐤𝐫\right)`$ (20)
$`\times \mathrm{cos}\left(kr+𝐤𝐫\right)k^2d^2\widehat{k}`$
$`=`$ $`\frac{r^2}{r_e^2}_k\mathrm{\Delta }k_{\sigma (k)}\left[\mathrm{cos}^2\left(kr\right)\mathrm{sin}^2\left(𝐤𝐫\right)\right]k^2d^2\widehat{𝐤},`$
$``$ $`\frac{r^2}{r_e^2}_k\left[\mathrm{cos}^2\left(kr\right)\frac{1}{2}\right]\mathrm{\Omega }_k\mathrm{\Delta }k`$
$``$ $`\frac{1}{2}\frac{r^2}{r_e^2}_k\mathrm{cos}\left(2kr\right)\mathrm{\Omega }_k\mathrm{\Delta }k;`$
The gradient is calculated by simulating the hologram $`\chi ^{\left\{n\right\}}\left(𝐤\right)`$ generated by a charge density distribution $`\rho ^{\left\{n\right\}}\left(𝐫\right)`$, and propagating back the hologram to the real space. Using $`\overline{U}\left(𝐫\right)=\left(\frac{r}{r_e}\right)^2U\left(𝐫\right)`$, eq. (8) can be written as:

$`𝐐`$ $`=`$ $`2\left[\overline{U}^{\left\{n\right\}}\left(𝐫\right)\overline{U}\left(𝐫\right)\right];`$
$`\overline{U}^{\left\{n\right\}}\left(𝐫\right)`$ $`=`$ $`\eta ^{}(𝐫,𝐤)\chi ^{\left\{n\right\}}\left(𝐤\right)`$
$``$ $`{\displaystyle \frac{2r}{\mathrm{\Omega }r_e}}\mathrm{Re}{\displaystyle _\sigma }e^{i\left(kr𝐤𝐫\right)}\chi ^{\left\{n\right\}}\left(𝐤\right)d^3𝐤;`$
$`\chi ^{\left\{n\right\}}\left(𝐤\right)`$ $`=`$ $`\eta (𝐤,𝐫)\rho ^{\left\{n\right\}}\left(𝐫\right)`$ (21)
$``$ $`2\mathrm{R}\mathrm{e}{\displaystyle \frac{r}{r_e}e^{i\left(kr𝐤𝐫\right)}\rho ^{\left\{n\right\}}\left(𝐫\right)d^3𝐫}.`$
These calculations can be very long as we need to sum over $`𝐫(x,y,z)`$ for every $`𝐤(k,\vartheta ,\phi )`$. Supposing that the number of points in $`x,y,z,\vartheta ,\phi `$ is $`N100`$ for every variable, we need to calculate the product in the integral of eq. (21) for every possible value, i.e. $`N^5`$ products and sums. The computation of $`\chi ^{\left\{n\right\}}\left(𝐤\right)`$ can be performed faster if we separate the sum over the variables $`x,y,z`$.

We begin by calculating $`\xi ^{\{n\}}\left(𝐫\right)=\frac{r}{r_e}e^{ikr}\rho ^{\left\{n\right\}}\left(𝐫\right)\mathrm{\Delta }r^3`$; we can express Eq. (21) as:

$$\chi _n\left(𝐤\right)2\mathrm{R}\mathrm{e}\underset{𝐫}{}e^{i𝐤𝐫}\xi ^{\{n\}}\left(𝐫\right)$$
(22)
which is a simple Fourier transform. However the holograms are usually measured in spherical coordinates, and to obtain the values of $`\chi _n\left(𝐤\right)`$ in the measured positions, we either need to perform an interpolation, or a convolution. Another approach is to use directly the spherical coordinates in the kernel:

$`e^{i𝐤𝐫}`$ $`=`$ $`e^{ikx\mathrm{sin}\vartheta \mathrm{cos}\phi }e^{iky\mathrm{sin}\vartheta \mathrm{sin}\phi }e^{ikz\mathrm{cos}\vartheta }`$ (23)
$`=`$ $`A_{k,x,\vartheta ,\phi }B_{k,y,\vartheta ,\phi }C_{k,z,\vartheta }`$
and Eq. (22) becomes:

$`\chi ^{\{n\}}(k,\vartheta ,\phi )`$ $``$ $`2\mathrm{R}\mathrm{e}{\displaystyle _x}A_{k,x,\vartheta ,\phi }{\displaystyle _y}B_{k,y,\vartheta ,\phi }`$ (24)
$`\times {\displaystyle _z}C_{k,z,\vartheta }\xi ^{\{n\}}(x,y,z).`$
Each sum requires $`N^4`$ calculation, reducing the time to perform the calculation by $``$N/3. The same trick can be applied for the holographic reconstruction.