Mesoscopic modeling of field evaporation on atom probe tomography

Reconstructions in atom probe tomography (APT) are biased by image distortions arising from dynamic changes of the specimen geometry that controls image projection. Despite the strong efforts to build realistic models for understanding and reproducing image artifacts, the current models are too slow or not adapted to be routinely used in image correction approaches. To understand the APT imaging process for real size samples submitted to realistic experimental conditions of electric field and temperature, we propose an alternative simulation tool based on a coarse-grained model of the sample surface. The surface electric field on a meshed surface is calculated by using continuous models describing field evaporation. The dynamic evolution of the sample surface and the image projection are predicted using materials properties. We show that the interplay between temperature and electric field is an important ingredient in predicting the ion projection, in pure metals and in more complex materials. This fast approach accurately reproduces the well-known local magnification and trajectory overlaps effects in the evaporation of small particles. By combining prior knowledge about the sample structure and properties, the model could be used to improve the reconstruction approaches for complex sample geometries.


Introduction
Atom probe tomography (APT) is a powerful technique to analyze materials in three dimensions (3D) at the atomic scale * Author to whom any correspondence should be addressed.
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. [1][2][3]. The main physical principle is the field evaporation (FE) of atoms which are located at the material surface. FE is the transition, in the presence of a sufficiently strong electric field, of a material directly from the solid phase to the plasma phase without passing through sublimation. For most materials, the critical field giving rise to FE without any thermal activation is in the range 10-60 V nm −1 . This intense electric field is easily achieved at the surface of a material using the classical electrostatic point effect. The material is prepared as a sharp needle with an end radius R in the range of 10-100 nm and is submitted to a high voltage V (a few kV) to achieve the required electric field. When FE is carefully controlled (by tuning V), atoms are liberated one by one from the surface, ionized, and projected in the diverging electrostatic field existing close to the tip surface. The sample acts itself as an electrostatic lens, projecting ions in a deterministic way. By collecting ions on a position sensitive detector placed in front of the sample, the initial positions of all atoms may be in turn be retrieved. In a single-phase material, FE is smooth and regular, and the ion projection is straightforward. APT tomographic reconstructions of these materials achieve excellent metrological performances [4]. However, in materials composed of several phases of different chemical composition, image reconstruction is often plagued by severe distortions.
The major cause of distortion is the complex, dynamic evolution of samples by FE. This complexity is linked to the interplay between the structural and chemical composition of the materials, and the intensity of the electric field on the sample surface that is governed by the 3D geometry of the sample surface, i.e. the atomic arrangement. We may note that the temperature of the sample surface can also be non-homogenous giving rise to additional distortions. Appropriate reverse projection algorithms must be developed with a correct understanding of the dynamic evolution of the FE at the sample surface. As a result, the evaporation process has to be understood, in order to become fully quantify and thus to enable a possibility to systematic accounting for and correcting of as many of the artifacts inherent in the current reconstruction protocols as possible [5].
A common approach is to use numerical models, which have proven to be useful to understand some known issues. These simple approaches using analytical sample shape are clearly insufficient [6]. Atomic scale approaches to FE were developed in the early 2000s [7,8]. These kinds of models were sufficient to give an overview of the experimental shortcomings, such as improving the understanding and providing correction methods for trajectory overlaps during evaporation of small precipitates. For instance, in [8], the artifacts that occur during FE of nano-oxides embedded in a metallic matrix were quantified and corrected with this model. More recently, the artifacts observed in the APT analysis of nanometric voids induced by irradiation in metallic materials were successfully explored. Simulated APT experiments demonstrate the local density variations near voids are controlled by the unique ring structures as the volume of voids gets exposed and the different evaporation fields of the surrounding atoms [9].
Atomic scale models for FE model have shown their value for understanding imaging artifacts. The basic principle is the description of the emitter geometry at the atomic scale via an adaptively, refined mesh. This discretization is used to solve the electrostatic Laplace equation to obtain the 3D distribution of the voltage and the electric field. The atom exhibiting the highest field is selected to be field evaporated [10]. The mesh is also useful to compute trajectories of desorbed ions towards a virtual detector. Since then, a lot of improvements have been introduced by different groups in the APT community to model emitter structures quite comparable with real samples [11][12][13]. Main modifications of FE models were mesh refinements which enabled computations on more realistic specimen shape structure and external boundary conditions for the electrostatic field and chamber [14]. Recently, a new model was developed by Rolland et al [15] disregarding the mesh structure. This model is based on a direct evaluation of the charge distribution over the conducting emitter surface as a solution of the electrostatic Robin's equation [16,17]. In practice, this is equivalent to solving an N-body problem (where N is the number of atoms at the emitter surface). An approximation for handling long-range interactions was introduced by Rolland et al [15], developed by Barnes and Hut [18], in order to reduce sequentially the time complexity for solving the Robin's equation from O N 2 to O (Nlog (N))).
Recently, FE simulations are no longer only used to understand the evaporation process, but also used to perform tomographic reconstruction [19].
The ultimate dream of using FE models to calculate reverse projection, atom by atom was proposed by several authors [19,20]. However, atomic models are still far away from giving perfect knowledge of the FE process. For instance, zone line artifacts, linked partly to rolling-up of atoms in the first steps of flight are not fully treated by current models [21]. The problem is that atomistic models, so far developed [14,15], do not take into account some complex atomic scale phenomena such as molecular dissociation [22], field induced surface diffusion [23], or rolling-up motion [21]. To achieve a better description, modeling at a sub-atomic scale is required, which can be done either by density functional theory (DFT [24,25]), or by molecular dynamics simulations [26][27][28][29]. Attempting to incorporate sub-atomic physics in current models drastically increases computational complexity, resulting in a reduction of model sizes that can be physically simulated. On the other hand, some necessary approximations are also used at the macroscopic scale, so that atom probe chamber details and limitations of the position-sensitive detector limitations are often ignored in models. We note several limitations in APT experimental conditions, such as detection efficiency and detector spatial precision (detection efficiency between 30%-70%, multi-hit capability, field-of-view ±20-30 • , spatial resolution ∼0.2 mm). As a result, the exact comparison of atomic scale FE model and experimental data is clearly not feasible.
However, it does not mean that a simulation tool cannot be used to improve the reconstruction algorithm significantly. Most of the image projection is governed by the sample shape, thus an accurate understanding of this shape would enable reconstructions with less distortion. The value of incorporating atomic details is unclear, since a lot of subtle physical artifacts may disturb ion trajectories in a non-deterministic way at this scale. For this reason, but also for the sake of saving computing time, FE models at the mesoscopic scale have been recently developed. Mesoscopic models assume an evolving sample surface with a defined curvature (tensor) field whose evolution is simulated as a function of time. Numerically analyzable in seconds or minutes on a laptop compared to days and months on cluster computers makes mesoscopic models useful tools for research on new methods of reconstruction [19].
Mesoscopic FE models were first developed by Haley and co-workers using level set approaches. First attempts using 2D approximations were found promising [30], and 3D extrapolation is in progress [31][32][33]. We propose in this paper an alternative approach to simulate FE at a mesoscopic scale (nm to µm scale). The model discretizes the evolving sample surface as a triangulated surface mesh. The evaporation sequence is simulated by estimating the evaporation rate of all the mesh primitives (vertices and facets of the triangles supporting the mesh). Note that the physical process of the FE is injected using a continuous model of FE, that is parameterized to handle the local nature of the material (FE constants, as a function of the specimen structure and composition), but also the local temperatures. Results of this modeling are presented for microstructures commonly encountered performing APT analysis. Results from this simulation are compared to experimental ones, demonstrating the validity and the relevance of this FE model at the mesoscopic scale.

