Impact of local electrostatic field rearrangement on field ionization

Field ion microscopy allows for direct imaging of surfaces with true atomic resolution. The high charge density distribution on the surface generates an intense electric field that can induce ionization of gas atoms. We investigate the dynamic nature of the charge and the consequent electrostatic field redistribution following the departure of atoms initially constituting the surface in the form of an ion, a process known as field evaporation. We report on a new algorithm for image processing and tracking of individual atoms on the specimen surface enabling quantitative assessment of shifts in the imaged atomic positions. By combining experimental investigations with molecular dynamics simulations, which include the full electric charge, we confirm that change is directly associated with the rearrangement of the electrostatic field that modifies the imaging gas ionization zone. We derive important considerations for future developments of data reconstruction in 3D field ion microscopy, in particular for precise quantification of lattice strains and characterization of crystalline defects at the atomic scale.


Introduction
Field ion microscopy (FIM) was the first technique to allow for the direct imaging of individual atoms at the surface of a sharp needle-shaped specimen [1]. Following the pioneering work of Müller, the unrivaled spatial resolution of the technique was exploited to investigate the details of interfaces [2], grain boundaries [3] and adatom diffusion on surfaces [4,5].
A key capability of FIM is its intrinsic capacity to observe defects at the surface of metals, allowing individual dislocations and vacancies to be imaged with atomistic resolution [6]. This ability was exploited e.g. in the investigation of radiation damage in a range of solid materials [7].
There is a strong interest in determining the structure of materials at the atomic-scale, because of the intrinsic microstructure-material property relationship. Therefore, comprehensive models explaining the link between atomic configurations and materials' properties are critical. Atom probe tomography (APT), a successor of FIM, has played a major role in the characterization of materials at the atomic sale in three dimensions.

Impact of local electrostatic field rearrangement on field ionization
APT has been instrumental in this endeavor to link atomic scale structure to the physical properties [8][9][10][11][12]. However, directly imaging all constituent atoms within a material, including vacancies and other types of defects, in 3D is still beyond the resolution limit of all current conventional microscopy techniques, including APT, and is still a daunting endeavor.
FIM relies on the ionization of inert gas atoms from protruding positions on the specimen surface, exploiting an intense electric field. Although each FIM image is intrinsically 2D, pulsed sequential evaporation of surface atoms as in APT, results in a series of images, each representing a 'snapshot' of the surface during the field evaporation process. As the depth of the specimen is investigated by the successive removal of the constituent atoms, the 2D images progressively provide 3D atomic information on the specimen [13]. Vurpillot et al [14] demonstrated a simple approach, whereby a series of digital field ion micrographs can be stacked together, enabling a direct visualization of the crystalline lattice as well as lattice defects. This simple approach was used to image and quantitatively analyze dopants in Si [9]. Unfortunately, this approach lacks the atomic resolution achieved by a previous method proposed by Beavan et al [15], where an observed structure is manually reconstructed atom-by-atom from FIM images. Recent work by Dagan and co-workers [16][17][18] aimed to combine these approaches by applying modern image analysis, data extraction and clustering methods to automate an atomby-atom reconstruction.
Application of these modern techniques enables the imaged positions of individual atoms to be digitally tracked across a sequence of FIM images, from the moment when they are first revealed at the surface of the specimen until their time of evaporation. Analysis of the images has shown apparent small shifts in the imaged atomic positions. While atoms appear to continually slightly change position between subsequent images, a more prominent shift was detected in the local neighborhood around the former location of an evaporated atom directly after its removal from the surface [16]. The magnitude of the displacements of the neighboring atoms left behind on the surface was found to correlate with their proximity to the evaporated atom.
We report an in-depth investigation into the mechanism of image formation in FIM, related to this phenomenon. This is the cause of intrinsic uncertainties in the imaged atomic positions that limit the technique's precision. First, we describe a new method to extract information from sequences of field ion micrographs. Second, we compare our observations with the results of molecular dynamics simulations that include the effects of high charge density and electrostatic field distributions derived from the works of Parviainen and Djurabekova [19] and Rolland et al [20]. We conclude that the observed apparent displacements of atoms by FIM are related to a rearrangement of the charges at the specimen surface, which causes a change in the field direction above the most protruding atoms. This change in field direction affects the trajectories of imaging gas ions, which in effect mimics atomic motion on the surface in FIM images after a field evaporation event. As these field effects limit the accuracy of FIM 3D mapping, they are also likely to affect the spatial resolution of its twin technique, atom probe tomography.

