Towards model-driven reconstruction in atom probe tomography

Reconstructions in atom probe tomography (APT) are plagued by image distortions arising from changes in the specimen geometry throughout the experiment. The simplistic and inaccurate geometrical assumptions that underpin the conventional reconstruction approach account for much of this distortion. Here we extend our previous work of modelling APT experiments using level set methods to three dimensions (3D). This model is used to generate and subsequently reconstruct synthetic APT datasets from electron tomography (ET) of an Al-Mg-Si multiphase specimen. Finally, we apply our model to the reconstruction of an experimental field-effect transistor (finFET) dataset. This model-driven reconstruction successfully reduces density distortions compared to conventional methods. By combining prior knowledge about the specimen geometry from sources such as ET, such an approach promises new distortion correcting APT reconstruction applicable to complex specimen geometries.


Introduction
Atom probe tomography (APT) is a three-dimensional (3D) mass spectrometry technique for materials characterisation in which surface atoms on a needle shaped specimen are ionised under a pulsed surface electric field via the physical process of field evaporation, and are accelerated towards a two-dimensional (2D) position-sensitive detector. While the time-of-flight chemically identifies the ionic species, the 2D detector position estimates the original positions of ions within the specimen. The detected ions form a specimen subvolume, within which the detected ions original locations are estimated Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. from their detector positions. This process is known as reconstruction.
Current reconstruction algorithms assume ion trajectories obey a mathematical point-projection law from a specimen apex that remains hemispherical throughout the field evaporation process [1][2][3]. However, this assumption is often highly unphysical and a dominant source of error. Image distortions arise in inhomogenous materials, which are often a focus of study [4][5][6]. Here differing material phases typically field evaporate at different rates under an applied electric field. These non-uniform evaporation rates over the specimen surface drive the apex away from the hemispherical geometry assumed by point-projection protocols, changing the surrounding electric field and altering the trajectories taken by field evaporated ions. The result is a break down in the assumptions of the point-projection method. This failure ultimately appears as positioning distortions in APT reconstructions, observable as density fluctuations within both the 2D detected ion positions and the 3D reconstruction.
In order to try and correct for these distortions, improved models for the evolving geometry of the specimen within APT experiments are required. Previous attempts have largely focused on atomic models [7][8][9][10][11], where the simulation grid is roughly inline with the underlying lattice structure, and analytical extensions [12,13]. However, these approaches have proved either computationally intensive preventing fullscale simulation [14], suffer from inherent instabilities, or are geometrically constrained (e.g. applicable only to multilayer specimens) and thus may not be suitable to drive a general and rapid reconstruction algorithm. Simulating a full scale APT experiment (typically on the order of ten million ions) with current atomic models would take 23-34 days given their typical evaporation rates of around 200-300 atoms per minute [10,14]. Our previous work [15][16][17] has demonstrated that level set methods could provide the speed and generality required, but have so far fallen short of fully simulating electrostatics in 3D. Another promising reconstruction approach aims to build up the specimen surface from the detected ions in reverse order [18,19]. However, this method currently fails to capture any concavity of the specimen surface (as seen in this study), does not physically simulate ion trajectories, and is missing information outside of the experiment field-of-view (FOV).
This study extends our 2D model [17] to full 3D and aims to show how this new model can be used to perform rapid modelling of realistic specimen geometries to improve understanding of phenomena in APT experiments, and to directly drive the APT reconstruction procedure.

Theory
Our model-driven reconstruction in APT requires three main components. First the model must track the specimen surface. Next the model must calculate ion projection up to a virtual detector. Finally, a calibration stage is required to match up simulation times with experiment times (detected ion indices). This is achieved by matching the simulated evaporated volume of material within the FOV to the detected ion count in the experiment. This final stage will be covered in section 4.4.

Field evaporation law
Throughout this work an Arrhenius' law is assumed for estimating the local evaporation rate R e , identical to that used in our previous work [17]. This law has been shown to be valid experimentally. The isotropic law without faceting is given by equation (1). (1)

R e (s) = Ae
Here A is a global constant with dimensions of distance over simulation time, β is a temperature dependent parameter Diagram showing how the specimen surface is embedded in the level set field ϕ(x). The gray plane represents the intersection with ϕ(x) = 0 (given by the black curve). This visualisation was generated using our 2D model [17].
ranging from approximately 40 (laser mode) to 400 (voltage mode), ||E(s)|| is the local surface electric field strength, F 0 (s) is the isotropic local evaporation field, and s ∈ Γ defines a point on the specimen surface Γ. As with our previous study [17], all simulations assume β = 40 and the field reduction rate F R , defined in equation (2), was held at F R = 1 by rescaling the voltage applied to the specimen. It is worth noting that in the continuum limit, it appears that such large values for β make almost no difference on the evaporation dynamics compared to using smaller values, with the exponential term in equation (1) effectively behaving similarly to a step function (see section 1 in the supplementary material (available online at stacks.iop.org/JPhysD/53/475303/mmedia)). The evaporation law used in this study makes the approximation that A, F R , and β are constants and independent of phase. Note that the exact value of A used in the simulation does not effect the simulated specimen evolution, but instead corresponds to a global rescaling of the simulation times. The evaporated material within the Field-Of-View throughout the simulation remains unchanged.