Three-dimensional geometry of the field emitter
The first step of the modeling process is to define the shape and the dimensions of the field emitter in 3D. As mentioned in the introduction, an APT sample is prepared as a thin needle with an initial radius of curvature at the apex in the range 10-100 nm. This allows, through application of an electrical potential (in the range of kV), the generation of a sufficiently intense electric field (few tens of V nm −1 ) to ionize and then to evaporate the atoms located at the apex surface. Whereas the definition of the sample was previously done with the coordinates of each atom of the emitter [13,15], the emitter surface in this approach is geometrically triangulated. The mesh size defines the minimum scale of features to be modeled as compared to an atomically refined approach. With this model, the sample information (and computing) is reduced to vertices of triangles of the surface mesh which discretizes the field emitter surface. For comparison, atomistic models need to store the information of atomic positions of the entire volume, to finally extract only surface atom positions [15]. The mesoscopic approach presents a true advantage in terms of memory allocated for calculations. At the atomic scale, the allocated memory is in the range of GB (at least) [8,13,15] while here it is, at best, one hundred times smaller . For the initial emitter surface (i.e. before evaporation), the coordinates of the vertices are defined, such as after meshing the surface, with a triangulation, each mesh has the same and the expected area (S at ) [34]. Figure 1 represents such an emitter with an initial radius of curvature (R c ) of 20 nm, a shank angle (γ) of 10 • , a length (L) of 100 nm and a mesh element area (S at ) of 4 N m 2 per triangle.
To characterize a surface and its evolution during evaporation (presented in the next section), the radius of curvature is commonly used. Indeed, the correlation between the electric field F and the radius of curvature R c is often considered using the standard expression [35]: with k f the field factor (linked to the specimen morphology [4,6]) and V the applied potential to the emitter. The radius of curvature is the reciprocal of the curvature, which is determined here using the approach of S. Rusinkiewicz [36]. Curvature quantification is done generally by calculating the mean curvature (average of the principal curvatures) or the Gaussian curvature (product of the principal curvatures). In this study, the mean curvature of each mesh is computed (as done by [36]) to estimate its associated radius of curvature (figure 2(a)). The values thus obtained (figure 2(b)) agree with the simulation input parameters (figure 1). Indeed, the radius of curvature distribution is a bimodal distribution with one mode centered to the set apex radius of curvature (i.e. 20 nm in figure 1) and the other mode with higher radius of curvature (>40 nm) corresponding to the shank (due to the shank angle, this last curvature radius continuously increases as a function of distance to the tip apex).

