3d Simulation of the Image Formation in Soft X- Ray Microscopes References and Links

In water-window soft x-ray microscopy the studied object is typically larger than the depth of focus and the sample illumination is often partially coherent. This blurs out-of-focus features and may introduce considerable fringing. Understanding the influence of these phenomena on the image formation is therefore important when interpreting experimental data. Here we present a wave-propagation model operating in 3D for simulating the image formation of thick objects in partially coherent soft x-ray microscopes. The model is compared with present simulation methods as well as with experiments. The results show that our model predicts the image formation of transmission soft x-ray microscopes more accurately than previous models.ray tomography of phenotypic switching and the cellular response to antifungal peptoids in Candida albicans, " Proc. Characterization of the resolving power and contrast transfer function of a transmission X-ray microscope with partially coherent illumination, " Opt. Validity criterion for the Born approximation convergence in microscopy imaging, " J. Numerical model for tomographic image formation in transmission x-ray microscopy, " Opt.


Introduction
Water-window soft x-ray microscopy allows imaging of intact cells in their near-native "hydrated" state with nanoscale resolution [1].However, the formation of the 2D images as well as their use in 3D tomographic reconstructions is not uncomplicated, primarily for two reasons: the depth of focus in these lens-based microscopes is typically smaller than the sample and the sample illumination is often partially coherent.Accurate description of the image formation process shows promise to alleviate this problem but previous computationally feasible models for thick objects are limited to 2D computations or incoherent imaging.Here we present an improved model for partially coherent microscopes based on 3D wave propagation and show that this model provides a superior description of the image formation.
Soft x-ray microscopes (XRM:s) operating in the water window exploit the natural contrast between carbon and oxygen in the E = 284-543 eV energy region (λ = 4.37-2.29 nm) and rely on zone-plate optics for the high-resolution imaging.XRM:s are operated at several synchrotron sources [2,3] and laboratory-source-based microscopes are emerging [4].The primary advantage of XRM:s is their unique ability to image the intracellular structure of whole intact cells, in 2D as well as in 3D via tomographic reconstruction.The last few years synchrotron-based microscopy has demonstrated significant progress and is now delivering results of high biological relevance [5,6].The challenge in XRM imaging is that the depth of focus (DOF) depends on the wavelength (λ) and numerical aperture of the zone plate objective (NA zp ) as λ/NA zp 2 [7], which does typically not extend through the whole specimen so that only part of the sample volume is imaged sharply.Furthermore, the illuminating condenser systems often have a smaller NA (NA ill ) than the zone-plate, resulting in a coherence parameter [8], m = NA ill /NA zp smaller than 1, i.e. partially coherent illumination.Such partially coherent illumination is often advantageous since it provides higher contrast for higher-spatial-frequency details [6].However, the combination of partially coherent illumination with part of the sample volume being outside the DOF can introduce non-trivial defocus-effects to the image, which, in turn, may cause non-obvious artifacts in the 3D reconstruction if not taken into proper consideration.The magnitude of these effects is sample dependent, where images of small and/or weakly scattering samples are less influenced while large and/or strongly scattering samples may produce significant imaging artifacts.With improved optics having higher numerical apertures the problem worsens, since the DOF is inversely proportional to the square of NA zp .Typically, state-of-the-art nanofabrication can facilitate zone plates that allow imaging down to 10 nm for thin objects but their DOF is short [9].
The complicated image formation process calls for a model that properly describes the resulting image.The purpose of the model is to make better use of resolution and contrast in 2D and 3D imaging by an improved understanding of the effects of defocus and the illumination settings.Transmission soft x-ray microscopes have been modelled previously, both for incoherent and partially coherent illumination.For partially coherent illumination Schneider modelled the image formation of thin objects, providing numerical solutions to 1D objects [10].Von Hofsten et al extended this method to 2D numerical calculations [11].For 3D samples, Knöchel applied the methodology of Streibl [12] to x-ray microscopy [13].However, this model is less applicable to thick 3D objects since it assumes weakly scattering objects (1st Born approx.)[14] and the numerical implementation for typical XRM objects would be computationally very demanding.Bertilson et al extended the partially coherent illumination model of [11] to thick objects by propagating the electric field through the object, but with a 2D treatment of the wave-propagation process [15].Oton et al studied the image formation process in the incoherent case for a thick 3D object subject to a point spread function [16].
In the present paper we describe a 3D wave-propagation-based model for simulating the image formation process in partially coherent XRM:s with thick 3D objects.The model takes into account relevant parameters of the illumination, the object and the objective.We compare our 3D model with the two computationally feasible thick-object 3D models [15,16] as well as with experiments and conclude that it provides an improved description the imaging formation in partially coherent x-ray microscopes.Such improved modelling can assist in optimizing experimental parameters in soft x-ray microscopes, contribute to improved qualitative interpretation of image data, and, potentially, allow retrieval of quantitative data for, e.g., tomographic reconstruction.

