Ab initio Description of Bond-Breaking in Large Electric Fields

Strong (10 V/m) electric fields capable of inducing atomic bond-breaking represent a powerful tool for surface chemistry. However, their exact effects are difficult to predict due to a lack of suitable tools to probe their associated atomic-scale mechanisms. Here we introduce a generalized dipole correction for charged repeated-slab models that controls the electric field on both sides of the slab, thereby enabling direct theoretical treatment of field-induced bond-breaking events. As a prototype application, we consider field evaporation from a kinked W surface. We reveal two qualitatively different desorption mechanisms that can be selected by the magnitude of the applied field.

The breaking of an atomic bond is one of the most fundamental phenomena governing materials transformation, reaction, and degradation. Phase changes, mechanical deformation, chemical reactions, corrosion, and many other important processes can be understood in very simple terms as a systematic and often coordinated sequence of bond-breaking events. Probing and controlling these processes is therefore only possible with a clear understanding of the underlying effects that stimulate bondbreaking.
Of the possible stimuli for bond-breaking, electric fields are among the most ubiquitous. A local 10 10 V/m electric field is of the same magnitude as the intra-atomic fields between electrons and nuclei [1] and is therefore perfectly capable of severing atomic bonds. Because the field at a material's surface scales inversely with the local radius of curvature, even moderate voltages can be locally enhanced into fields of this magnitude anywhere that sharp features exist, such as surface steps and kinks. [2] This field enhancement enables atom probe tomography (APT), a microscopy technique wherein nanosharp material samples are intentionally evaporated under strong fields. Ionized atoms that evaporate from the sample's surface are later collected at a counterelectrode ( Figure 1). [3][4][5] After the evaporation, the sample is computationally reconstructed by back-tracing each ion's trajectory using its time-of-flight and detected location at the counterelectrode. The accuracy of these projected trajectories, and consequently the accuracy of the image reconstruction, depends on our understanding of the mechanisms by which the original surface bonds were broken.
Density functional theory (DFT) calculations, which could enable a direct investigation of evaporation mechanisms, are hindered by the challenge of applying a finite electric field under periodic boundary conditions. [6] Here, we report on an efficient solution for this problem within a framework of standard DFT calculations, and demonstrate its usefulness for elucidating the complex evaporation mechanisms from prototypical kink sites on high-index surfaces.
Under three-dimensional periodic boundary conditions, a surface must be modeled as a two-dimensional slab; i.e. the model system will have a surface on either side. This slab must be sufficiently thick and have enough vacuum above and below it to prevent these two surfaces from artificially interacting, and any applied electric field must be accounted for when it crosses the boundary (Figure 2).
Previous DFT studies have accounted for the electric field by introducing a sheet charge in the vacuum as the counterelectrode in an overall charge-neutral setup, and enforcing symmetry along the z-axis to ensure that the  [7] with opposite fields on either side of the slab, c) the "dipole correction", [8] with a constant field across the cell, and d) the generalized dipole correction (this work), with a field-exposed and field-free side. The resulting potential φ (red) and field E (blue) are shown for each scheme.
surface charge on either side of the slab is well-controlled ( Figure 2 (b)). [6,7,9] Using this approach, the slab structure must not only be strictly symmetric, but also sufficiently thick to converge the potential and Friedel oscillations beginning from either side of the slab. The cell's vacuum region must also be enlarged to mitigate the artificial Coulomb repulsion between the two evaporating ions on either side of the slab. These constraints reduce computational efficiency and restrict the approach's feasibility to simple cases, such as an adatom evaporating from a flat surface. Experimentally, however, the most relevant sites from which field evaporation occurs are kink sites at the edges of terraces on the round emitter surface where the local curvature induces strong field enhancement. [2,10] In order to enable calculations for large surfaces that contain such low-coordinated features, like the (10 8 1) tungsten surface shown in Figure 1, the electric field's periodicity must be accounted for in a more general way. A starting point is the well-known "dipole correction", [8] in which an infinitesimally thin dipole sheet is added to create a discontinuous potential jump in the vacuum region of a DFT cell (Figure 2(c)). The magnitude of this dipole is chosen such that it exactly compensates the dipole of the slab, creating a constant-field condition even for asymmetric slabs. This formalism has been extended to introduce an additional finite field, [11] but above a critical field strength the vacuum potential is pulled below the Fermi level. This results in the spurious transfer of electrons into the vacuum (Figure 2(c)). [12] They key concept we propose in this letter is to augment the dipole layer with a charged monopole sheet. This charged layer acts as a counterelectrode in the vacuum, creating a discontinuous jump not only in the potential but also in the field. The result is a constant positive field between the counterelectrode and one side of the slab while the opposite side remains field-free ( Figure 2 (d)). The dipole correction must then compensate the dipole of the combined system (slab + counterelectrode). This approach, which we term the generalized dipole correction (GDC), may be conceptualized as a combination of the counterelectrode and dipole correction schemes, as it introduces a symmetric charge compensation without requiring that the slab itself be symmetric. This combination keeps the advantages of the respective approaches, while eliminating their respective disadvantages.
The generalized dipole correction that must be added to the standard electrostatic potential with periodic boundary conditions for a slab with charge Q reads (1) Here, z 0 is the cut position within the vacuum, A and c are the surface area and height of the slab supercell, respectively, and V 0 is an offset that brings the planeaveraged total potential V at z 0 to a constant value (we use V (z 0 ) = −E top z 0 where E top is the electric field on the top side of the slab). The z 2 term compensates for the implicit homogeneous background in the periodic potential, while the correction field is where µ is the charged slab's dipole moment with respect to z 0 , and E bottom the field on the bottom side, which is zero for the field evaporation calculations performed in this work. If the above V (z 0 )-alignment is used, a consistent electrostatic energy can be directly obtained from the total potential and the total charge density ρ es (including nuclear charges) as The GDC has been implemented in our DFT code SPHInX, [13] allowing us to directly investigate evaporation mechanisms from experimentally relevant surface sites using DFT.
To demonstrate the performance and applicability of the GDC approach, we consider evaporation from a kinked tungsten surface. The field-dependence of the activation energy for evaporation events in tungsten do not follow the behavior predicted by basic theoretical models, [14] which generally assume an ideal straight-line departure of the ion from the surface. This discrepancy has prompted the proposal of a number of nontrivial evaporation mechanisms, including possible out-of-sequence evaporation, [15] a roll-up motion of atoms onto neighboring step edges, [16][17][18][19], or diffusion across the surface prior to evaporation. [7] A combination of several of these effects is also possible.
To model field evaporation from this system, we model the kinked surface as six atomic layers with a (10 8 1) surface normal, resulting in a 98 atom structure that has semi close-packed (1 1 0) terraces with a single (0 0 1) step every 7 unit cells and a (1 0 0) kink every 3 unit cells along the step. [20] A representation of this slab, which is given 15Å of vacuum between its periodic images along the z-axis, is shown in Figure 1. All DFT calculations are performed using a local-density approximation (LDA) functional with a 3×3×1 k-point mesh, a 20 Ry energy cutoff, and 0.1 eV Fermi smearing to allow partial electronic occupations. The LDA functional is chosen based on its accurate reproduction of surface energies and work functions [21] and the energy cutoff and k-point resolution give forces converged to within 10 meV/Åof a 7×7×1 k-point mesh and 30 Ry energy cutoff.
We geometrically optimize the surface at a number of field strengths, and then pull an atom sitting at the kink site highlighted in Figure 1 from the surface by incrementing its z-coordinate above its original position. [7] The kink atom is chosen because it is the surface atom with the fewest nearest neighbors and should thus be the most weakly bound. Indeed, these are nearly always the first sites observed to evaporate in experiments. [22] At each incremental height the atom's x-and y-coordinates, as well as the top three layers of the surface, are reoptimized using a quasi-Newton algorithm based on the forces acting on each atom. This enables us to observe how the coordinates of the evaporating atom and the neighboring lattice change during the imposed evaporation event.
For a dense set of field strengths, we compute and plot the total energy as a function of the evaporating atom's z-coordinate, as shown in Figure 3. The discontinuities in the potential energy curves in Figure 3 are the result of discrete changes taking place along the evaporation reaction coordinate. For distances where the evaporating atom is very close to its original position (<1.5Å), the energy follows a smooth bond-stretching trend. There is an abrupt shift to a new minimum between 1.5 -2Å. At this height the evaporating atom shifts laterally into the nearest hollow site atop the neighboring (1 1 0) terrace. At lower heights this motion is sterically prohibited by the neighboring step atoms. The atom's evaporation then proceeds from this new minimum in a manner very similar to an adatom on a flat (1 1 0) surface. The second discontinuity in the energy vs. distance profile occurs when the original atom begins to pull up its nearest neighbor out of the step edge; at certain distances the two even form a dimer above the surface. The persistent interaction between the evaporating W atom and its nearest neighbor is likely responsible for the substantial number of spatially correlated co-evaporation events experimentally observed during tungsten evaporation. [23] At distances sufficiently high above the surface (>3.3Å), the bond between these two atoms becomes too weak to pull the neighboring atom up from its original position. Our calculations therefore reveal that the evaporation mechanism is effectively a two-stage process: a rollover event followed by the actual departure from the surface. Each of the two stages has its own respective energy barrier. The two barrier heights vary quite differently as a function of field, as shown in Figure 4. Since evaporation requires both mechanisms sequentially, the effective barrier observed in experiments belongs to the mechanism with the rate-limiting (i.e. highest) barrier. At low fields, the second step, in which the atom is forced to ionize, has a much higher barrier. However, the barrier in this step strongly decreases with increasing field. Due to this rapid decrease, its barrier drops below that of the rollover stage's barrier at a field strength close to 4.2 V/Å. We note that above 4.2 V/Å, the atom still to the initial rollover motion (blue) and the barrier to evaporation from the adatom site (red) are plotted separately. Experimental data previously reported for tungsten [14] is shown in black. For comparison, the barrier height for a W 2+ ion calculated using the image-hump model is also shown (dashed line). travels through the adatom site before evaporating, but experiences no barrier after the initial rollover motion. From an energetic standpoint, therefore, the mechanism effectively switches from two-stages to one-stage at high fields, although the evaporating atom follows the same pathway as for low fields.
If the local field is below the critical value of 4.2 V/Å, which can be the case on the shank of the emitter where the field is reduced, the kink atom may be thermally stimulated (e.g. by a nearby oncoming laser pulse) to roll over to the adatom position but still not have enough energy to evaporate. In this energetic trap atop the flat (1 1 0) terrace, the lateral hopping barrier for an adatom is reduced from 0.9 eV to around 0.7-0.8 eV by the field. [24,25] The suppression of surface diffusion barriers by the field can be even more pronounced depending on the material and nature of the diffusion mechanism. [26] Any net displacement of the atom's position before it evaporates, whatever the mechanism, is detrimental to the APT reconstruction's accuracy. [27,28] Therefore, experimental conditions should be chosen to avoid or mitigate diffusion from the adatom trap. Figure 1 displays the location of the transition state (where the evaporating atom experiences the barrier) as a function of field strength, showing that the transition state exists very near the atom's original location for fields above 4.2 V/Å. For fields below 4.2 V/Å, transition states for both barriers are shown, including the second one above the adatom site. The strong sensitivity of the location where the transition state exists renders an often employed approximation -that the zero-field barrier configuration can be used for all field strengths -invalid.
The ab initio-calculated barriers can now be compared with those derived from existing models. Historically, one of the most commonly used models to approximate field evaporation is the image-hump model, [29][30][31] which superimposes the field potential and the image potential to determine the barrier height for an ion leaving a flat surface. The main advantage of this model is its simplicity. The only material-dependent parameters to enter the formula for barrier height are the material's sublimation energy and relevant ionization energies, which in most cases can simply be looked up. However, this and related models [32] have been proven to predict severely inaccurate temperature-dependent evaporation fluxes. [14,16,[33][34][35] In Figure 4, for example, we compare the image-hump barrier for the evaporation of a W 2+ ion with experimentamatcheslly determined W evaporation barriers. [14] The model nearly matches the extrapolated critical field of 6.2 V/Å from our calculations, but predicts unphysically high barriers for all other fields. Because APT experiments are generally performed at fields below the critical field limit and often use lasers to thermally stimulate evaporation, [36,37] the model's predictions are invalid for exactly the experimentally most relevant range of fields.
The failure of the image-hump model to describe the experimental data in Figure 4 has led to the recent proposal of several new analytical models to calculate evaporation barriers. [22,38] One of the most commonly accepted is the so-called "charge-draining" model, in which the evaporating atom is considered to continuously donate charge to the slab and gradually ionize as it departs the surface. Previous DFT calculations on charged Al (1 1 1) adatoms support this nature of charge transfer. [7] Due to the continuous ionization, this model yields an evaporation barrier that decays linearly as the field is increased. Unfortunately, these models contain effective parameters that must be obtained from external sources. In practice, the slope of the decay is generally fit empirically to available experimental data. However, since this slope depends directly on the shape and size of the barrier encountered by the evaporating atom, [38,39] potential energy paths calculated using DFT with the GDC approach, as in Figure 3, are a reliable route to provide quantitative accuracy to these more conceptually sophisticated analytical models.
We conclude that conventional APT experiments in tungsten automatically probe the rollover response, which is ultimately detrimental to their 3D spatial resolution. The rollover response can be understood as a competition between the force of the evaporating atom's nearest neighbor bonds and the force of the field tugging on the ion. As a result, softer metals with weaker surface bonds are expected to exhibit a less pronounced version of this effect than what is observed here for tungsten.
The two-stage rollover evaporation mechanism provides a natural explanation for the experimentally observed evaporation barrier versus applied field in Figure 4. It clarifies that at very high fields, experimental evaporation events are dominated by a thermallyactivatable rollover barrier. Of course, observations from APT experiments also depend on several phenomena which occur at length and time scales inaccessible to DFT, including mesoscopic field and temperature gradients. The atomic-scale evaporation mechanism is therefore an important piece of the overall theory of field evaporation in APT, which requires considerations beyond DFT to account for these larger-scale phenomena.
The evidence provided in this study for a fielddependent, tunable evaporation mechanism is essential for accurately controlling and interpreting APT and field ion microscopy experiments on metallic systems. The generalized dipole correction developed here provides a computationally efficient and easily implementable approach to model the effect of strong electric fields in DFT calculations. The correction is universally applicable to other material systems in order to understand bondbreaking mechanisms in more complex materials systems, e.g. aqueous corrosion systems or catalytic surfaces. Using this technique to directly probe the response of materials and chemical reactions, such as bond-breaking, in extreme electric fields will provide a new tool to guide the interpretation and design of new experiments and applications.