Level set interface tracking
The level set method was used to track the specimen geometry by embedding the specimen surface, denoted by Γ, as a closed zero contour within a three-dimensional (3D) scalar field ϕ(x) (see figure 1). This definition is formalised in equation (3): This scalar field is initalised as a signed distance field from the surface (positive within the specimen and negative outside of the specimen) e.g. equation (4): where d(x, Γ) defines the mimimum Euclidean distance between domain point x and the specimen surface Γ.
Here Ω + is the subdomain inside the tracked surface, and Ω − is the subdomain outside of the tracked surface. It is this implicit representation of the specimen surface in the level set method that allows the model to capture surface discontinuities (e.g. facet corners) and internal holes (voids). The level set field is evolved under the level set equation (equation (5)) to propagate the motion of Γ according to a specimen surface velocity function R e (s), defined in equation (1). The level set equation is given by: where ν(x) is the velocity extension constructed from R e (s) over the entire level set grid via a closest point method [20]. Under this extension, the velocity ν at a domain point x is set to the velocity at the closest point on the specimen surface, as given in equation (6):

Solving the electric field
The APT evaporation problem can be treated electrostatically [21] by modelling the electric field surrounding the specimen via Laplace's equation with Dirichlet boundary conditions, given in equation (7): where u(x) is the local electrostatic potential, E(x) = −∇u(x) is the electric field, and x ∈ Ω − is the domain of points outside of the specimen. In this study, the specimen has been approximated as a conductor, corresponding to Dirichlet boundary conditions for the problem defined in equation (8).
In order to solve the electric field a direct constant potential boundary element method is used, similar to our prior work but here extended to 3D. The specimen surface Γ is discretised into k elements or panels Γ j of constant potential u j and normal flux q j . In order to solve the normal flux values on panels, the following system of k linear equations are constructed by enforcing k boundary integral equations at collocation points located at the centre of each panel x i ∈ Γ i . This system of linear equations is given in equation (9) and is indexed by i: where F ij and G ij are elements of the the k × k problem matrices, y defines a point on the boundary element Γ j , and n j is this boundary element's normal. The matrix elements defined in equation (10) and equation (11) are calculated by numerically integrating [22] over surface panels (see figure 2(a)), except for diagonal matrix elements for which slower analytical solutions are required to treat singularities. Unlike in 2D [17], the exact form of these integrals is complicated and can be found in [23]. Under the conductor approximation, the normal flux q i corresponds to the surface electric field strength ||E(Γ i )|| over the boundary element Γ i . The electric field at a particular point outside of the specimen is solved via the discretised hypersingular boundary integral in equation (12): K(z, y,n z ) = (z − y) ·n z ||z − y|| 3 (13) H(z, y,n z ) = − 1 ||z − y|| 3n z ·n j where z ∈ Ω − is the point at which the field is to be evaluated,n z is the unit vector along which to resolve the electric field, andn j is the panel normal. These field points are evaluated as required during the ion trajectory integration. Under the conductor approximation, the integrals in the first term of equation (12) are identically zero e.g. j´Γ j H(z, y,n z )d 2 y = 0 [24]. Therefore, for this study only the term associated with the K(z, y,n z ) kernel needs to be calculated. For every field point evaluation, equation (12) is calculated forn z = (1, 0, 0),n z = (0, 1, 0), andn z = (0, 0, 1). These give the electric field components in Cartesian coordinates E(x) = (E 0 , E 1 , E 2 ). A diagram showing the contribution of a surface panel in evaluating the electric field at a particular z is given in figure 2(b).

