SIMULATION OF MIRROR ELECTRON MICROSCOPY CAUSTIC IMAGES IN THREE-DIMENSIONS

A full, three-dimensional (3D) ray tracing approach is developed to simulate the caustics visible in mirror electron microscopy (MEM). The method reproduces MEM image contrast resulting from 3D surface relief. To illustrate the potential of the simulation methods, we study the evolution of crater contrast associated with a movie of GaAs structures generated by the droplet epitaxy technique. Speci¯cally, we simulate the image contrast resulting from both a precursor stage and the ¯nal crater morphology which is consistent with an inverted pyramid consisting of (111) facet walls. The method therefore facilities the study of how self-assembled quantum structures evolve with time and, in particular, the development of anisotropic features including faceting.

Thi s v e r sio n is b ei n g m a d e a v ail a bl e in a c c o r d a n c e wit h p u blis h e r p olici e s. S e e h t t p://o r c a . cf. a c. u k/ p olici e s. h t ml fo r u s a g e p olici e s. Co py ri g h t a n d m o r al ri g h t s fo r p u blic a tio n s m a d e a v ail a bl e in ORCA a r e r e t ai n e d by t h e c o py ri g h t h ol d e r s .

Introduction
Mirror electron microscopy (MEM) is a well-established technique which has a number of advantages for the study of nanostructure formation. [1][2][3][4][5][6][7][8] Electrons neither impact nor are emitted from the specimen surface. Instead, a normally incident beam is re°ected at equipotential surfaces just above the specimen. This is achieved by holding the specimen at a small negative voltage relative to the electron source. In the turnaround region, the electrons are traveling very slowly and are sensitive to spatial and/or temporal variations in electric¯eld which can be produced by changes in surface height or work function across the specimen surface. This results in the de°ection of electrons, which redistributes their position in the image plane, producing image contrast. Since the electron beam does not impact the surface, it can be applied to study sensitive specimens. In addition, the parallel nature of the technique facilitates the acquisition of real-time movies of surface evolution. [9][10][11][12] The re°ected electrons in MEM indirectly contain information related to the topography and/or the electrical and magnetic properties of the surface. This has stimulated signi¯cant e®orts to interpret MEM image contrast and extract quantitative information regarding the surface properties. 1,2,6,8,[13][14][15][16][17][18][19][20][21][22][23][24] In general, MEM image contrast can be highly non-intuitive since it arises from electric¯eld variations above the specimen. In the special case of weak electron de°ections, it can be shown that the image contrast results from the Laplacian of small height or potential variations across a sample surface, where the contrast is blurred slightly to account for the interaction of the electrons with the electrical potential away from the surface. 25,26 This is a signi¯cant simpli¯cation which allows rapid interpretation of MEM images.
However, in general, for larger de°ection of electrons, such as those arising from liquid droplets or quantum structures, the images consist of envelopes of electron rays or caustics which complicate the interpretation of the contrast. In order to link surface morphology and MEM image contrast, it has therefore been necessary to apply electron ray tracing methods to simulate caustic features. This so-called caustic imaging has been successful in explaining experimental MEM contrast from droplets 27 and nanowires. 24 In addition, it has been applied to simulate image contrast during droplet epitaxy, where cylindrical symmetry was assumed. 28 While the simulations explained the salient features of the images, image contrast associated with craters appearing in the quantum structure displayed approximately four-fold symmetry. This is clearly linked to surface energy anisotropy and faceting and cannot be reproduced by the cylindrical symmetry of the simulations. The purpose of this paper is, therefore, to fully extend caustic imaging theory to 3D. In principle, this method facilitates the simulation of low symmetry surface features. To illustrate the potential of the simulation methods, we study crater contrast associated with GaAs structures generated by the droplet epitaxy technique. [29][30][31][32][33][34][35][36] 2. MEM Imaging Geometry Figure 1 displays the imaging geometry associated with MEM. A converging electron beam of energy U passes through a grounded anode aperture A and emerges parallel to the optical axis z. The cathode of the immersion objective lens is formed by a quantum structure specimen located at z ¼ L and held at a negative potential V . The electron beam reverses in direction at z ¼ L M , a distance from the cathode surface such that, where Àe is the electronic charge. The electron beam is then reaccelerated into the imaging system of the microscope following de°ection by the electric¯eld surrounding the quantum structure surface. Contrast in the MEM image results from the redistribution of electrons on the virtual image plane at z ¼ Áf þ 4L M =3. 18,37 Here, the defocus distance Áf is controlled by the magnetic part of the objective lens.