FE
The FE process of this simulated system is based on the local FE rate at the surface of a material submitted to an electric field F. Assuming that the atoms evaporate from a surface at a defined rate, the surface center will move downward along its normal direction. From the evaporation rate equation, assigned to each primitive of the surface mesh (figure 1), the volume variation can be estimated. Using geometrical considerations, it is then possible to estimate the displacement vector associated with each vertex. A new surface is thus obtained. This operation is iterated until a specified end condition is reached (commonly an a priori defined evaporation depth). Experimentally, the evaporation process is triggered using pulses (electric or laser). Regardless of the pulsing mode, the evaporation rate (K) for a single atom is described empirically by an Arrhenius-type law. Extrapolating to our triangulated system, the evaporation rate of a mesh j (K j ) is then expressed as: with the activation energy barrier Q (since this process is usually considered thermally assisted), Q (F j ) is the energy barrier considering an external electric field (F j ), k B the Boltzmann constant and T j is the temperature. A linear approximation can be used to describe the dependence of the energy barrier and the electric field [37]: with C j as an energetic constant and F EV,j the evaporation field. The evaporation field is the value of the field at which the energy barrier is reduced to zero. Note that more refined equations defining Q (F) can be injected in the model without major modifications of the algorithm [1,38].
In the case of pulsed evaporation, the expressions (2) and (3) are time dependent, since both F or T may be time varying. If the evaporation process is triggered by electrical pulses, then the base temperature of the sample can be considered as a constant. Therefore, the temperature of each mesh (T j ) is set as a constant, defined by the user according to its APT analysis conditions (figure 4(a) with a constant temperature equal to 50 K). Most of the evaporation occurs at the top of electrical pulses. We can then model the evaporation rate per pulse by: With τ F being a fraction of the voltage pulse applied to the sample [1], and F j being the peak electric field induced by However, if the analysis is triggered by laser pulses, the temperature can no longer be considered as a constant [39,40]. Indeed, the interaction of the laser with a nanometric needle induces local enhancement of the temperature, which in turn triggers different local probabilities for the FE. During the pulse, each triangle is submitted to its own temporal evolution of temperature, with a temperature maximum T max reached in less than a nanosecond in the general case, followed by a temperature decay that can vary slightly on the surface. Assuming a simple shape of transient temperature pulse of width τ T (close to the decay time constant of a few ns), the atom flux ϕ T can be approximated by [41]: It is also reported for a mesh j that is a CMi, its evaporation rate (K j ) and its associated normal direction ( − → n j ).
Since the prefactor of the flux equation is linearly dependent of T max , or F, the exponential term is the main source of flux variation as a function of F or T. In a first approach, the use of a constant prefactor τ ′ The displacement vector of a defined vertex i ( − → D i ) results then from the evaporation rate of all the mesh (K j ) that are connected to it ( j ∈ CM i in figure 3), considering their outward normal direction ( − → n j ), since the electric field at the surface (F j ) is collinear to the latter: A constant parameter (ε) set by the user is added in the expression determining the vector displacement. Thus, it is considered that at each pulse (i.e. here each evaporation simulation iteration), the sample is evaporated by a certain amount. In order to be closest to the experiment, this quantity must be small (ε < 0.1 nm). This quantity is in fact the time step of the numerical scheme. Since the maximum evaporation rate is fixed, the displacement in depth is proportional to the time step between two iterations. The stability of the scheme, of first order in this simple approach, is strongly dependent on this quantity. Systematic tests were performed on different geometries to evaluate the propagation of errors of this approach for nanometer-size meshes with ε < 0.1 nm found experimentally as a good compromise between scheme stability and efficiency of the method.
Thus, at each iteration, it is necessary to determine four parameters for each triangle mesh j, in order to estimate the FE process of the sample, these are: T j , C j , F j and F EV,j . In fact, only three parameters are here considered, since the energetic constant (C j ) will be fixed as a constant for each mesh. This constant, as documented in [3,42] was only determined for very few materials, with values in the range 0.5-3 eV. It has to be set by the user in the range of tenth to few eV to correspond to experimental measurements [37,[43][44][45][46]].

Temperature and phase
In laser assisted atom probe, the maximum temperature is asymmetric on the tip apex, inducing then schematically an illuminated so-called laser side and a shadow side. The temperature distribution along the apex due to the laser interaction can be accurately reproduced with advanced models, such as proposed by Houard et al [47] (but also exhibiting in return a huge time computing). In this study, however, an analytical model is preferred, with a polynomial description of the distribution of peak temperature on the sample (expression (7)). This expression can reproduce the asymmetric distribution of the temperature (figure 4(b)), With X i and Y i the coordinates of the vertex i, α and β are parameters defined by the user to introduce different temperature evolution trend and T 0 the base temperature also the temperature in the shadow side of the sample. This last expression provides the temperature that is assigned to each vertex. The temperature that is associated to each triangle (T j ) is the average temperature of its vertices.
The different phases in a material are distinguished by their evaporation fields. In the case of a single-phase material, the evaporation field of all the mesh (F EV,j ) can be then considered as a constant (in the range of tens V nm −1 ). However, in multiphase materials, before each evaporation iteration, the evaporation field associated to each vertex needs to be defined. To evaluate these evaporation fields, at each evaporation step, the position of each vertex is compared to the position inside a 3D volume that defines the phases of interest. This volume contains the information of shapes, sizes, (particles, layers, …) and elemental nature (evaporation field) that are used as input data.