Panel clustering.
During the ion projection step, the electric field must be evaluated at many z along each ion trajectory. Practically this involves evaluating equation (12) thousands of times per ion trajectory integration. Solving equation (12) directly requires summing and integrating over every boundary element making up the specimen surface. Such direct evaluation is very likely far too slow for driving a practical reconstruction algorithm.
In order to accelerate the ion projection step, a panel clustering method has been implemented, similar to the method described by Hackbusch et al [25]. Such a method allows equation (12) to be rapidly evaluated by grouping together neighbouring panels, and approximating their contribution to equation (12) with a truncated multipole expansion [26]. Possible candidates for these clusters of panels are determined prior to the ion projection using an octree structure [27]. A function is then calculated, also prior to ion projection, for each cluster candidate that approximates the cluster's panels contribution to equation (12). This approximation holds if the cluster is sufficiently far from the field point being evaluated (currently determined by a tolerance). These approximating functions for the panels in clusters are derived by integrating and summing over a cluster's panels prior to ion projection. This way the approximating function only depends on the field point z. By using these functions, integration and summing over individual panels is only required when panels are close to the field point z, reducing the total number of summations required to evaluate the field at z.
For evaluation of the field at a particular z, appropriate panel clusters are selected depending on the location of z. These cluster candidates have been defined prior to field evaluation during the construction of the octree. Larger clusters can be selected further away from z without a loss of accuracy below the specified tolerance due to improved convergence of the multipole expansion. Panel clusters are selected from the list of precomputed cluster candidates so as to minimise the total number of clusters required. Only those panels close to z and are computed exactly without approximation. This approximation of the field solution for an example panel cluster is outlined in figure 3.
Our implementation proceeds by constructing an octree in space, similar to that in [10]. Boundary elements which have a centre point falling within a particular node of this octree are grouped into a cluster Γ c . The octree terminates once the leaf node cluster size falls below a user-defined number of boundary elements M c . This way every octree node has an associated cluster of boundary elements. This octree spatially groups panels into possible cluster candidates. During the field evaluation step, the octree will be traversed and a particular group of these clusters, which can also be considered a partition of surface panels P Γ , is determined to efficiently calculate the field at z. This grouping obeys the two conditions required by a set partition. The first condition is that every panel in the specimen surface must be in a cluster of the selected grouping. The second condition is that no two clusters in the selected grouping can contain the same panel (they must be disjoint).
Each cluster also has a centre y c and radius R c associated with it. Letting C Γ c be the unique set of corner points for boundary elements within a particular cluster Γ c , these quantities are defined in equation (15) and equation (16): where ν i is a corner point, and N Γ c is the number of boundary elements in cluster Γ c . Using these definitions, y c corresponds to the centroid of cluster panels and R c the minimum sphere radius including all panels in this cluster. For a particular cluster partition, the solution in equation (12) can be split into a number of partial sums associated with particular clusters. Such a splitting is given in equation (17). The aim is to find an approximate solution that holds in the case that a cluster is sufficiently far from the field point z, e.g. ||z − y c || >> R c . Here a first-order truncated multipole expansion is used to approximate the partial sum associated with the cluster, i.e. equation (18). This expansion M(z,n z , y c , R c ) is explicitly derived in appendix A.
Critically for a particular cluster and field direction vector n z , once the expansion coefficients for the cluster have been computed (see appendix A), the multipole series expansion M(z,n z , y c , R c ) is only a function of the field point z. Therefore, this same expansion can be used to compute the field at multiple points z without having to re-evaluate the partial sum and panel integrals in equation (18). To determine the validity of this approximation for a particular field point z, the following admissability condition ||z − y c ||ϵ t > R c should be satisfied, where ε t is a user-defined parameter. Throughout this study, ε t = 0.1 has been used.
The multipole expansion coefficients are calculated prior to ion projection via an upward pass of the octree. This starts at the leaf nodes of the tree and moves upward until reaching the root node. All series expansion integrals, given in appendix A, are solved numerically [22]. Such an operation has a combined computational complexity of O(k · log(k)), where k is the total number of surface boundary elements. This could be theoretically reduced to O(k) through moment-to-moment translations [28].
Once this octree has been constructed and coefficients calculated, the ion projection is undertaken. During this stage, equation (12) must be solved at many positions z ∈ Ω − , determined by the previous step of the trajectory integrator. The aim is to assemble an optimal group of clusters that defines a partition of the surface panels. In order to reduce computation, this partition should minimize the total number of clusters under the constraint that all clusters in the partition meet the admissability condition. For a particular field evaluation at a position z, a downward pass of the octree is performed. The algorithm recursively trials clusters in the octree, starting at the root node (top node). If a trial cluster is admissable, then it is selected for membership of the grouping P Γ . If not admissable, the algorithm proceeds to consider its children nodes for membership of P Γ . This downward pass terminates once all the specimen surface panels are in a selected cluster. Finally, equation (17) is evaluated, using the series expansion approximation in equation (18) for non-leaf clusters. Any leaf node clusters selected for the partition must, by nature of the algorithm, be close to z. Therefore, these leaf node clusters are directly evaluated by summing and integrating over boundary elements within these clusters, i.e. the left-hand side (LHS) of equation (18). This octree-based partitioning of surface panels for a particular field point z is illustrated using a simplified 2D case in figure 4.
The panel clustering approach reduces the computational complexity of the entire ion projection from O(k 2 ) to O(k · log(k)), while increasing the memory complexity from O(k) to O(k · log(k)).

Ion trajectory calculation.
Reconstruction in APT requires a mapping describing how ions at a particular spatial detector position and experiment time should be positioned in the reconstruction. To obtain these mappings ions are launched from specimen surface panel centres at fixed simulation times, trajectories calculated and final impact coordinates recorded on a simulated detector. Detector impact positions are then interpolated with respect to initial launch coordinates (figure 5). By also interpolating between each timestep, a function is constructed that maps particular detector positions and simulation times to particular positions within the simulated specimen volume. Note that only the initial and final points of ion trajectories determine this mapping-the ion behaviour during flight is not explicitly required. For this study, when generating the mappings, one ion was launched from the centre of each panel every 10 iterations.
To simulate APT data, ions were launched from the specimen surface with an ion launch density (ions per unit area of the specimen surface) proportional to the local evaporated volume (local evaporation rate multiplied by the iteration timestep). This density is determined by a user-defined ionic volume parameter. Each ion is labelled by the material phase defined at its launch position.