3D field ion microscopy (3DFIM)
In a FIM experiment, the sample is placed inside a cryogenically cooled vacuum chamber, into which a small amount (∼3 × 10 −5 mbar) of inert imaging gas, such as neon or helium, is admitted. Several kV of high voltage is then applied to the sample. FIM samples are produced in the shapes of very fine needles, with an apex radius of only tens of nanometers. Due to its sharpness, an intense electric field is generated at the apex of the needle, causing the polarization of imaging gas atoms. These atoms are then attracted towards the surface of the sample due to polarization forces. When the voltage is raised enough to create a sufficiently high field, the gas atoms are subsequently ionized in a process termed as field ionization. Once ionized, they are accelerated by the electric field away from the tip towards a microchannel plate/phosphor screen, creating a highly magnified image of the sample surface, as seen in figure 1(a). Each bright spot in figure 1(a) is the result of tens of thousands of imaging gas ions, projected onto the FIM screen, ionized above a variety of protruding atomic sites on the specimen surface. The pole structure evident in the figure is characteristic of a body-centered cubic (BCC) lattice.
In the 3DFIM experiment, as the voltage is further raised, surface atoms are removed in a process of field evaporation. As these atoms are removed, the next inner crystallographic layers of the sample are exposed and imaged, thus revealing 3D information. A segment of the evaporation sequence of a single atomic layer is seen in figure 1(b), focusing on a closeup view of the atomically resolved (2 2 2) pole. The terrace patterns corresponding to the crystallographic poles in the FIM image in figure 1(a), including the (2 2 2) pole shown in figure 1(b), appear bright because the locally higher electric field at the positions of atoms protruding most from the surface at kink sites. These sites will generate higher imaging gas currents, as opposed to inner plane atoms, and the brightness will, therefore, be higher. The higher electrostatic field above these atomic sites is also the reason why they will be preferentially field evaporated, enabling the inner plane atoms to be subsequently imaged more and more intensely. Figure 1 3 show the evaporation of one (2 2 2) layer, finally exposing the underlying (2 2 2) layer in figure 1(b) t 4 . Figure 1(b) t 3 shows the last three central atoms of the upper (2 2 2) layer. Imaged on the same FIM image are the outer terrace atoms of the (2 2 2) layer below.

Experimental details and procedure
A tungsten wire oriented along the [0 1 1] direction was shaped into a FIM specimen by electrochemical polishing at 5-8 V using AC with a 5% molar NaOH solution. The final specimen apex radius was estimated to be approximately 20 nm by the ring counting method [21]. FIM experiments were performed on a 3DAP-LAR FIM [22]. FIM images were captured on a phosphor screen and imaged with a CCD camera (AVT Stingray) at a rate 15 fps, with a resolution of 1280 pixels × 960 pixels and averaging every 8 images to optimize the signal-to-noise ratio (SNR mode).
Voltage pulsing is done to induce controlled evaporation. A voltage pulse on top of a DC voltage is supplied, such that total field on the tip is higher than the field necessary to stimulate field evaporation for tungsten. This pulse can be repeated at a certain rate in order to achieve controlled evaporation. The DC voltage was determined according to the best imaging conditions for tungsten with helium [23], and was raised throughout the measurement as the radius of the specimen increased, to maintain the field strength required for imaging. The measurements were performed at a temperature of 50 K. The temperature is measured at the sample stage. Ultra-high purity helium was used as the imaging gas. We maintained controlled evaporation using a voltage pulser set at a repetition rate of 100 kHz and an amplitude that was 25% of the immediate DC voltage. We continuously recorded the images of the arrangement of atoms on the evolving surface of the sample.