Electric field
It remains now the most important point, which is the calculation of the electric field for each mesh (F j ). In this study, the electric field is obtained from the estimation of the charge distribution on the emitter surface. This model, as is the one developed by Rolland et al [15], is based on the direct evaluation of the charge distribution over the emitter conducting surface, as a solution of the electrostatic Robin's equation [16,17]. As demonstrated by Rolland and collaborators, and extrapolated in this modeling approach, the charge q j (a positive non-integral partial charge) on a defined mesh j having an area (S at,j estimated at each evaporation iteration using the Heron formula) can be written as: With − → n j is the outward normal direction to the mesh j (as in expression (6)), − → r j,k the vector between the mass center of the mesh j and k, and N the number of meshes at the emitter surface. The mass center of each mesh is easily estimated from its vertex coordinates. It is then possible to construct a sequence q j,n for each mesh j to find the charge distribution over the field emitter surface: With q j,0 being an arbitrary positive number. It has been proven that this sequence tends toward the charge distribution at the equilibrium [48]. Moreover, this sequence has the advantage to converge very rapidly to the equilibrium (figure 5), even when starting from a random distribution of charge. Indeed, in the present case (figure 1), less than 15 iterations are required to reach the equilibrium (figure 5).
As expected, the highest charges are located at the apex of the emitter which is also in those regions where the electric field is the highest (electric field amounts to q j /S at,j ε 0 or σ/ε 0 with σ = q j /S at,j in accordance with basic electrostatic laws). We assumed that the equilibrium charge distribution is reached through equation (9), if the average relative deviation of the charge distribution between two iterations (⟨δq j,n+1 ⟩) is lower than 10 −5 % with It must be pointed out that the total charge carried by the conductor remains the same throughout the iterative process. To ensure the same behavior for the expression (9), q j,n /S at,j will be computed, at each iteration, such that:

Trajectory computation
Up to now, the model describes the evolution of the emitter morphology (and the electrical field which derives from the charge distribution), which is an important piece of data to understand APT analyses. By solving Newton's equation of motion for ions from vertices up to a planar detector situated a few centimeters in front of the sample for each step of the FE process, the dynamic and evolving projection APT law can be deduced. A numeric integration of this equation is conventionally done with an embedded fourth order Runge Kutta (RK4) method (at least for the last updated models [12,13,15]) or even a Leap Frog (LF) method for the oldest model [8]). For each method, there are pros and cons. In order to get the best compromise between speed and accuracy, the Verlet method was chosen in this model [49]. The Verlet method combines a good spatial accuracy (comparable to the RK4 method) and a low computation time (comparable to the LF method). Moreover, as it is conventionally done in FE modeling an adaptive time step calculation in solving the motion equation is introduced, since the emitter is at nanometric scale (≈100 nm) and the detector is located at metric scale from the emitter (≈0.1 m). Such adaptive time scale (dt) is determined using the velocity (V) deviation between two iterations: With t 0 the initial time calculation (⩽10 −16 for an accurate computation of the trajectory). The velocity deviation provides a more adaptive time step calculation than the velocity itself. Indeed, up to a certain distance from the emitter apex (few times the radius of curvature, figure 6(a)), the ion velocity increases drastically and then the velocity stabilizes (very low increase), while the velocity deviation (∆V) evolves continuously up to the detector. With the adaptive approach, it requires less than 100 iterations to calculate the trajectory up to a detector located at 10 cm from the emitter apex. In this way, most of the iterations are used to compute the trajectory close to the sample apex, where they are most prone to lead to trajectory estimation imprecision related to surface topography and the rest of the iterations are used for trajectory computation far away for the sample, where they are almost straight lines [1,2].
Experimentally, it is well known that due to the sample morphology and the analysis chamber setup (in particular an electrode, on which the electrical pulses are applied), the ion trajectories are not straight lines, but become compressed [50]. This compression is defined and quantified by a parameter, the so-called, image compression factor (ICF or ξ ). The parameter is expressed for each meshes j, considering small angle approximations: with L the flight length (distance between the sample and the detector), R c,j the radius of curvature of the mesh j and S D,j the mesh area on the detector. The experimental value of ξ is in the range of 1.3−1.6. Assuming an emitter of reduced dimensions (less than a micron in depth), the calculated ICF with the model is naturally reduced as compared to experimental values (blue histogram on the figure 6(b)). To achieve the exact experimental ICF value, it would be necessary to model the full geometry of the sample (needle on its coupon or capillary, local electrode atom probe (LEAP) sample holder, and local electrode ahead). We observed that increasing the length of the simulated sample, increases the average ICF, which corroborates this. Obviously, increasing the length of the tip, also increases the number of meshes and therefore the time computing. In order to maintain a reasonable computation time (sample length around hundreds of nanometers) and to have an ICF equivalent to the experimental one, we performed a post analytical compression of the trajectories, using an ad hoc coefficient λ. After a certain distance from the sample (in the range of few µm corresponding to the distance to an electrode and where the trajectories are straight lines), the trajectory calculation is stopped, and the position on the detector is obtained by a projection, using the velocity at this step, considering a compression of its polar angle, by θ = λθ ′ (with θ and θ ′ respectively the polar angle of the velocity where the calculation is stopped and on the detector), while its azimuthal angle remains unchanged φ = φ ′ (with φ and φ ′ respectively the azimuthal angle of the velocity at the same step as polar angle), as experimentally observed [1]. Using this approach with a coefficient λ equal to 1.5 (red histogram on the figure 6(b)), the average ICF (i.e. 1.53) is then comparable to experiment value. This coefficient is adjustable by the user in order to best adjust to the real analysis.