3D Electrical Potential
To simulate the MEM image contrast from a general surface morphology, it is¯rst necessary to evaluate passing through an aperture in the grounded anode Aðz ¼ 0Þ and is de°ected slightly as the beam passes through the anode, resulting in parallel illumination of the sample. 18,37 The electron beam turns around in the vicinity of the turning distance z ¼ L M ¼ L À , for some small distance from the cathode C. The electron beam is de°ected upon interacting with the electric¯eld above the cathode surface, and is reaccelerated away from the cathode passing back through the anode aperture Aðz ¼ 0Þ. The microscope is assumed to form an image of the electron positions as they would appear on a virtual image plane at z ¼ Áf þ 4L M =3, 18,37 with positions ðx; yÞ found by tracing back along the dashed blue lines.
Here Áf is the defocus distance from the plane z ¼ 4L M =3, and is controlled by the magnetic part of the objective lens.
The y axis extends out of the page as shown.
the electric potential in the region above the surface. This is accomplished by solving Laplace's equation using¯nite element methods with the specimen topography as one boundary and the grounded anode as the opposite boundary. Here, the anode is approximated as a planar, apertureless plate with the diverging e®ects of the aperture explicitly included in the ray tracing. 18,27,37 We utilize the FreeFEM++ nite element package with a 3D mesh. 38 Figure 2(a) displays a 3D height map of a quantum structure with cross-sections in the y ¼ 0 and x ¼ y directions shown in Figs. 2(b) and 2(c). The equipotential surfaces for the two cross-sections are shown Figs. 2(b) and 2(c) as a function of distance H from the cathode for the experimental parameters U ¼ 20 keV, V ¼À20000:4 V, and L ¼ 2 mm, and clearly take on the symmetry of the quantum structure.