Atom identification and tracking
Building on the work of Dagan et al [17] a new automated approach, was developed here to extract data from a 3D FIM experiment. A FIM image can be considered as a sum of Gaussian distributions of intensities centered around the atomic sites. A Laplacian filter was used on the image to highlight the regions with rapid intensity changes. A Gaussian smoothing filter was used to reduce the sensitivity of the Laplacian filter to noise. In order to detect the atomic coordinates an edge detection algorithm based on finding the zero-crossings in the second derivative of the image [24] was used, as demonstrated in figures 2(a)-(c). After these regions were identified, the position of the atom was assigned by the region's centroid.
The entropy of the image was used to classify the identified atoms to their corresponding terraces. Entropy of the image is a measure of similarity of information contained in an image and can be calculated according to the procedure described in [25]. The entropy (H) of the image is given as where k is the number of gray levels and p k is the probability associated with gray level k. Since the atoms at the ledge sites are imaged similarly the pixels in the images atoms are associated with a similar p k . And ledge atoms appear much brighter on the image than compared to the atoms sitting in the terrace's interior. This calculated entropy was used to distinguish the ledge from the terrace, which makes it possible to classify the atoms to their corresponding terrace. The entropy based classification is shown in figures 2(d) and (e). As the atoms are observed multiple times in successive FIM images until their evaporation, a hierarchical clustering algorithm [26][27][28] was used to track the identified positions of an atom throughout the sequence of images in which it appears. An agglomerative hierarchical clustering, with a single linkage criterion based on the Euclidean distances between the imaged positions was used to cluster the multiple images of the same atom. Figure 2(f) shows the result of employing such a clustering. A sequence of images corresponding to the field evaporation of a (2 2 2) pole were analyzed using the proposed routine.

2.3.1.
Quantifying the imaged displacements. The focus of this paper is to study the origin of possible displacements in atomic surface positions following evaporative events, hence results presented in this study characterize an atomicallyresolved and potentially defect-free region of the specimen. To quantify the effects of an evaporating atom on neighboring atoms on the same plane, the positions of all atoms were recorded on all images as exemplified in figure 2(f), tracking the displacements of atoms after the rate-controlled evaporation events. The image sequence of the evaporation along (2 2 2) was divided into three intervals, corresponding to the evaporation of the 7th atom from the terrace having 7 atoms, the evaporation of the 6th atom, and the almost simultaneous, or in a time interval not resolved by the experiment, evaporation of the 5th and 4th atoms. Each interval represents a series of images starting with the first image recorded after the appearance of the terrace configuration, and ending with the first image recorded after the evaporation events described.
Furthermore, three atomic configurations of the last remaining 6 atoms of three different (2 2 2) planes were analyzed to study the influence of the instantaneous atomic configuration and geometry on the measured displacements. It is important to note that a random error estimate cannot be given for these displacements as each figure corresponds to a different state of the tip, and recreating experiments with the exact same conditions is not possible.
We define two specific displacements: first the 'displacement after evaporation', the position difference measured between the first image of the terrace configuration and the first image after the evaporation event, and second the 'average displacement without in-plane evaporation' defined as displacements measured between consecutive images when an in-plane field evaporation event has not occurred. Both of these quantities were measured.