Specimen crystallography.
To correct for density distortions in reconstructions due to pole structure, our model must be capable of accurately simulating specimen faceting. The method for capturing shape anisotropy due to underlying specimen crystallography in the 2D model [17] can be transferred to the 3D model with almost no modification. Once again a set of facet orientations and evaporation fields can be provided for the particular crystal grains to be modelled where N is the number of major lattice planes present in the crystallographic structure, and gamma is an anisotropic constant. Particular facets are defined by normal vectors f i with a magnitude corresponding to the facet's steady state evaporation field. This allows a local evaporation field to be defined over the specimen surface via equation (19): tan θ min (s) = min Coloured nodes in the tree show the panel clusters used in the evaluation of the particular field point z shown. Panels contained in leaf node clusters are too close to z for the series approximation to be valid. Instead, these panels are directly evaluated by integrating and summing over them. f min (s) = argmin where θ min defines the smallest angle in 3D between the specimen surface normaln, and the closest facet normal vector f min . The closest facet is determined by calculating θ i for each facet and identifying the facet giving the smallest angle. Due to the 3D extension of the model, anisotropy functions are now given by polar plots of surfaces, rather than the polar plots of lines in the 2D model.

Implementation
A general outline for the algorithm is shown in figure 6. The level set field has been discretized via a finite difference method on a fixed N x × N y × N z Cartesian grid of cubic cells.
The specimen is currently modelled as disconnected and freely floating in space. On each iteration of the model, the specimen surface mesh is extracted from the level set grid via the Marching Cubes algorithm [29] and isotropically remeshed using the C++ library CGAL [30] to provide a well-conditioned mesh of triangular elements. The isotropic remeshing algorithm aims to produce triangular mesh elements with panel edges of the same length [31], although this condition is not strictly enforced. In this study, the target edge length for panels has been increased every ten level set grid cells away from the apex, up to eight times the initial panel edge length, to produce an adaptive mesh. This initial panel edge length has been set to the average edge size of the original mesh extracted via marching cubes. The electric field is then solved over this specimen surface mesh via the boundary element method. The local electric field strength and evaporation field on the specimen surface determines the evaporation rate. The level set field is iterated by integrating equation (5) through Euler's method and a firstorder upwind scheme found in [32]. This update is performed using an adaptive timestep dt, determined by the Courant-Friedrichs-Lewy (CFL) condition given in equation (22).
where h is the level set grid cell width and δ is the CFL fraction. Throughout this study a CFL fraction of δ = 0.2 has been used. The codebase has largely been written in Python 3, with program bottlenecks implemented in C++ and parallelised. The parallelised components include the BEM matrix construction, BEM field solution, extension velocity construction, and ion trajectory integration, parallelised by singlehost threading. The time consuming bottlenecks include the BEM system solution, e.g. solving equation 4, and the velocity extension construction. The BEM system has been solved via the GMRES solver in Scipy. Ion trajectories have been integrated up to a virtual detector using an adaptive-step fifth-order Runge-Kutta (RK5) integrator [33]. Throughout this study, a simulated flightpath of L = 0.12 m and detector radius of R d = 0.04 m was used, matching the instrument setup of the LAWATAP instrument. Simple geometries such as multilayer and spherical precipitate phases can be constructed analytically. However more realistic and general specimen geometries can be initialised from 3D imaging data, whether created by the user or generated directly from experiment.