3D Ray Tracing
In order to calculate the MEM image intensity, a family of electron ray trajectories is traced through the electric potential using a fourth-order Runge-Kutta method. 27,39 The incident electron paths begin at z ¼ 0 with a square grid of rays of spacing x 0 ¼ 6000=256 ¼ 23:4375 nm in both the x and y dimensions. These initially parallel rays are traced through the turn-around region in the vicinity of the quantum structure and back to the anode aperture. Treating the aperture as a diverging lens, 18,27,37 the emerging rays are then projected back along straight line paths to the virtual image plane at z ¼ Áf þ 4L M =3 as shown in Fig. 1. This image, demagni¯ed by 2/3, is then the object for the objective lens.
To simulate the image intensity in the two dimensional case, with the electrons positioned along the x-axis only, we were able to compare the spacing between two initially adjacent electrons on the virtual image plane to that expected for the ideal case of a°a t equipotential specimen. 28 A closer spacing implies an increase in intensity as more electrons are focused in that region, whereas a larger separation implies a reduction in intensity. In the 3D case, the electron positions are spread across the virtual image plane in two dimensions, but we can evaluate the intensity using similar reasoning. The Voronoi region of an electron position is the space on the virtual image plane that is closer to that electron than any other electron, as illustrated in Fig. 3. 39 In Fig. 3, panel (a) depicts a square grid of equally spaced electron positions that forms the input at z ¼ 0. Panel (b) shows an example of the electron positions on the virtual image plane for an ideal°at equipotential specimen, where the electron positions remain in a square grid but usually with a di®erent spacing. The shading represents the Voronoi region of the central electron position, formed by bisecting the lines connecting the electron position to all nearby electrons. The vertices that de¯ne this Voronoi region are indicated by the Â symbol, and the Voronoi region in this ideal case is a square for all but the outermost points. Panel (c) gives an example of the electron positions for a sample that scatters non-uniformly, with the Voronoi region for an example electron position labeled`2' shaded.
The area of the Voronoi region sðx i ; y i ;;ÁfÞ around electron i's position ðx i ; y i Þ provides an e®ective area occupied by that electron and no other, and can be compared to the expected area for an ideal°a t equipotential specimen to simulate intensity. Where s is smaller than the ideal case the intensity is increased, and conversely where s is larger the intensity is decreased. After simulating the electron positions on the virtual image plane for a particular defocus Áf, we used the program qHull 40 to compute the Voronoi regions around each electron's position. Each Voronoi region is de¯ned as a polygon with N vertices ðx 1 ; y 1 Þ; ...; ðx N ; y N Þ, and the area of a particular Voronoi region may be computed via a Jacobian giving 39 where N þ 1 also refers to vertex 1, and abs represents taking the absolute value of the signed area as we only need its magnitude. 39 With an input square grid of initial spacing x 0 in the x and y dimensions, the electron positions for an ideal°a t equipotential surface are separated by x 0 ð 2 3 À Áf 4L M Þ in the x and y dimensions. 27 In this case, the Voronoi region of each electron position is a square centered on each point (see Fig. 3 The intensity at a particular electron position for the non-uniformly scattering case with Voronoi region area s is calculated by comparing this area s with that expected for the ideal case, By repeating this process, an intensity value was generated for each electron position giving a 2D intensity simulation. Where electron positions become very close the Voronoi region can become very small (i.e. s ! 0), creating caustics in the image as electron trajectories overlap. This was treated by choosing a threshold area, e.g. 0:1, below which this threshold value was assigned to s therefore limiting I via Eq. (4). This is equivalent to specifying the saturation level of the detector.

3D Image Simulations
We now apply 3D caustic imaging theory to simulate MEM movies of central crater formation in droplet epitaxy. 9 The movies were obtained using a III-V epitaxy LEEM. 41 An undoped GaAs (001) epi-ready wafer was degassed at 300 C under ultrahigh vacuum for 24 h. The surface oxide was then removed by high temperature°ashing up to 600 C and annealing at 580 C for 2 h. Ga droplets were then formed by annealing above the congruent evaporation temperature at 650 C. The sample temperature was reduced to 460 C and images were recorded in MEM mode. [1][2][3][4][5][6][7][8]27 A Ga droplet was exposed to an As 4°u x of beam equivalent pressure (BEP) 1:45 Â 10 À5 Torr and the droplet crystallized into solid GaAs. Many of the morphological features associated with structure have approximate cylindrical symmetry but, of particular interest to this study, is the faceted crater which forms between 15 to 20 minutes of exposure to As 4°u x. MEM images taken from a movie of central crater formation are shown in Figs. 4(a) and 4(b). The precursor state in (a) occurs after 15 min and still contains a Ga droplet at its center. 9,28 The MEM contrast consists of four bright corner spots with faint streaks emanating from the corners. Note that the  40 The area of this region may be compared to that expected for the ideal case to calculate the intensity at that electron position. asymmetry in spot positions suggest some asymmetry in the way in which the structure is evolving with time. The image after 20 min in panel (b) corresponds to the¯nal central crater morphology and consists of a bright central spot with a four-fold pattern of radial streaks emanating from the center. 9 To aid interpretation of the MEM contrast, additional atomic force microscopy (AFM) experiments were performed in which Ga droplets were exposed to As 4 at identical°u x and temperature for 20 min and then quenched to room temperature. This resulted in the AFM height pro¯le maps (panels (c), (d)) and accompanying line traces (panels (e), (f)). 9,28 We note that quenching to room temperature may induce artefacts and the observed morphologies may not exactly re°ect the shapes undergoing droplet epitaxy at 460 C. The AFM data contained in Fig. 4 can therefore only be used as an approximate guide to the surface shape under actual growth conditions. Note that panels (a) and (b) were taken from a MEM movie of droplet epitaxy 9 so that the images correspond to the same structure.
To apply 3D caustic imaging theory to simulate MEM movies of central crater formation, we used surface pro¯les generated from the appropriately scaled AFM data in Fig. 4 as an initial starting point. MEM image simulations were used to further¯netune the surface features iteratively to obtain a best t to the experimental MEM data in Figs. 4(a) and 4(b). A defocus of Áf ¼À42 m was assumed based on simulation of droplet contrast before exposure to As 4°u x. 27,28 The resulting surface pro¯le generated for the precursor state is shown in Fig. 5. Note that in these simulations, liquid Ga is still part of the quantum structure and a work function di®erence of 0.3 V is applied between liquid Ga and GaAs. 27,42 The surface pro¯le generated for the¯nal central crater morphology is contained in Fig. 2.
The MEM image simulations generated from the surface pro¯les shown in Figs. 2 and 5 are shown in Figs. 6(a) and 6(b), respectively. These should be compared directly with the experimental MEM images in Figs. 4(a) and 4(b). In the case of the precursor state, the simulation reproduces the bright corner spots evident in Fig. 4(a) which are caused by electrons being focused into the corner depressions in the developing crater wall (Fig. 5). The fainter square feature linking the four corner spots is a result of the height discontinuity between the surface and the liquid Ga (seen in Fig. 5(b)). Although small, this height decrease is enough to weakly focus the electrons rays, and its square shape (with <110> sides) results in the square feature in the simulated image. However, this feature is too faint to be discerned in the experimental MEM image in Fig. 4(a) where only the bright corner spots are clearly visible.
Similarly, the simulation in panel 6(b) captures the appearance of a central bright spot as the crater acts as an electron lens, focusing the electrons to a caustic. 28 The four-fold radial streaks are also reproduced re°ecting the faceted nature of the crater. Based on the AFM data and image simulations, wē nd that a (111) plane facet angle (i.e. 54.7 from the vertical) is consistent with the observed contrast. This also agrees with previous work on faceting during quantum dot formation by droplet epitaxy. 43 From the MEM simulations in Fig. 6, it can be seen that the bright concentric rings are associated with discontinuities in the surface pro¯les in Figs. 2 and 5. Changes in the surface height pro¯le, such as discontinuities, modify the equipotential surfaces above the surface (see, for example Fig. 2). Although the equipotential surfaces resulting from surface discontinuities smooth with increasing distance from the surface, electron trajectories from either side of these discontinuities are de°ected in di®erent directions which causes electron paths to overlap in the returning beam. A projection of these emerging rays back to the virtual image plane at z ¼ Áf þ 4L M =3 results in bright caustic features in the images. 28 The ring structures associated with cylindrically symmetric discontinuities are particularly striking features of the simulations (Fig. 6). However, the outer bright concentric ring present in the simulations is not as pronounced in the images. This is probably due to intrinsic roughness of the surface in this region which lowers the experimental contrast. The four-fold patterns near the center of the images in Figs. 4(a) and 4(b) can also be related to the focusing e®ect of surface discontinuities shown in Figs. 2 and 5. This is clearly linked to surface energy anisotropy and faceting in the case of the t ¼ 20 min central crater. Such features can only be reproduced by full 3D simulation methods.
Additionally, one can see that the simulations in Fig. 6 reproduce the successive reduction in mean intensity levels from the planar surface to within the outer and inner rings of the structure. This can be explained by the increasing average slope of the surface morphology as the center is approached, until the crater is reached. This increasing average slope de°ects electrons through greater angles via the equipotential surfaces, thus locally reducing the intensity in the virtual image plane for the chosen defocus.
By adding appropriate shifts to the position of rays in the virtual image plane, it is possible to incorporate spherical aberration. 27 Similarly, a weighted average of a series of monochromatic intensity patterns for a spread of energy values can be applied to include chromatic aberration. 27 However, we¯nd both aberrations have a negligible e®ect for the relatively low resolution case considered here given C S % 0:1 m and a Gaussian energy spread of full-width half-maximum 0.3 eV.
Limitations of the current methods can be observed as weak background¯ne structure in Fig. 6 which results from errors in the electric potential calculated with¯nite element methods over a coarse 3D grid. The simulations can be improved at the expense of more computational power by using a¯ner 3D grid and/or adaptive meshing to more accurately evaluate the electric potential above the specimen. An improved mapping of the electric potential increases the accuracy of the simulated electron paths and hence image intensity. Nevertheless, the methods employed here are adequate to qualitatively simulate the salient features of the images.

Conclusions
The ability to perform full 3D simulations of MEM contrast opens up the possibility of studying the realtime evolution of faceting and growth anisotropy during nanostructure fabrication. This has been demonstrated by simulating crater contrast associated with GaAs structures generated by the droplet epitaxy technique. The simulated contrast is consistent with an inverted pyramid consisting of (111) facet walls. More generally, the simulation methods developed here can be used to model MEM contrast originating from surface topography, electrical and magnetic¯elds associated with low symmetry geometries. When applied to data sets obtained at several values of defocus (see, for example Refs. 24 and 44) and/or sample potential, this may form the basis for quantitative extraction of topology or¯elds from the images.