Atomistic simulations
In order to better understand the experimental results, complementary molecular dynamics simulations were employed to obtain detailed information about the atomic motion and the atomic-scale changes in the local electric field taking place due to evaporation events. For these simulations we used the HELMOD [29] hybrid electrodynamics-molecular dynamics code. The code dynamically calculates the induced partial charge on surface atoms, and the corresponding local electric field and electrostatic forces, on the sample due to the applied voltage, taking into account the current local surface geometry. The interaction between the charged surface atoms and the local electric field is included by calculating the resulting Lorentz force and adding it to the standard classical MD inter-atomic interaction. The atomic motion is then integrated using this total force. This approach allows the study of several important effects, such as atomic motion on surfaces, surface reconstructions, etc. However, the cryogenic temperature of the studied system, together with the limited computational resources, effectively limit the study to atomic relaxation and atomic vibrations (less than 1 ns d −1 on a typical quad-core workstation).
In this work the inter-atomic potential used is the EAM W potential by Marinica et al [30]. In the EAM formalism, the binding energy of an atom is given by the sum of pair potential terms, which makes it computationally simple and efficient, but means that it lacks the contribution of e.g. angular terms. Additionally, the chosen potential has not been optimized for surface effects, meaning that e.g. the binding energy of surface atoms may be somewhat inaccurate. However, as the potential has been confirmed to work well with various defect types [30], it is nevertheless a good choice in the absence of specifically surface-optimized W potentials.
HELMOD was modified to use the Robin-Rolland model [20] for charge and field calculation. In this model the partial charge induced on atoms is iteratively computed by finding the equilibrium charge distribution on a surface. The electric field is then only computed at the surface and along the investigated trajectory, instead of the entire vacuum, significantly reducing the memory consumption for the simulation as compared with the grid-based approach employed in previous works [19,29,31,32].
The simulation was carried out on a BCC W sample in the shape of a spherical cap, oriented along the [0 1 1]-direction, with a radius of curvature of 20 nm, and height of 16 nm (figure 3). On the (2 2 2) pole, atoms from the top layer were removed, leaving only the atoms matching the configurations observed in experiments. The bottom of the structure was fixed, while open boundary conditions were used elsewhere. Due to the truncated shape of the simulated specimen, geometric field enhancement occurs at the bottom corners of the structure. However, as this part was held motionless, and the region of interest was sufficiently far away from the corners, this did not affect the results.
The sample temperature was kept at 50 K using Berendsen temperature control [33]. In order to avoid creating surface oscillations, the applied voltage was increased linearly over 73 ps, starting from 0 V until reaching a value corresponding to a surface electric field of approximately 45 V nm −1 , after which the sample was equilibrated for 25 ps at a constant applied voltage. The final surface electric field was chosen to correspond to the known best image field (BIF) of W [23,34] as it should be reasonably close to the surface field in experiments.
Following this, 'virtual' He + ions, were placed at a distance of 3.5 Å [23] from the surface atoms of interest (in the direction of the surface normal, where the highest field is found) and their trajectories in the electric field were integrated using the velocity Verlet algorithm [35] until a distance of 53 mm from the sample surface, similar to the distance to the detector in experiments. All ions were assumed to be singly ionized and have no initial velocity. In order to simplify the calcul ations, these ions do not interact with the rest of the system, i.e. the trajectory calculation is not coupled to the interatomic interactions, and also do not affect e.g. the local electric field. The trajectories of these ions were calculated every 34 fs (corre sponding to an evaporation rate of 300 GHz) during a few picoseconds. While this rate is much higher than in experiments, it makes it possible to examine the deviations introduced by atomic vibrations, by collecting sufficient statistics, while still limiting the computational time required by the simulations.