Algorithmic optimization
One of the crucial points of any simulation model for practical application are computing time and memory consumption. In this section, we show some solutions added to the model to address this. We are not benchmarking the performance of them, as done in [51], but just give a brief overview of the improvements, in term of computing time (table 1).
As it was mentioned previously, the electrical field is estimated by solving the Robin's equation with an O N 2 complexity (with N the number of meshes). Thus, the first introduced optimization, is a parallelized computation (table 1). Observe that in equation (9), step n+1 is uniquely dependent on fixed values at step n, with demand for an exchange of pieces of information for each point of charge, parallelization is strongly efficient with a scaling close to the number of cores (Using C# parallel loop). The second improvement, as Rolland et al [15] at atomic scale, is a long range interaction approximation, developed by Barnes and Hut [18], in order to drastically reduce the computing time required for solving Robin's equation solving by reducing the complexity to O (Nlog (N)). The time computing is directly linked to the number of meshes N, so we also introduced logically a reduction of this number (table 1). Regarding the electrical field distribution (figure 5), the evaporation will mainly occur at the sample apex. During the whole process of FE, the model consists in moving the meshes as a function of the evaporation rate (the latter of which is field dependent). Obviously, a refined mesh at the bottom of sample is unnecessary as the electric field on these regions is much lower than the electric field on the apex. An average charge on these meshes is also sufficient to describe the longrange interaction for ion trajectory calculations. As a result, we introduced an evolving mesh area that depends on the sample Table 1. Time computing (second) to perform ten field evaporation iterations (≈0.3-0.6 nm depth erosion in this case) as a function of the number of cores used (the computation time is estimated using a quad core processor). All the sample have a radius of curvature of 30 nm, a length of 100 nm and the meshes at the apex have an area of 2 N m 2 . The computation is done considering if the meshes area is constant (i.e. α = 0 and β = 0) or evolving according to the equation (7) (which defines the coefficients α and β).

Single core
All cores α = 0 and β = 0 400 139 α = 10 and β = 2 50 18 α = 50 and β = 2 12 5 depth, using a polynomial law (as for the temperature (7) but this time with a depth dependency) ( figure 7(a)). With this approach, a small mesh area at the apex is defined, so that an accurate estimation of the electrical field is achieved, where the evaporation process prevails, and a large mesh area is imposed on the shank, where the evaporation rate is very low, regarding the charge distribution ( figure 7(b)). This approach significantly reduces the number of meshes and therefore the computation load.
The input of all these optimizations makes it possible to simulate the evaporation over a depth of 100 nm on a realistic description of an APT sample keeping a fine scale description of the sample apex (with a mesh area at the apex of 2 N m 2 ) with computation time at the hour scale, on a standard computer. The same result on atomistic models takes days [8,15] or weeks [12,13].

Comparison with experimental results
The advantage of atomic models [8,12,15] is the possibility to directly compare the results with the experimental ones. To perform such comparison with our model, we must transform meshes into 'atoms'. Between two iterations, we associate to each mesh an evaporated volume (V ev ), which corresponds then to a number of evaporated atoms N ev (V ev = N ev V at /Q with V at the atomic volume and Q the detection efficiency). Then, we randomly place this number of atoms N ev in the regions where the meshes impact the detector. Some examples of this transformation will be shown in the next section. Figure 8 represents results of the simulation of FE of a singlephase emitter. The initial morphology (before evaporation) was defined with an initial apex radius of curvature of 20 nm, a length of 100 nm and a shank angle of 10 • . The mesh elements area at the apex were set to 4 N m 2 with the evolving mesh method (see section 2.6). From this model, 20 nm has been evaporated to test the simulation. More than considering only a constant temperature ( figure 8(a), T = 50 K), also nonhomogeneous temperature has been tested to check the ability to reproduce laser assisted FE ( figure 8(b)). A gradient of temperature was imposed along the X axis of the tip following a polynomial law (expression (7) with β = 2). The peak temperature on the side of the tip was set to 200 K (α = 3).