Results and discussion
In the following section, various case studies regarding model accuracy and applications are considered. The first case study (section 4.1) investigates the numerical stability of our 3D model through a sphere collapse test. The next two cases (Sec- how the model can simulate realistic emitter geometries, and how these models can be used to drive reconstruction.

Numerical stability
Too low a reinitialisation rate and the level set field begins to shift from its original signed distance form, increasing numerical error due to an uneven spacing between level sets [20]. However, too high an initialisation rate also results in the numerical drift of the interface.
In our previous 2D implementation reinitialisation was performed by rebuilding the level set field using a closest-point method. However, this approach was found to be relatively costly in 3D, taking 90.5 s per reinitialisation for the simulation to be described in section 4.5, or 95.6 µs per voxel. Instead, reinitialisation is performed iteratively via the relaxation scheme first proposed in Sussman et al [34]. This method was found to take only 1.1 s for the same simulation, or 1.1 µs per voxel. For a particular instance of the simulation at simulation time t, the level set field ϕ t (x) = ϕ(x, t) is relaxed according to equation (23): where ζ~h is a numerical perturbation to ensure convergence, h is the level set grid width, ϕ t (x, 0) is the level set field prior to relaxation, sgn is the sign function, and τ is a pseudo-time used to parameterise ϕ during the relaxation procedure. Like with equation (5), this system was solved via the Euler method and a first-order upwind scheme. In order to investigate the numerical accuracy of the implementation, a collapse test for a sphere held at a constant voltage u c was undertaken, similar to our previous studies [16,17]. Given that for a conducting sphere ||E(s)|| = uc R , where R is the sphere radius, the Arrhenius law defining the solution to the sphere collapse problem is given in equation (25): where for the subsequent collapse test, A = 1 m s −1 , β = 10, u c = 1 V, F 0 = 10 Vm −1 , and an initial radius of R 0 = 30 m was used. Due to the lack of an analytical solution, this has been solved numerically. The results of the collapse test using a reinitialisation rate every 30 model iterations is shown in figure 7(a). The test terminated once the sphere volume had reached 5% of its initial volume. It was found that a reinitialisation rate every 30 iterations was sufficient to minimise this drift and this rate was used throughout the study. However, even without reinitialisation the drift in the specimen surface is minimal, and is almost certainly negligible compared to the error introduced due to  physical assumptions e.g. equation (1) and equation (2). The error in surface panel fluxes for the surface field solution, shown in figure 7(b), likely arises due a combination of the surface disretisation and the relatively poor conditioning of the implemented BEM. Like with our previous work [17], it was found that using higher-order time integrators will provide little benefit until the BEM conditioning and quality of the field solution close to the boundary is improved.

Grain boundary simulation
In order to drive new distortion correcting algorithms in APT, our model must be capable of capturing the wide range of physical phenomena arising within APT experiments including specimen crystallography. Specific material grains can be simulated by defining anisotropic evaporation fields. The equilibrium emitter shape for a homogeneous specimen with a given anisotropy function ( figure 8(a)) is shown in figure 8(b). These functions can be locally assigned at different spatial locations within the model, with an examplar 15 o tilt grain boundary shown in figure 8(c). To simulate the grain boundary a level set grid size of 52 × 52 × 158 has been used. Here both Grain 1 (G1) and Grain 2 (G2) have the same evaporation function form (given by figure 8(a)), but the G2 function is rotated 15 o clockwise around the y-axis. The continuum model successfully reproduces both low density poles and zone lines, with the projection of the grain boundary appearing as an additional zone line-like structure down the centre of the simulated detector hitmap ( figure 8(d)).
Although the model can successfully reproduce pole structure, additional work is required before such a model can accurately capture the crystallographies of realistic APT specimens. Exact shapes for the anisotropic evaporation field functions for real crystallographic structures remain unknown and would have to be inputted as pre-determined model parameters. Controlled experimentation, density functional theory (DFT), or molecular dynamic models may assist in providing this data, although such methods are not explored here. However, this work demonstrates that continuum models can reproduce large scale crystallographic features in APT, raising the possibility for model-driven reconstruction to at least partially recover global crystallographic information within reconstructions. However, due to the continuum nature of the model, atomic scale effects such as roll up cannot be captured.

Void simulation
By implicitly tracking the specimen surface using a level set method, the merging separation of the specimen surface with void surfaces can be handled without having to explicitly join or separate meshes. Voids can be incorporated into the model as internal contours with a zero evaporation rate.
A proof of concept simulation for a void within a homogenous specimen (F 0 = 30 V nm −1 ) is shown in figure 5(a). The specimen geometry was constructed with a 25 nm tip radius, shank angle of α = 7.1 o , and a 10 nm void radius. The model containing the void had an initial shank length of 300 nm, with simulated evaporation on a 84 × 84 × 97 level set grid taking 9 minutes to run (see section 8 for system specification details). A calibrated point-projection reconstruction in constant shank mode was performed on APT data simulated using this model (figure 9(a)), with an ICF = 1.3 determined via correlation with the initial specimen geometry. Through spatial voxelisation, a density function was calculated for the reconstructed data ( figure 9(b)). This revealed a region of increased density due to the void, surrounded by a region of lower density. These density fluctuations arise due to trajectory focusing and largely agree with what has been observed in experimental studies [35][36][37] and simulation results from atomic models (see figure 8 in the supplementary material of Wang et al [37]). Trajectory calculations also revealed a spatial crossover region in the reconstruction, with ions originating from the void surface mixing with those from the bulk region.

Simulated reconstruction of an experimentally derived specimen geometry
So far all geometries considered have been constructed via simple analytic functions. However, geometries in experiment are rarely this perfect. In order to simulate and drive the reconstruction of real specimen geometries, electron tomography (ET) [38] can be used to provide the input for the model specimen geometry. This includes both the specimen surface and internal material phase boundaries. In order to investigate the accuracy and numerical stability of a model-driven reconstruction protocol, the same model is used to both generate synthetic APT data and then drive its reconstruction.
The dataset considered here is of a needle shaped Al−Mg−Si alloy containing MgSi precipitates, initially MgSi 2 that have undergone some dealloying process. HAADF-STEM imaging of the specimen was taken at 5 o tilt intervals, and reconstruction of the tilt series was then performed via the SIRT algorithm in the ASTRA toolbox [39]. The resulting tomogram is given by a 3D scalar field of contrast values. Thresholding is used to extract each phase, and edge artifacts removed by taking only the largest contour defining the specimen surface.
As effective evaporation fields for these ternary phases are unknown, substitute field values have been used, with precipitates assumed to be low field (F 0 = 13 V nm −1 ) compared to the bulk (F 0 = 17 V nm −1 ). While these substitute values are very likely unphysical (given that the evaporation field for pure Al is 19 V nm −1 and pure Si is 33 V nm −1 ), they serve in demonstrating how ET can be used to setup a realistic specimen geometry in the model. However, the results from this specific simulation should not be directly compared with experiment.
The result of simulated evaporation for this model is shown in figure 10. The model had an initial shank length of 590 nm, with simulated evaporation on a 39 × 42 × 297 level set grid taking 30 minutes to run (see section 8 for system specification details). The spreading of detector isolines (shown in red) surrounding precipitates captures the local trajectory focusing occurring in these regions. As each detector isoline corresponds to a constant position on the detector, a spreading of these lines in a region of the specimen implies local demagnification (as more specimen volume now corresponds to the same volume of detector space).
APT data was also simulated using the model and reconstructed via a calibrated constant shank mode point-projection reconstruction. An initial cap radius of R 0 = 22.5 nm and shank angle of α = 1.9 o were determined via a least squares fit to the initial model geometry. By manually sweeping through image compression factor (ICF) values, an ICF of 1.15 was identified as giving a good correlation between the initial precipitate microstructure and the reconstructed phases (figure 11(a)). The reconstruction reveals a compression of the low-field precipitate phases within the APT reconstruction typical of low-field microstructure.
To perform model-driven reconstruction, the model must be calibrated to experiment. This requires each experiment time (detected event index) to be matched to a simulation time. To obtain this mapping, landmarks are placed at specific experiment and simulation time points known to coincide e.g. microstructural landmarks. Between these landmarks, experiment and simulation times are related by the volume of evaporated material within the FOV. In this work, simulation frames are matched to microstructural landmarks (shown in blue) within the temporal histogram over experimental detection events ( figure 11(c)). An appropriate ionic volume for reconstruction between landmarks is calculated to ensure that for a particular landmark, the ion exactly detected at this landmark's  experiment time is placed in the reconstruction using the model trajectories calculated at its corresponding simulation time. Using these calculated ionic volumes and considering the evaporated volume throughout the simulation, a piecewise linear monotonic mapping is defined matching simulation and experiment times. These landmarks effectively constrain the depth positions of microstructure within the reconstruction to match those present in the model. Relaxing the constant ionic volume constraint ensures the model does not drift from experiment, while any discrepency in the calculated ionic volumes provides an additional metric for assessing the quality of model calibration and APT data ranging.
The resulting model-driven reconstruction is shown in figure 11(b), with precipitate ions shown in purple and bulk ions in green. Even from visual comparison with the pointprojection reconstruction ( figure 11(a)), it is clear that the model-driven reconstruction provides a significant improvement in the correlation of its chemical mapping with the original positions of microstructure. However, numerical errors are present in the calibration process. This can be seen by the marginal offset and distortion of the first encountered precipitate, as well as compression in the final precipitate: possibly due to breakdown of the trajectory mapping arising from trajectory crossover. Note that due to initial ion launch positions being known, direct quantification of ion placement accuracy is possible although omitted from this current study.

Experimental finfet reconstruction
Reconstruction distortions are particularly problematic when analysing semiconductor devices. Such systems typically consist of phases of varying evaporation fields and dielectric properties which can result in distortions within point-projection reconstructions that make key measurements of limited use ( figure 12(c)). It is worth noting that some of these distortions are so extreme that trajectory overlaps can occur in the central channel of the finFET [6]. These distortions cannot be corrected for.  We extend our prior 2D model for the finFET geometry into 3D, examining the same structure as in [6]. From SEM micrographs following specimen preparation ( figure 12(b)), an initial needle geometry with shank angle of α = 15 o , and initial apex radius of R 0 = 20 nm was estimated. As previously discussed, there is a lack of information regarding evaporation fields of compound phases. However, a previous study of the same dataset yielded an estimate for the evaporation field ratio of the SiO 2 phase to the central SiGe phase of 1.8 [6]. The evaporation fields obtained using the image hump model for Si (F 0 = 33 V nm −1 ) and Ge (F 0 = 29 V nm −1 ) are relatively close [40]. Therefore, for the purposes of this study the phase Si 0.25 Ge 0.75 is assigned an F 0 = 30 V nm −1 , giving SiO 2 a value of F 0 = 54 V nm −1 . The authors recognise that these field values are likely inaccurate. However, as only relative field differences between phases are important, any calibration errors should be confined to the region surrounding the Ge capping layer at the reconstruction apex. The model for the finFET had an initial shank length of 440 nm, with the simulated evaporation on a 84 × 84 × 134 level set grid taking 24 minutes to run (see section 8 for system specification details). In order to aid in the model calibration with the experimental data, maps showing the local phase over the detector at different experiment times were generated via k-means [41] (figure 12(e)). These experimental phase maps have then been compared with phase maps generated using the simulated data ( figure 12(f)). Comparing such maps after the emitter had formed a steady state allowed a value for the orientation of the SiGe central channel to be determined for the model geometry (estimating a 3.5 o rotation around the analysis axis).
As per the reconstruction in section 4.4, experimental landmarks (marked in blue) have been placed at specific times within the experiment detection event sequence. The experimental landmarks used for the finFET reconstruction in this study are marked on the histogram of the detection event sequence ( figure 13(a)). Matching simulation markers correlating to particular evaporation events within the ion sequence are shown in figures 13(b)-(d). With no features available to constrain the model start and end time points, simulation markers for landmarks A and D were chosen by manually selecting time points providing ionic volumes consistent with the rest of the reconstruction. Compared with the point-projection reconstruction, there is a clear improvement in the correlation of microstructure morphology between the model-driven reconstruction and electron microscopy data. The visualisation of the local density also indicates some significant improvements compared to the point-projection approach (figures 12(c) and (d)), especially within the central fin region. However, high density artifacts remain close to the interfaces due to imperfect calibration between the model and experiment. This calibration error can be directly observed by the offset in the position of the central fin interfaces between the generated reconstruction and model ( figure 14(e)).
The quality of the model-driven reconstruction can be gauged by comparing the density profiles of the model-driven reconstruction and the conventional point-projection reconstruction. To generate these density profiles, the ions within both reconstructions have been binned on a 100 × 100 × 100 grid. Histograms of the ion counts for all ranged species were generated for both reconstructions ( figure 14(f)). Edge voxels have been filtered out by discounting any voxel adjacent to a voxel containing zero counts. The tail of the point-projection reconstruction density profile (red) extends out to high densities away from the main peak centred around 5 counts per bin. In the case of the point-projection reconstruction a second smaller peak is also present, corresponding to voxels located within the central fin region. This spread in densities arises due to increased voxel densities within the central fin and is highly unphysical. In comparison, the tail of the model-driven reconstruction density profile tapers off significantly quicker, with no second peak being identifiable. The absence of the second peak supports what figure 14(c) implies, that the density distortions in the central fin region have mostly been corrected for.

Computational performance
Unlike our previous study, many of the major algorithmic components of the 3D model have been parallelised. This parallelisation extends to the BEM matrix construction and solution, extension velocity construction, and the ion projection. The only major step currently not parallelised is the isotropic remeshing -which is performed by the CGAL library [30] and currently does not support this functionality. All benchmarking has been performed on a standard laptop (i7-7600U CPU 2.8GHz system with 2 real cores). The times per iteration for the major algorithmic components are given in figure 15, obtained by averaging 10 iterations of the finFET model, as given in figure 12. This performance analysis reveals that for large system sizes, i.e. the asymptotic limit, the extension velocity construction, the ion projection, and field solution, currently calculated via the Scipy GMRES solver, become the algorithm bottlenecks. Typically ion projection is not performed every iteration of the model, and so in practice it is usually not the limiting step. Future work should focus on algorithmic optimisation of these steps.
Algorithmically, only the panel clustering and the isotropic remeshing step have been introduced over our previous 2D implementation [17]. Panel clustering has provided a significant speedup to the ion projection step. For example, measuring using the Python's time.time() function, the projection time for 1000 ions off the finFET geometry given in figure 12(h) has been reduced from around 60 minutes to only a few seconds. The result is the total time for ion projection, with one ion per surface panel, being reduced from hours to only a few minutes. Similarly, the isotropic remeshing step produces an adaptive remeshing of the specimen surface. The area of boundary elements increases with distance from the apex, reducing the total number of boundary elements present and accelerating all other algorithm components listed in figure 15. This isotropic remeshing also serves to improve the conditioning of the linear system in equation (9) by removing poor quality elements (see section 2 in the supplementary material). It is worth noting that the isotropic remeshing step is the only major step that is currently not parallelised.
A rough comparison to atomic models such as TAPsim [9] can be obtained by estimating and comparing the evaporation rate. The total evaporated volume of material in the finFET simulation over the 24 minutes runtime was 1.03 × 10 −20 m 3 . Assuming an ionic volume of 6.4 × 10 −29 m 3 (corresponding to a typical lattice spacing of 0.4 nm [40]), a real time evaporation rate for this specific simulation can be determined to be around 100, 000 atoms per second. In general this rate will depend heavily on the simulation parameters used.

Limitations and improvements
This study improves our previous 2D continuum modelling efforts by successful extension of the conductor approximation to 3D. Our previous work [17] demonstrates how this simulation framework could be extended to include dielectric materials. The potential of this model to drive distortion correcting reconstruction algorithms for general specimen geometries has also been clearly demonstrated, improving both microstructure morphology and local density. However, spatial calibration issues between the model and experimental data remain. Remaining density distortions are observable at the interfaces of the model reconstruction in figure 14(c), and while some of this could be due to trajectory overlap [6], the authors suspect a large part of these residual density distortions are due to model miscalibration. At the interface of the finFET there is a high density region in the detector data caused by trajectory focusing and demagnification (see figure  4 in Melkonyan et al [6]). Interface misalignment between the model and the experimental dataset during the model-driven reconstruction, as implied for the finFET by the overlay of the model and reconstruction in figure 14(e), would result in a lower than predicted magnification being used for the ions where demagnification should be at its maximum, causing an increased density within the reconstruction.
These interface calibration issues could potentially be resolved by a slight transformation of the experiment detector data to enforce spatial overlap between the high density interface focusing regions of both the experiment and model. The magnitude of the required transformation would provide an additional metric to assess the suitability of the emitter model used. One class of transformation methods we are currently considering for addressing this registration problem are diffeomorphic mapping-based methods [42]. These methods are commonly used in the field of computational anatomy for matching up experimental data with geometric models. Such transformations could potentially be minimised or avoided through shape optimisation of material phase boundaries within the model. However, shape optimisation would require significant acceleration of the model.
The current model implementation does not include a local electrode or instrument chamber, which is known to effect the electric field solution and resulting ion trajectories [43]. The truncation of the shank will also introduce some error into the trajectory compression [43]. As using a boundary element method (BEM) removes many of the computational issues arising due to differences of scale, future implementations could directly incorporate these components into the electrostatic solution although further optimisation of the electrostatic solver will likely be required. Another possible solution could be to apply some kind of trajectory transfer function, perhaps a projection law, that accounts for this additional compression. However, the authors are currently unsure about the validity of such a method to specimen geometries with concavities.
Currently panel clustering is only being applied to the domain electric field calculation step. Future models could further accelerate model electrostatics by solving equation (9) via panel clustering or a fast multipole method [25,28]. Ion projection could also be further accelerated by a fast multipole method, using the method to evaluate electric field points close to the emitter.
The requirement for isotropic remeshing is likely to be a simulation bottleneck due to its current lack of parallelisation. Further optimisation or avoidance of the isotropic remeshing step will be required for increased automation of the reconstruction procedure. Possible solutions include an adaptive grid level set algorithm [43], alternative isosurface extraction algorithms [44], or by tracking only changes in the specimen apex using the level set method.
Future implementations could automatically constrain model parameters directly from the experimental data e.g. by maximising overlap between simulated and experimental phase maps (section 4.5) or correlation of the reconstruction with the model geometry or electron tomography data [45,46]. The best data features to use for calibrating models with the experimental data remain unknown. Combined APT-ET datasets will be critical to the further development and testing of this algorithm beyond multilayer geometries.
Finally, it is worth noting that the laser has so far been omitted from our models. Experiment has shown that the laser has a major effect on the specimen geometry [47] Therefore incorporating its effect into the model is critical in order to obtain accurate reconstructions and remove distortions. While the level set method can handle arbitrary evaporation rate functions, e.g. equation (1), the exact physics behind laser-assisted field evaporation, and therefore the correct evaporation function for the laser, remains unknown. One option is to define a spatially anisotropic temperature distribution over the specimen surface [17,48]. Hatzoglou et al estimated such a distribution through simulation [48]. However, Vella et al have suggested the exact effect of the laser is likely more complicated [49].

Conclusions
In this study we presented a full 3D continuum model for specimen field evaporation in APT. By tracking the specimen surface using a level set method, our model can capture a wide range of APT phenomena including crystallographic faceting and evaporation of voids. At practical simulation resolutions (level set grid cell widths of 1-4 nm), times on a standard laptop or desktop typically take tens of minutes compared to the hours or days required by atomic resolution electrostatic models.
This work demonstrates the key components required for a model-driven reconstruction algorithm. APT data was simulated using the model for a realistic specimen geometry initialised from ET data. This data was subsequently reconstructed using the same model in a round trip test, confirming that numerical errors are sufficiently minor to provide an improvement over a calibrated constant shank mode point-projection reconstruction. Finally, a model-driven reconstruction for an experimental finFET dataset was carried out using a manually calibrated level set model. Once again the model-driven reconstruction showed a quantifiable reduction of unphysical high density distortions compared to the constant shank pointprojection reconstruction. Given combined correlative ET and APT datasets, and further improvements to interface calibration, such a method could be used to drive distortion corrected reconstruction for very general and complex specimen geometries.

Acknowledgment
Charles Fletcher acknowledges support from CAMECA ® for financially supporting this research. Michael P Moody and Daniel Haley acknowledge support from the ESPRC project Advanced Nuclear Materials EP/P001645/1. We would like to thank Thomas Slater (Diamond Light Source) for providing the Al−Mg−Si alloy ET data, and Claudia Fleischmann (IMEC) for providing the finFET APT dataset and accompanying electron micrographs. Charles Fletcher would also like to thank Brian Geiser and Joe Bunton for useful discussions regarding the problem.

Source code
The source code for the 3D simulation tool can be accessed and downloaded from https://gitlab.com/fletchie/casra-3d. All simulations were performed on Ubuntu 18.04 using system Python 3.6, with Python library versions installed via apt and linked against OpenBLAS version 0.2.20 and LAPACK version 3.7.1 supporting multithreading. All library versions used are fully specified by the Ubuntu distribution and can be checked at packages.ubuntu.com. All simulation benchmarking was performed on a standard laptop (Dell Latitude 7480, i7-7600U CPU 2.8GHz system with 2 real cores)