Theoretical background
Figure 1 illustrates the principal experimental arrangement for present operational soft x-ray microscopes.Here the object is illuminated by condenser optics, which images the x-ray source onto the object plane.The x-rays are absorbed and scattered by the object and collected by zone-plate optics which forms the intensity image in the detector plane.It is the task of any image formation model to predict this signal given the object and the microscope parameters.

Projection-based models
In the simplest model for thick objects it is assumed that the image intensity at the detector I D (Mx,My) is formed by rays travelling parallel with the optical axis.Each ray carries an intensity I(x,y,z a ) which is attenuated by the local absorption µ(x,y,z) of the object following Beer-Lambert's law where M is the magnification, and z a and z b define a range outside of which µ(x,y,z) equals 0.
The projection-based model is simple and applicable where the wave-nature of the light can be ignored and the whole object is inside the depth of focus.These are valid assumptions when the blur is negligible, e.g., in classical lens-less projection imaging and in lowresolution microscopy, but are generally not applicable for soft x-ray microscopy, where the optics introduce a resolution limit as well as a limited DOF.Furthermore the projection-based model does not include the effects of coherence in the illumination, which has been shown to have significant influence on the image formation [15].

Point-spread-function-enhanced projection models
In the point-spread-function-enhanced projection methods the modelling of the imaging system is improved by including a defocus-dependent point spread function (PSF).The method was first introduced for soft x-ray microscopy by Oton et al [16].In their interpretation, hereafter referred to as the PSF-enhanced method, the intensity I D (Mx,My) at the detector plane is given by (in our notation), ( where I zi (Mx,My) is the intensity image without an object at the detector, z a and z b are zpositions before and after the object, µ is the object's local absorption, h is the PSF and D(z) is the amount of defocus at z.The PSF-enhanced method has the advantage of being analytic, albeit only for special cases such as where the PSF is constant in the z-direction or the object is considered thin in comparison with the DOF.In the present paper we compare our 3D wave-propagation model (Sect.3) with a numerical version of the PSF-enhanced method with a z-dependent PSF.

3D wave-propagation model
Our image formation model consists of the three steps illustrated in Fig. 1: decomposition of the illumination into plane waves, propagation of each plane wave through the object, and imaging by the objective.As is common practice, we assume that scalar theory is applicable [11,[15][16][17][18]. Also, since the bandwidth of XRM:s is typically small, the first two calculation steps are performed monochromatically with wavelength λ 0 , since the effect of a typical experimental bandwidth would only be noticeable in the imaging step.The model is an extension of the numerical model of Bertilson et al [15], hereafter referred to as the 2D wavepropagation model.The major difference is that in our model all calculation tasks are performed in 3D, while the model of [15] is limited to 2D in the computationally intense wave-propagation step.In the wave propagation step the electric field E k is propagated forward in the z-direction for each incident plane wave E k (z a ).In the imaging step the output E k (z b ) from the wave propagation step is imaged coherently using fourier optics methods to the detector plane.The final partially coherent image is formed by taking the sum over all coherent images weighted with the illumination intensity.

Illumination
XRM:s have different condenser schemes (zone plates, capillary optics, normal incidence mirrors etc.) to collect the light and illuminate the object.Following [10] and [11] our model assumes that the condenser illumination consists of secondary source points.Each source point independently produces an electric field E k , resulting in a plane wave with wave vector k, which is directed from the source point to the object and has the magnitude k 0 = 2π/λ 0 .This method is applicable to any of present condenser schemes.Assuming quasimonochromaticity, it is convenient to decouple the stationary envelope u k from the fast harmonic space-and-time dependent oscillations E k as ( ) where j is the imaginary unit and ω = 2πc/λ 0 is the angular frequency.At the input plane z = z a of the wave-propagation step, the stationary envelope is u k (x,y,z = z a ) = exp(j(k x x + k y y)) for the plane wave with wave vector k = (k x ,k y ,k z ).

Propagation through the sample
The modelling of the propagation through the sample is based on the parabolic wave equation (PWE) [17].In addition to scalar theory and quasi-monochromaticity, we assume that the medium has no free charges and is isotropic and non-dispersive.These are the same assumptions and approximations as are made in Fourier-optics theory [18] and also by previous models of XRM image formation [11,[15][16][17].For each plane wave incident on the sample volume the PWE is solved in the form ( ) where χ is the space-dependent complex electric susceptibility of the object.The plane-wave illumination defines the boundary condition at z = z a .At the x and y boundaries of the calculation volume, we employ the transparent boundary condition (TBC) [19], which allows energy to flow out but not to be reflected.Equation ( 4) is solved numerically to obtain the electric field at the z = z b -plane.The solution is calculated using the finite differential method (FDM) with the 4th order Runge-Kutta (RK4) implementation [20].In each iteration the solution u n to the PWE is advanced with a step size ∆z in the z-direction as , where k 1 -k 4 are intermediate evaluations of ∂u/∂z, n is the iteration step, and TBC n is a term that manipulates the border elements in order to apply the transparent boundary condition.Defining f(u,z) = ∂u/∂z, k 1 -k 4 are calculated as The explicit RK4 method is chosen in favor of an implicit method such as the Crank-Nicholson method used in the 2D wave-propagation method [15].Implicit methods normally have the advantage of higher numerical stability and accuracy, but require that a linear equation system is solved in each propagation step, resulting in long computational times.A general linear equation system may be written Ay = x, where A is an N × N band matrix, y and x are vectors of length N, and N is the number of elements in u n .Solving that equation system requires O(N × B 2 ) operations [21], where B is the bandwidth (the width of the tightest diagonal band that contain all non-zero elements) of A. For the 2D case, where u n is an Melement vector, A will be a tridiagonal (i.e., bandwidth B = 3) matrix with N = M.The number of operations required to solve that system is thus O(N × 3 2 ) = O(N), which has proven itself to be a manageable complexity given the problem size in many calculation tasks [15,22].For the 3D case, where u n is an M × M matrix, A will increase in size to N = M 2 .In addition A will have two additional diagonals with non-zero elements at offset ± M from the main diagonal that represent how the PWE now imposes coupling along an additional dimension.Thus, the bandwidth B in that case equals 2M + 1, which yields a complexity O(N ) on the number of operations required to solve the equation system.Both the increase in the size of the data set (from N = M to N = M 2 ) and the higher order in complexity (from O(N) to O(N 2 )) results in a calculation time that makes it infeasible to solve the PWE using an implicit method for the problem sizes considered in this work.
The RK4-method performs a fixed number of operations for each element in u n at each propagation step, yielding the complexity O(N) on the number of operations per propagation step.With this the PWE is solvable with a reasonable amount of computational time, despite the increase from M to M 2 in the data set size.Finally, we note that we do not observe any unstable behavior for this method and that the step size can be determined by the desired sampling in the z-direction rather than by stability issues.

Imaging by objective
The electric field from each illumination direction that exits the object at z = z b is independently imaged coherently by the objective as described in [18].In short, the electric field E D,k that would be at the detector with no other contributions to the illumination than the illumination point with wave vector k is given by the electric field after the wave-propagation filtered by an aperture function A and magnified with a factor M as where K k is the transmission function defined by the combined action of the object and the objective on the plane wave originating from the condenser with unit strength and wave vector k, and  is the fourier transform.The aperture function takes into account the zone plate's numerical aperture (NA zp ) and the defocus (∆f) relative z b via where ν x and ν y are the spatial frequencies, and f is the focal length of the zone plate.The final intensity image at the detector I D is the incoherent sum of the electric field resulting from each illumination point, 2 0 , where defines the intensity distribution of the plane waves incident on the object.With Eq. ( 9) the image formation can be calculated for all intensity distributions of the condenser illumination where the condenser is incoherently illuminated.In particular an arbitrary NA ill is set by defining a circularly symmetric illumination shape where the largest angle of incidence is arcsin(NA ill ).Thus, the degree of partial coherence characterized by the coherence parameter m can be set arbitrarily.For a more general treatment of partial coherence that also includes the degree of coherence at the condenser we refer to the appendix.However, for transmission x-ray microscopes this effect is small and, thus, not included in the present work.Further possible refinements to our model include relaxing the monochromaticity condition to a narrow-bandwidth illumination condition in a similar manner as in [15], and adapting the aperture function to take into account arbitrary aberrations of the optics in addition to the defocus considered here.

Computational issues
To speed up the calculations and have better control of the memory usage the algorithm was implemented in C/C++ , instead of in an interpreted language.The time requirement on a standard desktop computer when simulating the image formation of a 10 µm × 10 µm × 10 µm object with 10 nm sampling is typically 1 minute/illumination point.We have so far used ~200 illumination points, giving a total time of about 3 h for a full simulation.Note that in a typical simulation task the majority of the time is spent on the wave-propagation step.This step, however, normally only needs to be performed once for any given object.After the first computation the image formation for that object can be quickly reevaluated with different microscope parameters.
All parts of the simulation are accessed through a graphical user interface which provides an overview over all available XRM parameters.This simplifies the setup of the simulation task.The graphical user interface also aids in the design of realistic phantoms with detail and complexity that would otherwise not be manageable.
Finally, we note that the code can be run as a batch job, allowing it to be processed in parallel on a computer cluster.The parallel computing is especially suited for generating a full tomographic data set in the same time as it would take to generate the data for a single sample orientation.

Results and discussion
Here we first compare our 3D wave-propagation model with previous models in order to assure consistency as well as to identify where the results of the 3D-model differs from previous models.Secondly, we compare the models with experiments.It is concluded that our 3D model provides a significantly improved description of the image formation.

Comparison with previous models
For a quantitative comparison with previous models we simulated the XRM image formation with the same basic phantom as object.We used the 2D phantom that was employed in [15] (Fig. 2(a)) and extended it to "extended 2D" (Fig. 2(b)) and "full 3D" (Fig. 2(c)) phantoms.The phantoms are made of mylar discs (Fig. 2(a)), cylinders (Fig. 2(b)) and spheres (Fig. 2(c)) with diameters ranging from 66 nm to 280 nm and concentration ranging from 100% to 78%.The mylar features were placed in a water-containing 5-µm diameter 170-nm thick cylindrical shell with 6.3% mylar concentration, and it was furthermore assumed that the 5µm diameter phantoms were embedded in a 10-µm thick water/amorphous-ice layer.
The simulations were performed with the same microscope parameters as in the aperture matched (m = 1) case in [15]: λ = 2.48 nm, zone plate outermost zone width 30 nm, and condenser NA 0.04; parameters corresponding to the Stockholm laboratory microscope [23].For each simulation a focus series from −9 µm to +9 µm was produced.Figure 2(d) shows the result of the 2D wave-propagation model [15] and Figs.2(e) and 2(f) show the same focus series for our model in the x = 0-plane on the extended 2D cylinder phantom and the full 3D sphere phantom, respectively.Finally the PSF-enhanced method [16] was applied to the full 3D phantom.Here we follow [16] by using the defocus dependent PSF of an ideal lens with the same NA as the zone-plate optics.Figure 2(g) shows the resulting focus series in the x = 0-plane.
Already a visual inspection of the focus series in Figs.2(d)-2(g) provides valuable qualitative information.The similarity of Figs.2(d) and 2(e) indicates that the difference in numerical method between [15] and our model does not influence the result.Comparing Fig. 2(f) with 2(e), however, shows significant differences, indicating that the 2D or extended-2D calculations do not provide a sufficient description of the image formation for full 3D objects.Fig. 2. Comparison of computational models.(a) shows the original 2D phantom of [15] with mylar discs, (b) shows the extended 2D phantom with cylinders, and (c) shows the full 3D phantom with mylar spheres.(d) depicts the result of a defocus series when the algorithm of [15] is applied to the original 2D phantom.(e) shows the defocus-dependent intensity at x = 0 when the present 3D model is applied to the extended 2D phantom.In (f) the same calculation is performed for the full 3D object, indicating significant differences compared to the 2D cases.For comparison, the method of [16] is applied to the full 3D phantom and plotted at x = 0 in (g).Finally, in (h) and (j) the four methods are compared quantitatively.The magenta lines in (h) shows how the intensity varies in the center of a mylar feature as a function of defocus, while the blue, green and red in (i) are y-profiles at different defocus positions."Normalized intensity" here refers to the ratio between the image intensity with and without object in the calculation volume.
For a more quantitative comparison the intensity was plotted along different lines (Figs. 2(h) and 2(i)).Figure 2(h) depicts the quantitative intensity as a function of defocus, along the magenta-colored lines which pass through one of the mylar features.Following the qualitative observations above, we note also an excellent quantitative agreement between the 2D wavepropagation method [15] and our method on the extended 2D phantom.As expected, the quantitative differences for the full 3D-phantom are larger, yielding a significantly shorter DOF for the mylar feature of interest.The PSF-enhanced method [16] exhibit a larger modulation in the intensity than the other methods.Figure 2(i) shows quantitative intensity profiles in the y-direction at three defocus positions given by the blue, green and red lines.Also here our method, when used on the extended 2D object, shows good agreement with [15], while the perturbations decay significantly faster with out of focus position when the full 3D object is used.The PSF-enhanced method again results in a higher modulation in the intensity plots.
We believe that the differences between the PSF-enhanced method and our 3D method are due to contradictory assumptions in [16].The general theoretical framework of the PSFenhanced method is based upon the image formation process being incoherent, but Eq. ( 5) in [16] implicitly makes the assumption that the electrical field only propagates parallel to the optical axis.Consequently, in this step the image formation process is treated as fully coherent since the illumination cone would be infinitesimal (NA ill = 0).

Comparison with experiment
Finally we compare the computational methods with experiments.The experiments were performed at the TXM-U41 beamline at HZB Berlin [2].As test objects we used the silicarich frustules of diatoms, since they have the size of a typical mammalian cell as well as a rich 3D structure with high spatial frequencies.The frustules were prepared from a sample of mixed freshwater diatoms and placed on ordinary TEM grids.The XRM imaging was performed at 510 eV (λ = 2.4 nm) with a 40 nm zone plate, yielding the numerical aperture NA zp = 0.03 for the zone plate.The DOF is 2.7 µm while the size of the object was 14 µm.
The numerical aperture of the illumination is NA ill = 0.02, resulting in a clearly partially coherent imaging case since the coherence parameter m = NA ill /NA zp = 0.67 is well below the limit m = 1 where the imaging can be considered incoherent.To allow for a proper comparison with simulations the actual intensity distribution in the hollow cone illumination was measured.Figure 3 shows the measured brightness from the HZB-TXM capillary condenser.The second step was to develop a phantom similar to the frustule-pair that could be used in the computational models.For this purpose we recorded a full tomographic data set at z = 0 defocus of the frustule-pair and performed a tomographic reconstruction using weighted back-projection assuming there was no DOF problem.This reconstructed volume was combined with a priori knowledge about diatom frustules and electron-microscopy images of detailed substructures in frustules of the same species to define a similar (in structure and material composition) 3D phantom.Figure 4(b) shows the result.The 3D phantom was then used as an object in our wave-propagation model, with microscope parameters identical to the values used in the experiment and an illumination equal to the measured brightness (Fig. 3). Figure 4(c) shows the resulting images from this simulation at defocus positions −5, 0, and +5 µm.
The qualitative agreement between the experimental images of Fig. 4(a) and the corresponding simulation of Fig. 4(c) is excellent.In the experimental data as well as in the simulated data different parts of the object move in and out of focus at different defocus positions.This is consistent with the 2.7 µm DOF.Furthermore we observe considerable fringing at the parts of the object that are out-of-focus.This is especially visible at the highspatial-frequency structures like the pores.Fig. 5. Focus series of a pore in the diatom frustule from the experimental data (a) and our 3D wave-propagation method's simulation of the same pore (b).(c) shows the simulated focus series of the pore using the PSF-enhanced method and (d) shows it using ideal projections.Note that the defocus is relative to the diatom center and not the pore, which is why the more focused images are found at the more negative defocus positions.
In Fig. 5 we make a detailed comparison between the experimentally observed intensity pattern from a 120 nm pore and the results from different computational models that allow 3D modelling.Figure 5(a) shows the detailed experimental focus series of the 120 nm diam pore which is marked by rectangles in Fig. 4(a).Figure 5(b) shows the result of our present 3D wave-propagation model.It provides reasonably similar intensity pattern also at this very detailed level.For comparison we include results from the PSF-enhanced method [16] (Fig. 5(c)) as well as the ideal projection method (Fig. 5(d)).It is clear that ideal projections (Fig. 5(d)) provide images very far from experiment.The focus series from the PSF-enhanced method [16] in Fig. 5(c) improves the situation but still falls short of a proper description of the fringing.

Conclusion
We have presented a 3D simulation model based on wave-propagation for the image formation in soft x-ray microscopes operating, e.g., in the water window.The key difficulties are the limited depth of focus and the partially coherent illumination of these microscopes.Our method is benchmarked against previous image-formation models by comparative simulations on similar phantoms.Furthermore, the method is compared with actual experimental images of diatom frustules measured at the HZB-TXM.
The 3D wave-propagation model shows excellent agreement with the previous 2D wavepropagation model when the object has 2D geometry, indicating that the computationally effective Runge Kutta algorithm implemented here performs equally well as the implicit method used before.However, the discrepancy compared with a simulation performed on a 3D object demonstrates the necessity for a model operating in 3D.When comparing the results from different simulation methods with the actual experimental XRM images, our method produced the most accurate predictions of the image formation.For clarity we have used strongly scattering objects in this demonstration.
The 3D wave-propagation method presented here provides an improved and more complete understanding of the image formation in transmission soft x-ray microscopes than previous models.This shows promise for several important applications.Parallel modelling and experiments will allow for a reduced ambiguity in the interpretation of experimental images.In addition, by modelling an experiment before the actual experiment, data acquisition parameters and illumination properties can be chosen for optimal results.In the future, we plan to include the model in an iterative 3D tomographic reconstruction algorithm to improve the results of such 3D imaging.
A special case occurs when the coherence length is shorter than the sampling distance at the condenser, i.e., µ(P k ,P h ) = δ kh .Equation ( 14) then reduces to the sum in Eq. ( 9) used in this work.

Fig. 1 .
Fig.1.Schematic representation of the simulation model.The illumination is decomposed to plane waves that are individually defined by their wave vector k.In the wave propagation step the electric field E k is propagated forward in the z-direction for each incident plane wave E k (z a ).In the imaging step the output E k (z b ) from the wave propagation step is imaged coherently using fourier optics methods to the detector plane.The final partially coherent image is formed by taking the sum over all coherent images weighted with the illumination intensity.

Fig. 3 .
Fig. 3. Measured brightness after the capillary condenser and the central stop at the XRM at HZB.The areas outside of the dashed line were outside the field of view during the measurement and their contents have been extrapolated.θ x and θ y are the angles the illumination makes with the optical axis in the x-and y-direction.

Figure 4 (
Figure 4(a) shows images of two overlapping frustules with different defocus (−5, 0, and +5 µm).The effect of a limited DOF in combination with the partially coherent illumination is obvious, especially at the pores which exhibit significant fringing.

Fig. 4 .
Fig. 4. (a) Experimentally acquired images on two overlapping diatom frustules for defocus position: −5, 0 and +5 µm.The red rectangles indicate the pore which is studied in greater detail in Fig. 5.(b) Renderings of a phantom designed to have similar coarse and fine structure and material composition as part of the frustule-pair used in the experiments.(c) Images from simulations of the image formation using our 3D wave-propagation model.