Single phase
At a constant temperature (T = 50, figure 8(a)), model predictions agree with those from classical experimental observations. During the evaporation process, the tip apex maintains its quasi-hemispherical shape, with an increase of the radius of curvature (initially 20 nm and at the end 27 ± 2 nm, figure 8(c)), due to a non-null shank angle as a function of depth erosion. The meshes, that impact the detector (R = 4 cm, at 11 cm from the tip, with λ equal to 1.6), are coming from a small spherical region from the apex (R = 10 nm) and a number of impacts distribute centrally around it.
The introduction of a gradient of temperature on the apex temperature, in order to reproduce the laser assisted evaporation, modifies the evaporation sequence ( figure 8(b)) and thus reproduces some experimental observations. The evaporation rate is obviously locally increased by the local high temperature (expression (2)), thus the tip shape gradually becomes asymmetric, as observed experimentally by [52]. The impact distribution on the detector is no more centered (highest density of impact is shifted to the direction of the laser side). By back projecting the detector onto the tip surface, it is interesting to note that the meshes impacting the detector are located on the tip on a deformed elliptical region shifted to the shadow side of the sample. This kind of image projection deformations are experimentally ignored. Considering the dimension of this last region (X ∈ [−15,5 nm]), the variation of radius of curvature on either side of the sample, which therefore impacts the detector, is equal to 5 nm (from 22.5 to 27.5 nm, figure 8(c)) and the electric field variation is equal to 5 V nm −1 (from 35 to 30 V nm −1 , figure 8(d)). In experiment, this surface electric field variation generates variations in the post ionization of the emitted ions. In turn, this results in an evolution of the charge state ratio from the laser to the shadow side [53][54][55]. It is important to note that this asymmetry (and the effects resulting from it) become even more pronounced when the temperature gradient increases (figures 8(c) and (d)).
We may note that the empirical relationship F ∼ V/R is not verified through this simple simulation which may be surprising. However, this empirical expression is generally given for a sample of cylindrical symmetry, which is not the case here. Our model could be then used to understand composition measurement fluctuations, induced by the laser assisted evaporation [56] or even to estimate the temperature rise at the sample apex induced by its interaction with the laser [57], introducing an accurate evolution [41,47].

Multiphase FE: planar layers
The first layer case that we studied with our model is a horizontal layer ( figure 9(a)), which is also the most studied case with FE simulations [30][31][32][58][59][60][61][62][63]. A key advantage of our model is its ability to model samples with dimensions comparable with experimental ones. Thus, we simulate FE for a sample with a radius of curvature of 50 nm, a length of 200 nm and a shank angle of 5 • , on which we introduce a layer with a thickness of 30 nm. The mesh elements area is set to 4 N m 2 at the apex, using the evolving mesh method along the shank. The simulations assume constant sample temperature of 50 K. Two evaporation field values for the layer (F Ev,B ) have been studied: 0.7 and 1.4 times the matrix one (F Ev,A ). We will focus on describing the first case (F Ev,B = 0.7 ×F Ev,A ) as the second case is similar, except for the observation that the evaporation sequence is reversed, as we will be discussed later.
Regarding the evolution of the field at the surface during the evaporation ( figure 9(b)), we define two regions of interest (ROIs): the interface A/B (figure 9(c), ROI 1) and the interface B/A ( figure 9(d), ROI 2). During the evaporation of the interface A/B, we observe an increase of the field on the A phase, due to its protruding position. This expands the magnification locally and decreases the field of view ( figure 10(a)) since only a limited region of the tip apex is collected on the finite size detector, (here detector radius is R = 4 cm, at 11 cm from the tip, with λ equal to 1.6). At this step, we also observe a slight increase of the electric field at the interface B/A, induced by a local decrease of the radius of curvature, because the evaporation rate of A phase is lower than the B phase. During the evaporation of interface B/A, the reverse phenomenon is observed. The B phase becomes flatter (high radius and low field), due to its higher evaporation rate compared to the phase A below,  inducing then a drop of the magnification and an increase of the field of view ( figure 10(a)). The radius of curvature of phase A decreases regarding its high evaporation field.
For the second case (i.e. F Ev,B = 1.4 ×F Ev,A ), the situation is the exact opposite. The evaporation of the interface A/B corresponds to the evaporation of the previous B/A, and reciprocally for the B/A interface. We obtain then an increase of the field of view at the A/B interface and the opposite at the B/A interface ( figure 10(b)).
Note that these results and observations agree with those which have already been obtained in the literature, but they have been obtained (i) on a real size tip and this (ii) with a relatively short computation time (2 h for each evaporation field).