Results and discussion
We initially measure and analyze the displacements of imaged positions on a (2 2 2) pole when evaporation proceeds from 6 to 5 atoms. Figure 4 shows FIM images before and after a field evaporation event in column (a) and (b) respectively, and the displacements measured from experimental impact positions in column (c). Three different configurations of the same (2 2 2) pole were examined, each reported in an individual row in figure 4, to study not only the magnitude of displacements but also the directions of the changes. To allow for direct comparison of the measured experimental displacements with those predicted by the simulations, the image coordinates were converted to an absolute length scale based on the resolution (pixels × pixels) of the captured FIM images and the physical dimensions of the active detector area. Column (d) of figure 4 shows the simulated impact positions on the screen for the corresponding experimental configurations, as explained in section 2. 4. Column (e) shows a 2-dimensional view of the predicted real space atomic displacements of the atoms in the simulation cell following the field evaporative event. A comparison of the measured impact position's displacements to the simulated impact position's displacements show that they agree well. This indicates that it is the change in the direction of the surface normal of highest field-which manifests itself as a change in the imaged positions-rather than atoms themselves relaxing post-evaporation. To further prove this, we also analysed the cases where the evaporation proceeds from 7 atoms to 6 atoms on the surface and then also for 5 surface atoms to 3, shown in figure 5. The same trend as in the first case is also observed in these cases. Columns (c) and (d) in the figures 4 and 5 show that the measured impact positions and the simulated impact positions show a similar trend. Especially the nearest neighbor's displacement matches quite well with simulated impact positions of the nearest neighbors to an evaporation event. This strongly suggests that it is the change in the direction of the most intense electrostatic field line that manifest as apparent displacements.
Evaporation usually starts from the ledge of a outermost terrace and proceeds towards its center. As each terrace evaporates the lower terraces are exposed to the surface. Apart from the top terrace, atoms from terraces below also keep evaporating in the same order as mentioned above. These evaporation events occurring on the terraces (excluding the outermost terrace) are known as 'out-of-plane' evaporation events. The measured impact positions in experiments should also be affected by the field evaporation of atoms not only exposed at the surface but also due to out-of-plane evaporation events. To illustrate this, figure 6 shows images of the displacements of imaged positions with and without field evaporation from the second terrace, whilst there is no atom field evaporating from the first terrace. Atoms marked with 0 displacements in column (c) of figure 6, start showing displacements after adjacent atoms from the second terrace, labeled with blue arrows, have evaporated.
Further, the overall magnitude of the displacements of remaining atoms on the terrace appears to increase with increasing proximity to the evaporative event, with the nearest neighbors showing maximum displacement. This can be seen in figures 4 and 5 column (c), where the highest displaced imaged positions are the first nearest neighbors of the evaporated atom. Additionally, the displacement magnitude also depends on the number of atoms that evaporate from the plane. The displacements from measured impact positions become larger, as seen in the examination of figure 5 column (c) and (d) for the case of 5 to 3 atoms, when two atoms have evaporated from the terrace. The nearest neighbor's displacement is 0.57 mm (measured) and 0.41 mm (simulations).
Depending on the electrostatic field conditions and the type of surface (chemical composition, surface crystallographic plane, specific surface motif, vicinity to terrace edge), the ioniz ation of the imaging gas occurs at a certain distance above the surface, known as the critical distance. For the current scenario, W (2 2 2) plane at 50 K, with He as the imaging gas, the critical distance is about 3.5-4 Å [23]. The surface normal to the atom along the maximum electrostatic field direction intersects this critical surface at a 'critical point' [36]. This is interpreted as the atom position or center. From the theory of field ionization, it is known that the gas atoms through a series of hops become thermally accommodated with the tip surface. In this process, they make multiple passes through the critical surface. As the probability of ionization is highest at the critical point, most gas atoms ionize there. When an atom field evaporates, the partial charges of the remaining atoms rearrange which causes a change in the local electrostatic field direction thus changing the location of critical points within the critical surface. The changed field direction also effects the trajectories of imaging ions. This manifests as the change in imaged positions of the remaining atoms on the plane. This is schematically shown in figure 7; an atom evaporates, the critical distance decreases, the resolution changes and the critical point or the maximum field direction changes.
The molecular dynamics results made it possible to discriminate between two plausible explanations for the observed displacement of imaged atom positions after a field evaporation event: the actual displacement of atoms, or a change in electric field direction. In figures 4 and 5, columns (d) and (e) present simulated impact positions and simulated atomic  movements respectively. As shown before, all trends identified in the experimental results are reproduced in column d. Simulations also show displacements being larger in magnitude with closer proximity to the evaporative event, and with number of simultaneous evaporation events. In other words, while tungsten atoms do slightly shift their positions on the surface of the sample following field evaporation events, the displacements are smaller over an order of magnitude than deduced from the measured impact positions. While this is now confirmed to be contributing to the observed effect, it is also found to be small. A more significant effect is found to be the field re-distribution after each evaporation event. The true contribution of the atomic relaxation to the observed phenom enon can also be shown as insignificant, by observing the change in positions due to field evaporation events from other terraces.
Field evaporation of atoms from other terraces cause a similar displacement as the ones seen in higher order neighbors due to an in-plane field evaporation event. In figure 6 the displacements shown in second row, after field evaporation from the second terrace are comparable with figures 4 and 5 column (c), higher order neighbor displacements. As the measurements shown here are due to field evaporation from other terraces, the true atomic relaxation would be even smaller, thus making field redistribution argument a major contributor to the observed effect. The electrostatic field redistribution is due to change in the partial charges of the remaining atoms. It is known that under high electric fields the surface atoms are either depleted or are in excess of charge densities [37]. The amount of charge density distribution is dependent not only on the external field but also on its nearest neighbors. The removal of an atom thus causes the partial charges to be rearranged and changes the maximum electrostatic field direction. Hence, the critical point displacement manifests as change in the imaged positions of the remaining atoms after a field evaporation event and the observed effect can be described very compellingly with the local electric field rearrangement.
We also briefly considered other viable explanations for the observed phenomenon. Previous FIM [38] and field desorption studies [39] have indicated atomic movement of the field evaporated atom prior to field evaporation, often termed as a rollup event. In [39] and references therein, a broad review about this type of phenomena can be found. Discrepancies between FIM and desorption images of the (1 1 0) pole in tungsten were explained by these roll-up motion that kink ions perform prior to their field evaporation, assuming a short-lived metastable position. Unlike those works, this study focused on the displacements of nearest neighbors to the evaporated atom. And such roll up events are too short-lived to be observed in FIM. There is another possible phenomenon that could explain the observed effect. For instance, the imaging gas statistics that contribute to every spot on the image are expected to be a factor, as each spot is generated by a large, yet finite number of helium atoms ionizing from the vicinity of the same surface atom. Other than statistics, it is also important to acknowledge that the precise theoretical field adsorption behavior of the gas is not fully understood. Multiple sites for gas adsorption representing the same atom, as well as the presence of helium atoms above the directly adsorbed layer could potentially influence the final imaged coordinates of atoms [40]. The displacements in the imaged positions can originate from the arrangement of the imaging gas on the surface, i.e. the actual configuration of atoms on the surface remains the same, however, the image is changing. Imaging gas atoms adsorbed on atomic sites on the surface will have induced dipole moments, repelling one another, that could cause the adsorbed gas atoms to be slightly displaced outward of the center of the pole. It is possible that as one of the atoms evaporates, gas atoms adsorbed on the remaining surface atoms will relax their positions, and may result in such displacements. While the surface binding energy of tungsten atoms at kink sites is 8.6 eV [41] (likely to change as atoms are field evaporating), the desorption energy of a field adsorbed helium atom is 128 meV [42]. These numbers suggest that the tungsten atoms are kept in their place more strongly than the imaging gas atoms and so could support the case that these displacements are indeed caused by image gas relaxation. But this explanation seems unlikely because the imaging is primarily due to ionization of thermally accommodated gas atoms rather than the desorption of an adsorbed gas layer. However, we cannot completely rule out the interaction of the adsorbed gas layer and the thermal accommodation of gas atoms. This would require a thorough quantum mechanical treatment of the phenomenon which is outside the scope of the present article. However, by all the numerical demonstrations shown above, a more significant, larger scale effect is occurring with each field evaporation event which is the electrostatic field redistribution.

