Finite-difference time-domain simulation of cathodoluminescence patterns of ZnO hexagonal microrods

The Finite-Difference Time-Domain (FDTD) numerical simulation method has been applied to interpret cathodoluminecence patterns observed for ZnO nanorods grown by a hydrothermal method. The 3D FDTD simulation reproduced the radial electromagnetic field pattern in the hexagonal resonator, corresponding to the CL emission maps of real ZnO microrods. The simulation result for the H z (TE) polarization—the intense field distribution along edges of the structure, in particular in the corners, but weak in the centre—matched the CL pattern particularly well. Since the experiment was not polarization sensitive, we suppose that polarisation sensitive transmission of electromagnetic field through the ZnO/air interface leads to such an observation. The results of the simulation show also that the lack of axial Fabry-Pérot-like resonances in the CL experiments is caused by leaking of the electromagnetic field from the ZnO resonator into the GaN substrate.


Introduction
In our everyday life, light emitting diodes (LEDs) or laser diodes (LDs) made of wide band gap semiconductors widely replace conventional light sources. Nevertheless, they still suffer from limited internal quantum efficiency, low extraction efficiency etc. Many of those problems are related to structural imperfections (like high dislocation density) occurring at the interfaces between the layers of structurally dissimilar materials (with different crystallographic structure, different lattice parameters etc.). One of considered solutions of at least some of those problems is replacing continuous thin layers of semiconductors with ensembles of quasi-1D nanowires (NW) grown perpendicularly to the substrate. The NWs may contain axial or radial heterostructures necessary to create efficient light-emitting devices, lasers etc [1,2]. Merging together advantageous physical properties of semiconductors, very well developed technology of semiconducting structures and benefits of nanostructurization paves the way for fabrication of nano-devices potentially better than those made with use of planar technology. So, optical and electronic properties of quasi one-dimensional objects are investigated in view of optimising their preparation conditions and crucial parameters.
In quasi-1D structures the strains at the interfaces occurring due to the lattice mismatch can be relatively easily accommodated. This facilitates growth of complex 1D structures consisting of layers made of different materials and becomes an important advantage for fabrication of the devices embedded in 1D structures. Moreover, a quasi-1D object with regular shape defined by the crystalline structure and the size comparable with the wavelength of light can confine electromagnetic radiation and play a role of a resonator for a nanodevice embedded in a microrod. In such a system, for example, extremely strong exciton-photon coupling may occur, even at room temperature [3,4]. Therefore, nanostructures with various shapes and sizes are grown and investigated in view of selecting the optimum preparation method of nanoresonators suitable for possible direct use in photonic devices.
Among many fields of investigations of ZnO nanoobjects, a considerable interest is focused on their behaviour as resonant cavities supporting whispering gallery modes (WGMs). Such effects were observed in microdisks [15], micro/nanoneedles [16,19], nanowires [18], etc. These studies are aimed at the enhancement of luminescence from the structures.
A ZnO microrod growing along the hexagonal axis of its wurtzite crystal structure has the regular hexagonal cross-section and may naturally serve as a small-sized optical resonator. The light can circulate around this microcavity in a hexagonal loop path due to the multiple total internal reflection (TIR) from the six side walls of the microstructure acting as reflective mirrors and supporting the optical gain. In this geometry an angle of incidence is 60°, much exceeding the critical angle of ZnO (25°). The light is confined inside of the ZnO microstructure and the optical losses are considerably reduced. The single ZnO microrod with well-defined cavity geometry can support a set of transverse cavity modes and can also act as a Fabry-Pérot resonator both along its c axis and between its opposite side walls [20]. So, the proper observation and interpretation of the electromagnetic field distribution in the nano/microrod is crucial for assessment of its usability as a part of an optical system.
Cathodoluminescence is a continuously developing experimental technique which enables us to observe, map and analyse optical excitations in materials excited with a focused electron beam [21][22][23][24][25]. Retaining the submicrometer lateral resolution from the scanning electron microscopy and high spectral resolution from optical spectroscopy techniques, it gives us information about spectrum of luminescence from various spots on an individual nano/microobject. In particular, it turned out to be suitable for revealing optical resonances in such systems [15,[18][19][20]26].
Proper interpretation of the images obtained with use of cathodoluminescence mapping needs a method for reproducing the observed pattern by first-principles calculations or at least by numerical modelling. Tracking the light rays reflected from the side walls of the resonator by the methods of the ray optics turned out to be an effective way of modelling luminescence spectra inside such a system if the size of the resonator is much bigger that the light wavelength. The resonance wavelengths corresponding to WGMs can be analytically calculated assuming the total internal reflection from the side walls [15,[27][28][29][30][31]. Some roughness of the walls can also be taken into account [32]. In the opposite situation, for systems much smaller than the radiation wavelength, quasi-static approximation gives acceptable solutions at a reasonable computational cost. Description of intermediate situations needs more sophisticated methods, like numerical solving of the full wave equations used in [33] or numerical solution of the 3D Poisson's equation, current continuity equations and rate equations using finite element method [34].
The FDTD method can also be used in modelling optical systems with basic size parameters comparable with the light wavelength. It is conceptually simple, relatively easy to implement, although usually computationally expensive. Nevertheless, for relatively small and simple systems application of this method becomes reasonable and gives useful results [35]. Relatively easy adjustment of the parameters describing the studied system makes this method more convenient for simulation of real experimental data than the results of more advanced calculations available in literature. Some ZnO nanoobjects ware also simulated with use of FDTD method [28,[36][37][38][39].
In this paper, we report on using the FDTD method to interpret the results of cathodoluminescence study of hexagonal microrods of ZnO obtained in a hydrothermal process. The features revealed in the cathodoluminescence spectrum can be interpreted as a manifestation of whispering gallery modes in a hexagonal resonator, while their energies can be calculated using simple ray-optics model assuming total internal reflection from the side walls [40]. However, the patterns obtained in the CL imaging mode (cathodoluminescence intensity maps) cannot be directly reproduced within such a simple model. Therefore we applied the FDTD method to model electromagnetic field distribution in a 2D hexagonal resonator as well as in a 3D hexagonal prism. We assumed that excitation with the electron beam of the sample regions corresponding to the maxima (antinods) of electromagnetic field in a resonator leads to stronger cathodoluminescence emission. So, the CL intensity maps correspond to the electromagnetic field distribution in the system. We applied the FDTD method implementation developed and published by John B. Schneider [41]. This approach enabled us to reproduce the cathodoluminescence intensity pattern for the wavelength/microrod-diameter ratio corresponding to that observed in the CL experiment. The simulation was also used to find the system feature responsible for the lack of axial Fabry-Pérot resonances. The results suggest that disconnecting the microrods from the GaN substrate would reduce the leaking of the electromagnetic field from them. Then, the axial resonances should be more probable, at least for the wavelengths corresponding to the low absorption of ZnO.

Cathodoluminescence of ZnO microrods grown by a hydrothermal method
The ZnO microrods are grown by a variant of the microwave-assisted hydrothermal method [42,43] which enables us to obtain the structures with a high structural quality (confirmed by TEM studies [43]), after a rapid process (1-3 min.) carried out at relatively low temperatures (below 100°C). The ZnO nanorods can be grown on variety of substrates (like Si, GaN, ZnO, glass) while conditions of the chemical process and substrate properties determine size and orientation of the nanowires. Growth of the ZnO microrods proceeds in the aqueous solution. The reaction mixture consists of deionized water, zinc acetate dihydrate and sodium hydroxide (used as pH regulator) from Sigma-Aldrich Co. Zinc acetate was used as a zinc precursor by other groups for nanorods growth [44,45], however, in contrast to our method, other phases of zinc (Zn(OH) 2 and ZnO(OH) 2-) were employed. In our approach, a more reactive Zn(OH) + phase is used, thus dynamic chemical reaction proceeds instead of slow crystallization process from saturated solution. This gives us so high growth rate. Since uniform heating of the solution is very important, a microwave-assisted version of hydrothermal reactor is used. On GaN substrates (commercial templates from Saint Gobain Lumilog) ZnO nano/microrods grow as hexagonal prisms with the (111) axis perpendicular to the substrate surface. The lattice constant misfit between ZnO and GaN is about 1.8%. The microrods shown in figure 1(A) are obtained with the optimum pH value of reaction solution equal to 7 and the zinc ion concentration of 0.13 mol dm −3 . The growth is carried out in 2-3 minutes at the temperature of 70 0 and under atmospheric pressure.
The morphology and size of the so-grown ZnO microrods are assessed by scanning electron microscope (SEM) Hitachi SU-70 using secondary electron detector. Light emission properties of the ZnO microrods are investigated by CL spectroscopy and imaging with use of SEM Hitachi SU-70 equipped with the Gatan Mono CL3 system. The monochromatic CL map of a typical ZnO microrod is shown in figure 1(B). The measurements are carried out at temperature of 5 K and for the accelerating voltage from 5 to 15 kV, the beam current from 2.4 to 14 nA. The emission intensity is not distributed homogeneously across the microrod, but is concentrated near the hexagonal boundary, especially at the six corners of the nanorod. Figure 1(C) shows a CL spectrum collected for the ZnO microrod with a diameter of 1.61 μm. The electron beam excites the microrod at the center of its top face. The spectrum covers the band of defect-related emission. [40] contains more examples of CL spectra collected for the ZnO nanorods with different diameters. The spectra corresponded to the near-band-edge emission. The multi-peak structure was interpreted in terms of ray optics as whispering gallery modes of the hexagonal resonator. An alternative interpretation taking into account phonon replicas of band-edge transitions was also considered and regarded as unsatisfactory.

FDTD simulations
The finite-difference time-domain method enables us to model electrodynamic systems described by Maxwell's equations. It is based on the Yee algorithm [46]. We replace the differential Maxwell's equations with corresponding finite-difference equations and we solve them in the discretized space and time. We define the system by ascribing relevant parameters describing material properties (like relative permeability μ, relative permittivity ò) to each point of the discretized space. Then we use the finite-difference equations to calculate alternately the electric and magnetic fields for next future time moment starting from the values known for the past moment. As a result, we observe time evolution of the fields in the investigated system.
The procedures and codes applied to simulate electromagnetic fields in ZnO nanorods are based on those given in [41] for 2D and 3D simulations. We use second-order absorbing boundary conditions on the edges or walls of the finite-size box in which simulation was carried out to mimick distribution of electromagnetic field in the ZnO nanorod surrounded by an infinite space. To excite the system, we apply a harmonic source or the Ricker wavelet. The later enables us to simulate propagation of the fields in a wide frequency range.

Results and discussion
4.1. 2D resonator Taking into account the axial symmetry of the ZnO nanorods, we can try to simplify the simulation assuming that the radial distribution of the electromagnetic field is governed only by the shape of the nanorod crosssection. This enables us to reduce the dimensionality of the problem and start the analysis of the electromagnetic field distribution in ZnO nanorods from simulation of a 2D ZnO resonator with hexagonal shape.
The simulation space is a square of the size of 2000 × 2000 points with the second-order absorbing boundary conditions at the edges. The test simulation of the empty space with a harmonic point source at the center proved that such configuration effectively prevents reflection of the electromagnetic waves from the edges of the simulation space. So, no parasitic pattern interferes with that caused by the presence of the ZnO resonator.
The ZnO resonator is defined by ascribing the high frequency relative permittivity (dielectric constant) of zinc oxide ε ∞ = 3.75 to the points inside the hexagon with the sides of 200 points (figure 2). ZnO is birefringent and ε(∞) ⊥ is lower than ε(∞) P by 0.05 [47]. However, such difference does not manifest itself in the results of the simulation. The relative permeability μ is kept equal to 1.
The CL patterns of ZnO nanorods were measured with the electron beam scanned over the whole top surface of the nanorod. Therefore, the simulation is carried out with all points of the resonator excited in-phase by the harmonic (sinusoidal) sources. The frequency of the sources is chosen so that the D/λ ratio is equal to 2.9. This value corresponds to the ratio of the parameters describing the investigated real system ( figure 1(B)). The diameter of the ZnO nanorod was equal to D = 1.61 μm and the CL wavelength for the maximum of defectrelated luminescence was λ = 550 nm. The results of the simulation are shown in figure 2. The patterns correspond to E z or H z field component squared and added up over a full period of exciting harmonic source, in order to asses the distribution of nodes and antinodes of the standing electromagnetic waves appearing in the resonator. In order to facilitate comparison of those patterns with the results of the CL experiment the patterns are blurred by convolution with the Gaussian function. The images shown in figures 2(B) and (D) are obtained for the σ parameter corresponding to approx. 7% of the resonator diameter. Neither of these patterns matches the CL image with the emission intensity concentrated at the edges and in the corners of the hexagonal cross-section of the nanorod. The patterns resulting from the 2D simulation have the antinodes of the electromagnetic field distributed over the whole resonator surface, consistently with the symmetry of the system. Therefore, an attempt to simulate electromagnetic field distribution in a hexagonal prism, more accurately reproducing the structure of ZnO nanorods, has been undertaken.

3D resonator
The 3D simulation space is in the form of a square prism with the height of 302 points and the 402 × 402 points cross-section ( figure 3). The ZnO nanorod is defined as a hexagonal prism with the values for relative permittivity and permeability of ZnO ascribed to all points of the simulation space inside it. Non-zero values for the loss factor are also introduced to mimic absorption of ZnO. The prism height equals 230 points and the hexagon edge is of 110 points. The GaN substrate, on which ZnO nanorods have been grown, is modelled in a similar way with use of the values for relative permittivity and permeability of GaN. The scheme in figure 3 shows also two cross-section planes (A and B) at which the electromagnetic field distributions are shown below to visualise the results of the simulations.
In order to check if the size of the simulation space may influence the results of simulations, some of them are also carried out for the same system embedded in a larger simulation space: 302 × 500 × 500 points. As no meaningful difference has been observed, the smaller system is used for all simulations discussed below to make them computationally less expensive.
The first aim of the 3D simulation is to check if any resonance-like features appear in the electromagnetic field distribution of investigated system. The CL experiment, which gave the result shown in figure 1(C), consisted in spectral analysis of luminescence collected above the nanorod excited in the centre with the electron beam. So, we start with a simulation of the ZnO nanorod excited in the centre with a pulsed source of the Ricker wavelet characteristics. The integrated field intensity above the top face of the nanorod is recorded as a function of time. Use of the Ricker-wavelet pulsed source allows us to model the response of the system in a wide range of frequencies.
The H z time dependence is converted into the frequency (wavelength) domain by a Fourier transform ( figure 4). The wavelength scale is constructed assuming that the D/λ ratio in the simulation equals to that of the real microrod of D = 1.61 μm. Analysis of other field components (e.g. E z ) gives qualitatively similar results with the same peak positions, but some smooth background may appear.
Since luminescence bands correspond to spectral ranges of considerable absorption (interband or defectrelated one), the loss factor has been introduced to the equations describing electromagnetic field evolution. The simulations are carried out with increasing ZnO loss factor in order to assess its influence on the shape of spectral and spatial field distribution. The overall shape of the simulated spectrum does not reproduce the maximum of CL shown in figure 1(C) because the exact shape of the spectral dependence of ZnO luminescence is not taken into account. The relative permeability μ, the relative permittivity ò and the loss factor are independent from the wave frequency.
The shape of the H z spectrum shown in figure 4, consisting of the set of distinct maxima, suggests that resonances occur in the simulated system. The overall shape of the simulated spectrum does not reproduce the maximum of CL shown in figure 1(C) because the exact shape of the spectral dependence of ZnO luminescence is not taken into account. The relative permeability μ, the relative permittivity ò and the loss factor are independent from the wave frequency. The number of the maxima for the wavelength range of 480-700 nm (9 in the simulation) is reasonably close to 11 resonance features observed in the experiment even if their wavelength values cannot be precisely reproduced. An increase of the loss factor leads to blurring of the spectral features. They are gradually coalescing into wider features dominated by every third maximum. To confirm that those features correspond to the subsequent resonances, we determined the E z 2 (TM) and H z 2 (TE) patterns for the ZnO loss factor of 0.001 (shown in figure 4 above the corresponding maxima). This value for the loss factor is chosen after comparison of the results for 0.00, 0.001, 0.002, and 0.003. It gives the clearest patterns, in particular for H z 2 . The patterns were taken along the plane marked in figure 3 as the cross-section A. The patterns are obtained with harmonic excitation of the wavelength corresponding to the selected resonance peak. The spatial profile of the excitation assumed in the simulation mimics the excitation profile induced by the electron beam impinging the top face of a ZnO nanorod during the CL experiment (the profile of energy deposited in the sample by scattered primary electrons). The maximum excitation is set at the height of 80% of the nanorod, than it decays exponentially towards the top of the nanorod. The whole cross-section of the nanorod is excited in-phase. The presented patterns are distributions of the squared field component accumulated over a simulation 'time' period corresponding to the period of the harmonic source of stimulation.
The observed substantial differences between E z 2 and H z 2 patterns can be interpreted as resulting from the differences between the phase shifts of the TE and TM components at the internal reflection from the walls of the resonator. So, the two components meet different resonance conditions and this leads to different field distribution patterns. The maxima correspond alternately to the TE an TM resonances, like in the model based on the ray optics [40]. It is particularly clear for the E z 2 patterns exhibiting regular, systematically changing shape of radial distribution for λ = 825 nm, 485 nm, and 363 nm, with increasing number of antinode shells.
An overall comparison of the patterns obtained by FDTD simulations with the results of CL experiments ( figure 1(B) and [40]) shows that only the distribution of the H z (TE) field component reproduces satisfactorily the experimental data. It seems to indicate that the experiment is sensitive mainly to the TE component of the electromagnetic radiation propagating in the studied ZnO nanorods. This can be the result of stronger leaking of this component from the resonator but this effect cannot be assessed convincingly within the reported simulation. Figure 5 shows the H z 2 pattern (at the cross-section A) obtained for λ = 405 nm and broadened by convolution with a Gaussian profile, assuming the value for the σ parameter corresponding to approx. 13% of the resonator diameter. Such broadening is intended to cover the most important factors determining blurring of the CL patterns: primary electron scattering in the nanorod material and excitation diffusion. The resulting pattern reproduces well the CL pattern obtained for the ZnO nanorod ( figure 1(B)).
Apart from the resonances based on reflection from the side walls of nanorods, one may expect that axial Fabry-Pérot-like resonances occur in nanorods of the length of 2-3 μm (due to reflection from the top and bottom hexagonal walls). Unfortunately, the Fabry-Pérot-like resonances do not manifest themselves in our CL experiments. In order to check if such phenomenon can occur on ZnO nanorods grown on a GaN substrate, the electromagnetic field axial distribution for λ = 485 nm is reproduced for a ZnO microrod attached to the GaN substrate and for a free-standing one. Figure 6 shows the results. The intensity of the field along the microrod is markedly stronger in the free-standing object. It can be interpreted as a result of the reduced ε difference between ZnO and GaN. It facilitates leaking of the electromagnetic field from the microrod into the substrate. Such a leakage of the field prevents formation of the standing wave along the microrod and can account for the impossibility to observe Fabry-Pérot-like resonances in our CL experiments.

Summary and Conclusions
We applied FDTD method of numerical modelling to reproduce the electromagnetic field distribution in hexagonal resonators made of medium with optical parameters simulating ZnO. The results were used to interpret the cathodoluminescence maps acquired for ZnO microrods grown by a hydrothermal technique.  A simplified model of 2D ZnO resonator failed to account for the shape of the radial distribution of CL emission in ZnO microrods. It was not possible to reproduce the intense field distribution along edges of the structure, in particular in the corners, but not in the centre.
The 3D FDTD simulation enabled us to reproduce the radial electromagnetic field pattern corresponding to the CL emission maps of real ZnO microrods. The simulation result for the H z (TE) polarization matched the CL pattern particularly well although the experiment was not polarization sensitive. We suppose that polarisation sensitive transmission of electromagnetic field through the ZnO/air interface may lead to such an observation.
The results of the simulation also account for the lack of axial Fabry-Pérot-like resonances in the CL experiments. They suggest that disconnecting the microrods from the GaN substrate would reduce the leaking of the electromagnetic field from them. Then, the axial resonances should be more probable, at least for the wavelengths corresponding to the range of low absorption of ZnO.
The reported FDTD modelling proved that resonant effects in hexagonal ZnO microrods can be regarded as a reason for the observed regular, symmetric but not uniform patterns of CL emission from such nano/ microstructures and should be considered in line with other interpretations, like uneven distribution of defects or impurities.