Nanoparticles
The FE process of second-phase particles has been widely studied with atomic modeling [8,10,13,64,65], for almost 20 years now [7]. APT reconstruction particles are known in the literature to be strongly biased by artifacts named local magnification effects and trajectory overlaps [8,66,67]. Models enabled to obtain better quantitative understanding of these phenomena at least in idealized geometries. However, the cumbersome nature of atomic scale simulations often limited the result to a small range of particle size, or FE properties.
Local magnification occurs when different phases with individually different evaporation properties and fields are evaporated at the same time. Those phases (for example A and B) will develop locally different radii of curvature, (respectively R C,A and R C,B ), governed by their individual evaporation field and evaporation properties, which leads then to different local magnification (respectively M A and M B ). In a simple approximation, the following expression correlates their magnification, according to the local radius of curvature: This leads to trajectory aberrations. These are responsible for morphological distortions and in the worst case cause trajectory overlaps (when trajectories from different phases impact the detector at the same place), inducing then chemical biases. Those aberrations lead then to erroneous reconstructed atomic densities as was demonstrated analytically in [5,64].
We have simulated the FE of a spherical particle, with a radius of 4 nm that is embedded in a sample with a radius of curvature of 20 nm, a length of 100 nm and a shank angle of 10 • (figure 1). The meshes area is equal to 1 N m 2 at the apex, with an evolution along the shank. The simulations operate at a constant temperature equal to 50 K. Different evaporation fields of the particles (F Ev,B ) have been studied, from 0.4 to 1.4 times the matrix one (F Ev,A ) (by step of 0.1). After a transformation of meshes into atoms, associating an elemental nature corresponding to the phase from which they originate, a 3D reconstruction is performed, from their positions onto the detector, using standard reconstruction protocol [68,69]. Concentration profiles were computed: one cutting through the particle along the evaporation direction and another one pointing in a direction perpendicular to the latter. Whatever is the evaporation field of the particle, the size of the particle (Z), along the evaporation direction is found almost constant and equal to the expected one (i.e. a radius of 4 nm). This observation matches previous observations from atomic modeling [13,67]. For each profile in the perpendicular direction, we calculated the reduced density profile (ρ, atomic density normalized by the matrix one), a measurement of the average core concentration in atoms of the particle, and a calculation of the shape factor (S, square of the dimension along the perpendicular direction (i.e. X) divided by the Z dimension, i.e. (X/Z) 2 ). The shape factor is a descriptor for the morphological distortion [5]. All those measurements are reported in the figure 11.
As expected, no local magnification, no concentration bias, and no morphological distortion occur when the two phases have the same evaporation fields (F Ev,B = F Ev,A , X B = 100 at.%, ρ = 1, S = 1) ( figure 9). Otherwise, as the results show, we can distinguish three different principal cases, according to the evaporation field of the particle. The origins and the interpretations of these three behaviors will be based on the simulation results reported in the figure 12.
This case is referred in the literature as a high field particle case, the particle has a higher evaporation field than its surrounding matrix. Those particles exhibit no chemical bias in the core (X B = 100 at.%), but a shape dilatation (S > 1) surplus atomic density fluctuation, with an higher one specifically located at the interface (ρ > 1) and a lower in the core (ρ < 1) (figure 11). Those observations have already been made experimentally for MgO particles in a Au matrix [70], nitride dispersion strengthened steel particles [71] or even Al 3 (Zr 1−x Ti x ) particles in an α-Al matrix [72]. When this kind of particle merges at the surface, its evaporation rate is lower than the matrix. Consequently, the particle is protruding (with a small radius of curvature and so a high electric field compare to the matrix) and the matrix at the interface is exhibiting a very high radius of curvature. This interface can also have negative radii of curvature, since the concavity can be modified, when the evaporation field difference is very important. This surface topography induces divergent trajectories from the particle, as well as of the matrix at the interface, but to a lower extent. In effect, a size expansion of the particle is observed due to trajectory aberrations of the core (then a low atomic density of impact on the detector) and a trajectory overlap between atoms located at the interface (matrix and particle with then a high atomic density).
This case is referred to in the literature as the low field particle case. Those particles exhibit no chemical bias in the core (X B = 100 at.%), but a shape compression (S < 1) and a high atomic density (ρ > 1) (figure 11). Those observations have already been done experimentally for an α' particle in FeCr [73,74]. Their evaporation rate is higher than the matrix. In effect, the matrix at the interface is then becoming protruding and the particle develops a convex interface (negative radii) Figure 11. Reduced atomic density, measured through the particle (with an initial radius of 4 nm) in the direction perpendicular to the evaporation direction, according the distance from the particle core (nm), average core concentration in particle atoms (at.%) and the shape factor, as a function of the evaporation field ratio of the particle with the surrounding matrix (from 0.4 to 1.4 by 0.1). The dash lines in the reduced atomic density map represents the interface between the particle and the matrix, define at half height of the concentration in the particle core and the one in the matrix. and a 'flat' core (high radii), thus a low electric field. This topography induces on the one hand convergent trajectories from the particle and on the other hand divergent trajectories from the matrix at the interface. This divergence is not enough to overlap with the trajectories for ions originating from the particle (so no chemical bias [75]). In this simulated case, it should be noted that the difference in topography within the particle, between the core and the interface, causes a variation in the atomic density, both exhibiting an over atomic density, but a higher one at the interface (especially visible for F Ev,B = 0.9 × F Ev,A ). In both cases, the observed high atomic density is thus only induced by the morphological compression.
This latter case is not very widespread in the literature [1], since it is difficult to probe in experiment. Indeed, there is a chemical composition bias (X B < 100 at.%), an over atomic density (ρ > 1) and a morphological variation depending of the evaporation field: compression if 0.5 × F Ev,A < F Ev,B < 0.8 × F Ev,A (S < 1) and a dilation otherwise (S > 1). The common point in this evaporation field range is the particle trajectories crossing-over due to an intense convergence, caused by an entire convex particle (negative radii). Up to now, this phenomena has been observed experimentally in oxide dispersion strenghtened (ODS) steel particles by [64] and in SiGe layer embedded in SiO 2 [76]. In all these cases, the matrix interface divergence (induced by it protruding) is enough to overlap with particle trajectories (resulting in chemical bias [5]). If F Ev,B > 0.5 × F Ev,A , the cross-over occurs far from the tip. Consequently, the defocusing which is induced by the crossover is insufficient to compensate for the initial focusing. This results in a compressed particle (S < 1). This is the case of Cu particles in Fe [77]. As the FE decreases, the cross-over occurs closer to the tip leading to an apparent divergence of the trajectories (S < 1) which increases the overlap (so the chemical bias).
Our model allows to reproduce not only experimental observations but also results that were obtained by atomic simulation [1]. With our model, we observe that ion cross-over occurs for particle evaporation field lower than 0.7-0.8 × F Ev,A (figure 11). With an atomic model, Vurpillot observes that this phenomena occurs for particle evaporation field values lower than 0.5-0.7 × F Ev,A [1]. It is obvious that the particle size ratio with the tip size [5,64] is correlated with the evaporation field at which this crossing over occurs. In this case, our model presents the following advantages: (i) The model reproduces experimental and atomic scale simulation results in reasonable computing hour (less than one hour per case) and (ii) to simulate real sample dimension (with a radius of 50 nm), also in a reasonable time computing but mainly without any memory allocation issue (200 MB compared to 30 GB in atomic model [12]). The investigation of the size effect (of the particles and the tip) will be the subject of a future study.