Conclusions
We have shown that experimentally observed displacements correspond to the change in the direction of highest electric field that result in a change in the position of the critical point within the critical surface, rather than to a true relaxation of the atoms on the surface. This was done by comparing experimental results with high electric field molecular dynamics simulations. These calculations show that the corresponding true atomic relaxations after field evaporation are typically less than 4 pm, much less than the displacements measured from the FIM images. The imaged positions of atoms in the first terrace are also shown to be affected by field evaporation taking place in the terrace below. This finding further affirms that the observed atomic displacements in FIM are due to the field redistribution rather than due to atomic relaxations.
This work has provided more insights into the physics of image formation in FIM. Some more inferences can be drawn with regards to data extraction from FIM images. It is shown that the critical point moves away from the surface normal of the true atomic center as the surrounding atoms field evaporate. Hence, the more intense the imaged atom is, the farther it is from its true atomic position. Hence, the least erroneous positions of the atoms would be the ones which are extracted from the atoms imaged in the middle of terrace. Although not much can be done in terms of resolution of FIM, this study helps to improve the precision of the FIM. This will lead to more precise measurement of atomic strains around the lattice defects, when using FIM, for instance.
To summarize: • The apparent motion of imaged positions of atoms correspond to the change in the position of critical point on the critical surface, rather than relaxation of the atom itself. • Simulations show that the corresponding actual atomic relaxations due to evaporation are quite small (less than 4 pm) to be mirrored in the imaged positions of atoms in FIM. • The imaged atom positions in the first terrace are also shown to be affected due to the evaporation happening in the terrace below. • The least erroneous positions of the atoms would be those extracted from the atoms imaged in the middle of terrace. • Applying these corrections would lead to more precise measurement of atomic strains around the lattice defects.