Vertical layer
Our last studied case is a setup with a vertical layer studied for laser assisted evaporation. Our objective is to see if the model is capable of reproducing observations made by Lee et al [78] as well to confirm their interpretations if possible. They made TEM observations of a vertical layer Si embedded in SiO 2 (figure 13(a)) and its opposite (SiO 2 embedded in Si, figure 14(a)), after APT analysis. In the first case, they observe that the SiO 2 layer radius of curvature on either side is not the same (from 25 nm in the shadow side to 35 nm in the illuminated side). The authors attribute this to the laser assisted evaporation, considering that the laser side develops higher radii [60]. To test this, we simulate a vertical layer of Si (thickness of 20 and a sample radius , it is also reported the results obtained for another temperature evolution from 50 to 150 K (with α = 2 and β = 3). Reprinted from [78], Copyright (2014), with permission from Elsevier. of 25 nm [78]) with an evaporation field of Si 0.7 times lower than the surrounding SiO 2 . The simulations have been performed at constant temperature (T = 50 K, figure 12(a)) and with an evolutive temperature (from 50 to 200 K, with α = 2 and β = 3, figure 13(b)). Clearly, we note that the temperature simulation makes it possible to reproduce a surface morphology very similar to that of the experiment, which is not the case for constant temperature (same radii in both side for SiO 2 layers). We measure a variation of 11 nm in the radii of curvature on either side of the tip (at X = −20 nm, the radius is equal to 16 nm, figure 13(e), and at X = 20 nm, the radius is equal to 27 nm, figure 13(f)), thus confirming the conclusions of Lee et al. We also reproduce the difference of the evaporation depth for both SiO 2 layers (about 3-4 nm). Our results, however, are not quantitatively similar to the experimental ones. The variation of radii is quite the same (about 10 nm), but the simulated radii are lower (in the shadow side, 16 nm in simulation and 25 nm experimentally). In order to be quantitative, it will undoubtedly be necessary to increase the size of the tip, even to modify the evaporation field of the Si layer, but more fundamentally to model more accurately the temperature distribution in the tip apex (figures 13(d)-(h)).
For the opposite case (SiO 2 in Si, figure 14), we make similar observations. Indeed, we note that the temperature simulation (figure 14(c)) makes it possible to reproduce a surface morphology very similar to that of the experiment ( figure 14(a)), which is not the case for constant temperature ( figure 14(b)). However, in this case, our results are far from being quantitatively comparable with experiment. They measure an increase of 10 nm of the SiO 2 layer radius (from 35 to 45 nm, figure 14(a)), whereas in simulation, we measure an increase of 4 nm, in the worst case (i.e. at the highest temperature in the figure 14(d), for X = −6 and X = 6 nm). Moreover, we observe a difference in the Si layer top depth (about 3-4 nm), as experimentally ( figure 14(a)), but in the opposite direction (the uppermost layer is the laser side experimentally , it is also reported the results obtained for another temperature evolution from 50 to 150 K (with α = 2 and β = 3). Reprinted from [78], Copyright (2014), with permission from Elsevier. and inversely in simulation). This difference is due to a noncentered SiO 2 layer between the two layers of Si, as observed by [78] before evaporation, in this case. However, it is obvious regarding the simulation results that the initial sample morphology (and even the evaporation field and temperature) must be modified to quantitatively correspond with experimental observation, even if, the main observations are well reproduced with the mesoscopic model. The mesoscopic simulation model allows to reproduce the evaporation of layers (but also of particles). The model also reproduces the laser assisted evaporation with a simple approach (polynomial evolution). The next optimization of this model will be to introduce a temperature evolution which can take into account [41,47] the thermal diffusivity variation between different phases.

Conclusion
A new model of FE was presented. The approach is based on the classical relationship between electric field (estimated from the Robin equation) and an evaporation rate formalism that is applied this time on a triangulated surface mesh of an emitter. The model is compared to predictions from a previous model at the atomic scale. This model, despite its simple approach, has many advantages: (1) low memory allocation requirement (around 100 MB), (2) low computation time (around 1 h), (3) real sample dimension can be introduced (regarding previous advantages), (4) can simulate various microstructural objects (particles, layer, …), with results comparable to the experiments and (5) also to reproduce laser assisted evaporation. This last point remains to be optimized in order to quantitatively reproduce the laser assisted experimental results.
The model accurately solves the electric field at the surface, and the ion trajectory from the sample to the detector. Given its performance in computing time and memory allocation, this model could be used as a basis for a 3D reconstruction algorithm [19,20].

Data availability statement
The data cannot be made publicly available upon publication because they contain commercially sensitive information. The data that support the findings of this study are available upon reasonable request from the authors.