Field enhancement of multiphoton induced luminescence processes in ZnO nanorods

The near-ultraviolet photoluminescence of ZnO nanorods induced by multiphoton absorption of unamplified Ti:sapphire pulses is investigated. Power dependence measurements have been conducted with an adaptation of the ultrashort pulse characterization method of interferometric frequency-resolved optical gating. These measurements enable the separation of second harmonic and photoluminescence bands due to their distinct coherence properties. A detailed analysis yields fractional power dependence exponents in the range of 3–4, indicating the presence of multiple nonlinear processes. The range in measured exponents is attributed to differences in local field enhancement, which is supported by independent photoluminescence and structural measurements. Simulations based on Keldysh theory suggest contributions by three- and four-photon absorption as well as avalanche ionization in agreement with experimental findings.


Introduction
ZnO nanorods are high aspect ratio nanostructures, which are interesting for their salient linear and nonlinear optical properties and relative ease of sample preparation [1]. Grown for example on glass or silicon substrates, ZnO nanorods have been shown to be an effective medium for both second-harmonic generation (SHG) [2][3][4][5][6][7][8][9] and third-harmonic generation (THG) [7][8][9][10][11]. Their photoluminescence (PL) bands [12][13][14][15][16][17] in the ultraviolet (UV) and visible (VIS) spectrum, due to near bandedge excitonic and deeper defect level emission, respectively, are accessible even at room temperature (RT), in the former case due to a high exciton binding energy of E x = 60 meV [18]. Furthermore, the PL bands are interesting for applications, e.g. RT UV-lasing on a chip level [19] or photodynamic therapy [20]. Nanorod structures can exhibit relatively large conversion efficiencies for nonlinear optical processes. These efficiencies often exceed the obtainable conversion in bulk films of equal or larger thickness. For example, [1] reported efficient SHG for ZnO nanorod samples. While the coherent SHG process is relatively well understood there is a concomitant photoluminescence contribution deeper in the ultraviolet. This contribution is induced by multiphoton absorption of near infra-red (NIR) pump pulses, and its origin has been explained by a multitude of different nonlinear optical processes, ranging from two-photon [1, 4-6, 21, 22] over three-photon [3,23] to four-photon [24] absorption processes (#PA). As a further complication, power dependence measurements often reveal fractional exponents, which indicate the presence of at least two simultaneous nonlinear processes of different order. Moreover, it is puzzling that exponents corresponding to 2PA Journal of Physics D: Applied Physics

Field enhancement of multiphoton induced luminescence processes in ZnO nanorods
processes are sometimes observed, even though the pump photon energy was found insufficient to excite carriers across the bandgap E g = 3.37 eV in a 2PA process. As a possible explanation for these findings, various mechanisms such as band gap renormalization (BGR) [21], bandgap shrinkage due to heat accumulation [4,21], or Rabi oscillations [5,21] have been proposed.
Here, we report a method for reliable extraction of the power dependence components from interferometric frequency-resolved optical gating (iFROG) measurements [25][26][27][28]. The method is largely immune to power fluctuations and includes intrinsic consistency checks that are based on simultaneous measurement of the coherent SHG contribution and the multiphoton-induced PL signal. Moreover, iFROG allows the characterization of our pulses in situ by employing the nanorod samples as the nonlinear medium for THG, allowing us to precisely estimate the peak field intensity in the samples. Using ZnO nanorod samples that were manufactured with three different methods, exponents ranging from 3 to above 4 are found, indicating simultaneous three-photon and four-photon excitation, or possibly other nonlinear processes, with varying degrees of each contribution. Tunneling ionization is excluded based on a Keldysh parameter γ ≈ 0.04 1 for our experimental conditions, clearly indicating operation in the multiphoton absorption regime [29]. Independently determining the luminous efficiency of the three different samples together with morphological and structural characterization of the nanorod films, we are able to reconstruct the local field enhancement [30]. We observe a clear correlation between strong field enhancement and large exponents, which are supported by numerical simulations based on Keldysh theory in the multiphoton absorption regime [31,32]. Furthermore, the Keldysh simulations also suggest that avalanche ionization contributes significantly to the generation of carriers, allowing a stronger power dependence than could be expected from the efficiency of multiphoton absorption processes alone. Our results suggest that large enhancement factors can have a beneficial influence in obtaining the highest conversion efficiencies near the damage threshold of the films. Moreover, we find that different growth techniques can produce nanorod films with significant differences in SHG conversion efficiencies approaching one order of magnitude. We cannot, however, find any evidence for a prominent role of 2PA, or 2PA enhancing mechanisms in the UV-PL process. Instead, our experimental findings are fully explainable within the framework of Keldysh theory, as evinced by simulations.

Sample preparation and characterization
Three differently grown samples of ZnO nanorods were used in our experiments to ensure that our findings are not simply a feature associated with a single growth technique. There exist numerous growth techniques to produce ZnO nanorods, with varying results in achieved optical properties [33]. Our chosen techniques produce high quality nanostructures, as evidenced by the strong UV-PL emission spectra found in all of the selected samples. Some of the nanorods were grown using a vapor-liquid-solid (VLS) technique [34,35] with (i) Au nanoparticles or (ii) an Au thin film as a catalyst on substrates of fused quartz. Alternatively, (iii), we used vapor phase transport (VPT) on a ZnO buffer layer with no catalyst to grow ZnO nanorod films [36,37] on a fused silica substrate. Finally, a chemical bath deposition [36] sample was also tested but ultimately discarded as no UV-PL emission was observed. Significantly weaker (often by 2-3 orders of magnitude) PL emission from chemical bath grown ZnO samples compared to samples grown by high temperature processes is a very common observation [38]. Figure 1 shows scanning electron microscope (SEM) images of the chosen nanorod samples. Although the samples were grown using different techniques, the intrinsic material properties of all the samples are essentially identical according to previous studies [35,37]. To verify this, x-ray diffraction measurements were conducted (not shown), corroborating the high quality ZnO wurtzite structure with similar lattice constants in the samples, all close to the value obtained for bulk ZnO. Moreover, the Au catalyst used for samples (i) and (ii) was not detected in these measurements or any other of our characterization measurements. Whatever undetected Au remains in the vicinity of the nanorods, it is not expected to significantly affect the relevant properties of the ZnO because the Au resides almost entirely on top of the nanorods, as shown in previous studies [35]. Furthermore, the unavoidable Au impurity concentrations remaining within the nanorods do not significantly affect or degrade either the low temperature, or room temperature, UV-PL emission originating from excitonic recombination  [39]. Rather, the essential difference between the samples lies within the geometry of the nanorods, and the possible gold deposits on top of the rods, as these can affect the local electric field strength within the nanorods [40]. Nevertheless, due to the higher initial Au content in the thin film catalyst used for sample (ii) in comparison to the Au nanoparticles used for sample (i), it seems probable that more residual Au remains on the substrate, or in the vicinity or on the walls of the nanorods of sample (ii), leading to slightly different plasmonic properties between the two Au catalyst samples.
The VLS samples show similar nanorod diameters of ≈90 nm, and layer thicknesses of ≈1 μm, but while the ZnO nanostructures grown via Au nanoparticles appear to be randomly oriented in sample (i), the sample grown with an Au thin film (ii) shows ≈1 μm wide patches of similar orientation apparent in larger scale images (not shown). The nanorods in (iii) are larger, ≈170 nm in diameter with a layer thickness of ≈40 μm, and are oriented consistently along the substrate normal. For the subsequent experiments, where the samples are illuminated with NIR pulses at (or close to) normal incidence, the in-plane orientation of the nanorods with respect to the pump polarization can be considered as random because the spot size of the beam (≈5 μm) exceeds the transverse width of any observed features on the samples. This also explains why sample rotation with respect to the pump polarization was not found to significantly affect the measurements. The samples are visibly different in appearance. The VLS samples have a number of regions of variable translucence, most probably due to variance in the thickness of the nanorod layer. Sample (i) has a white and sample (ii) a pink hue with a dark ring, while sample (iii) is white, opaque, and uniform. We note that the observable pink hue of sample (ii) is indicative of plasmonic resonance effects visible in this sample due to the larger amount of gold present, mentioned further below, due to Au nanoparticles on top of the nanorods or elsewhere. Considerable fluctuation of the nonlinear response and the susceptibility to damage was observed in the visibly different regions of the samples, most prominently in sample (ii). A single region of each sample was therefore consistently selected for the power dependence measurements.
We further analyzed the SEM micrographs in figure 1 to obtain an estimate for the 2D fill factors, which range from 3.5% for sample (i) to 9% in sample (iii). Using these numbers and accounting for the film thickness, we compute the relative content of ZnO in the focal volume of the laser employed in the experiments discussed below. These estimated volume fill factors are listed in table 1 together with all other relevant parameters of the samples.

iFROG measurements
In order estimate the number of photons involved in the excitation of the UV-PL process, we adapted the pulse characterization technique iFROG for spectrally resolved power dependence measurements capable of distinguishing coherent and incoherent emission. The iFROG employs a pair of collinearly propagating, variably delayed but otherwise identical copies of an input pulse, and measures the spectrum of a nonlinear mixing process (typically SHG or THG) as a function of the time delay, producing a 2D spectrogram known as an iFROG trace. Two differing examples of iFROG traces are presented here. The first example, shown in figure 3(b), contains both SHG and UV-PL signals and is used for measuring the power dependence. A second example employs THG and is used to characterize our pulses, see figure 4(a). A more detailed explanation of iFROG is given in [25,27,28]. Further details of the underlying physics of the SHG and PL processes in ZnO nanostructures can be found e.g. in [16,41].
In order to detect possible 2PA enhancing mechanisms previously suggested [4,5,21], we aim to operate with photon energies below, but as close as possible to the 2PA threshold. To this end, an RG780 glass long pass filter was used to limit photon energies of the incident field to E < (E g − E x )/2 or >750 nm such that no direct 2PA can take place, to free excitonic states or between the valence and conduction bands (CB). Furthermore, the spectral overlap between the SHG at ≈400 nm and UV-PL at ≈380 nm is greatly reduced due to the suppression of the short wavelength tail of the SHG spectrum. The filtered fundamental spectrum spans the 750-970 nm range with the maximum at 810 nm; see figure 2. This close proximity to the 2PA threshold could not have been utilized with the time-resolved PL measurements described below, as the SHG and UV-PL contributions would have been inseparable. We proceed to explain why this is not the case with iFROG.
iFROG has been previously used by Schmidt et al [41] to measure the VIS-PL defect emission alongside SHG for a ZnO thin film. In particular, it was shown that the temporal sampling of iFROG readily allows identification of nonlinear processes with different coherence properties. At increasing time delays, the coherent harmonic generation, e.g. of SHG, results in an increasing spectral tilt of the interferometric fringes with respect to the central fringe at zero delay, which is parallel to the wavelength axis. This tilt effect can be clearly seen in the wavelength range >390 nm in figure 3. In contrast, incoherent spontaneous emission (PL) gives rise to a markedly different fringe pattern, in which all of the fringes are parallel, see the spectral range below 390 nm in figure 3. The ability to discern and isolate these two processes is particularly important under our experimental conditions where the UV-PL and SHG spectra can partially overlap. Previous accounts have resorted to identifying the spectral components, e.g. by assigning multiple Gaussian distributions to the recorded spectra, and then extracted power dependences from these curves [1]. This somewhat arbitrary assignment can be avoided with our method. An influence of the spectral wings from the neighboring spectral component is mitigated by applying Gaussian filters centered at the maximum emission wavelength of each component. Aside from its spectral resolving capability, the collinear geometry of iFROG allows tight focusing, which, in turn, enables us to reach high peak intensities close to the damage threshold of the material, which proved essential for achieving efficient UV-PL. As iFROG is essentially a spectrally resolved interferometric autocorrelation (iAC, [42]), a simple spectral integration of the trace provides the iAC signal. The contrast ratio r between the iAC signal at infinite delay and at zero delay can then be used for measuring the power dependence of the nonlinear process involved [41]. The number of photons involved m is connected to the ratio r via In the case of SHG, a two-photon process, one expects a contrast ratio of 1:8. For an underlying THG process, in contrast, a 1:32 contrast ratio is expected. Minor experimental imperfections can, however, easily induce deviations from these ideal values. In our measurements, we can readily obtain contrast ratios above 1:7 in the SHG component of the iFROG trace, yet observe ratios well above 1:32 in the UV-PL signal, which is clearly indicative of quite a high number of photons involved in the excitation of the associated UV-PL process; see figures 3(a) and (c). In addition to SHG and UV-PL, THG was also measured using the same experimental conditions except for the detection wavelength range. As the THG emission band is spectrally completely isolated from any other band, we can use it to characterize our pulses in situ using iFROG [11]. An iterative algorithm [43] was used to retrieve the complex electric field of the pulse, giving us a full width at half maximum pulse length of 39 fs; see figure 4. The power dependence exponents derived for the THG measurements were consistently in the range [2.9, 3.0], i.e. very close to the ideal value of 3 for a third order process, further corroborating the good performance of the measurement system.
Knife edge measurements gave a beam waist parameter of w 0 = 2.6 μm, setting the peak intensity at focus to I 0 = 0.5 TW cm −2 for the filtered pulse, which is an order of magnitude below the reported damage threshold of 5TW cm −2 for bulk ZnO [44]. Nevertheless, we observed optical damage in all of the nanorod samples in our experiments. We attribute this to local field enhancement and to heat accumulation. Moreover, a redshift of ≈50 meV of the UV-PL emission peak was observed for sample (iii) when subjected to a pump beam power ≈20% below the empirical damage threshold. This redshift is attributed to bandgap narrowing due to heat accumulation in the nanorods. The observed ≈50 meV shift corresponds to a ≈100 K increase in sample temperature [45]. To avoid damage and heat accumulation, we reduced the peak intensity in the focus by placing a neutral density filter (ND 0.6) in the beam path when necessary. The optical density of this filter was chosen such that no catastrophic damage or emission redshift was observed anymore.
In our experiments, we observed large fluctuations of up to an order of magnitude of the nonlinear emission intensity with respect to transverse position on the nanorod samples. In addition to variance in the amount of ZnO in the focal volume, this finding suggests a variance in the achieved local field intensities within the nanorods. This fact also manifests itself in the measured exponents for the UV-PL process, as we will discuss later on.

Photoluminescence measurements
We further characterized our ZnO nanorod samples by time-resolved photoluminescence (TRPL) measurements using a streak camera and multiphoton excitation. A Ti:sapphire oscillator was tuned to a center wavelength of 830 nm in order to have the least possible amount of overlap between the SHG and UV-PL emissions while still maintaining similar exper imental conditions to the iFROG experiments discussed above. The spectral separation here could not be made as small as with the iFROG measurements, because (unlike with iFROG) a possible overlap of PL and SHG would be difficult to detect with the streak camera measurements. Pulses with a duration of ≈ 80 fs and 570 mW average power were focused onto the sample at approximately 20 • angle of incidence, and the emission spectrum was recorded by a combination of a spectrograph and a streak camera synchronized to the oscillator repetition rate. These measurements were repeated for samples (i)-(iii) and, additionally, with a 1 mm thick single crystal ZnO sample (figure 5).
Two spectrograph gratings with different groove densities were used: a sparse grating (300-580 nm) for figure 5(a), capable of capturing the entire emission spectrum, and a dense grating (340-400 nm) that captures only UV-PL, allowing exposure adjustment independent of SHG emission, for figure 5(b). The former measurements were spectrally integrated for a comparison of the luminous efficiencies of the samples. For the latter setup, the higher spectral resolution decreased signal levels such that it was necessary to vary the incident intensity to ensure adequate signal-to-noise ratio for each sample. Thus, the luminous efficiencies are not comparable in figure 5(b). All samples including the bulk sample show a PL maximum in the range from 382 to 391 nm, but    . Similar fast decay times have previously been attributed to plasmonic effects due to Au nanoparticles in close proximity to ZnO nanorods [46]. This observation corroborates the suggestion above that, in comparison to sample (i), a larger amount of Au nanoparticles is present in sample (ii), which the clearly visually observable pink hue of the sample also indicates. Separating the SHG response at 415 nm, we used these spectra to compute the luminous efficiency η relative to the bulk material, see table 1. The differences between the nanorod samples and bulk may not appear overly impressive at first sight. However, one has to consider that, e.g. sample (ii) only hosts ρ = 0.07% of the ZnO in the bulk sample within the focal volume of the laser. This allows us to estimate the field enhancement factor [30] γ from the structural and PL data according to where m is the number of photons involved in the multiphoton excitation process, obtained in turn from the iFROG measurements.
A third set of time-resolved PL measurements was conducted to directly measure the power dependence of the PL emission, i.e. the order m of the multiphoton excitation process. The incident pulse was attenuated via neutral density filters with optical densities ranging from 0.1 to 1.0, allowing almost an order of magnitude variation in incident average power. The applicable power range was limited from above by the onset of material damage, and from below by a decreasing signalto-noise ratio. The isolated PL emission was integrated over time and wavelength to obtain a single value for the luminous efficiency as a function of the incident power. A linear fitting procedure for the log-log presentation of the measured data sets yields the values m, presented in the bottom of table 2.

Experimental results
The top part of table 2 shows the exponents m with 95% confidence intervals (assuming normal distributions) for the three nanorod samples, measured with iFROG. all of the samples display fractional UV-PL exponents between 3 and 4. The ideal value of m SHG = 2 is within error margin for samples (i) and (ii), but an absolute deviation of ≈0.1 is found for sample (iii). This observation is explained by the spectral proximity of UV-PL and SHG emission maxima for the latter sample. As the two emission spectra partially overlap there is resulting interference with the exponent measurement.
The values obtained with time-resolved PL measurements (table 2, bottom) are in good agreement with the iFROG measurements, with the results from the two techniques overlapping within the computed error margins. The most important outcome here is that the relative ordering of the exponents m PL for the different samples is the same for the two different measurement techniques, i.e. highest for sample (ii), lowest for sample (iii) with sample (i) in the middle, but found closer to sample (ii). The largest deviation from the value measured with iFROG is found for sample (ii), where the streak camera measurements yield a very high exponent of 4.3. This large value is accompanied with a very large error margin of 0.8 units, mainly due to the fairly weak luminous efficiency of sample (ii) compared to the other samples. As the exciting fields for each measurement differ significantly (iFROG: 40 fs pulse, broader spectrum down to 750 nm; time-resolved PL: >80 fs pulse, narrower spectrum down to 810 nm), the derived exponents m PL are not directly comparable. Avalanche ionization in particular can play a much more prominent role in exciting carriers for an increased pulse duration-a possible explanation for the high exponent derived for sample (ii) with the time-resolved PL measurements. Avalanche ionization becomes more efficient with an increasing conduction band occupation, and it may be the case that a significant population is only reached with sample (ii). Further analysis of the exponent results is conducted with the iFROG values, as their error margin is significantly lower.
Using the values m PL in table 2 for the iFROG measurements, the field enhancement factors γ in table 1 are computed using equation (2). Enhancement factors ranging from 1.5 to 2.5 are found, corresponding to 5-39 times higher SHG efficiency compared to bulk ZnO. The thinner VLS deposits with their finer needles clearly showed stronger field enhancement than the VPT samples, even when considering the conservatively estimated error bars; see figure 6. We believe this significant impact of nanorod structure on the nonlinear conversion efficiency can be of practical importance in application design. Even though it exhibits the lowest field enhancement γ, the strongest UV-PL emission was detected from sample (iii), followed by (i) and (ii). The order is reversed when considering the highest sustainable intensity with no significant material damage (of type 1 as discussed below), e.g. sample (iii) has the lowest damage threshold. This trend follows the 3D fill factors in table 1, such that the thinnest layer is least susceptible to optical damage. The thickest layer is most likely least efficient in diffusing heat into the substrate. This accumulation of heat leads to positive feedback, where the heat first decreases the band gap, which then leads to increased absorption, which in turn allows more heat to be deposited from the optical field, until excessive absorption causes a catastrophically high CB population and material damage. Consequently, samples (i) and (ii), which can sustain higher intensities without damage, also show the highest mean exponents at m PL ≈ 3.8, while sample (iii) gives a lower value of m PL ≈ 3.4, suggesting proportionally more 3PA with respect to 4PA or avalanche ioniz ation in sample (iii) than in the others. This finding is explained by the MPA simulations below (figure 6), showing that a higher field intensity leads to a higher observed power dependence exponent. The fact that sample (iii) is the strongest emitter regardless of its lowest γ of the three samples is explained simply by its highest ZnO concentration in the focal volume: a weaker local field can still result in a stronger emission if there is a larger volume of emitter material.
While the measurements for samples (i) and (ii) were conducted such that the pulses reach the nanorods before the glass substrate, sample (iii) was reversed so that the beam is focused through the substrate onto the nanorods. This reversal shifted the UV-PL emission wavelength to shorter wavelengths or, more likely, less of the short wavelength photons were reabsorbed and scattered by the nanorods before reaching the detector. Thus the overlap of UV-PL and SHG was reduced, making exponent extraction easier. This would also decrease the achieved peak intensity because of reflection losses from the substrate, contributing to the observed decrease in m PL .
As noted above, considerably higher values of power dependence exponents were measured for the Au catalyst samples (i) and (ii) in comparison to sample (iii), grown without a metal catalyst, with sample (i) showing the highest individual values. One plausible, albeit speculative, explanation for this difference is that the local field achieved in the Au catalyst samples is assisted by plasmonic effects in residual Au nanoparticles on top of the ZnO nanorods or elsewhere. Field enhancements due to a combination of plasmonic effects [40] and 1D nanorod morph ology [47] could significantly exceed the field enhancement solely due to the 1D morphology of the nanorods. Comparing the Au catalyst samples, the slightly higher power dependence of sample (ii) is likely connected to its higher Au nanoparticle content. In addition to a higher field enhancement, more Au nanoparticles can lead to enhanced bandedge emission via coupling of Au surface plasmons and ZnO excitons located close to the nanorod surface and hence to nearby Au nanoparticles [46]. Two different types of damage incurred by the nanorods were identified: (1) quickly occurring damage within <1 s, resulting in a change of the measured exponent, and (2) gradual degradation of the UV-PL emission intensity over 1 s with no apparent change in the exponent. The former damage mechanism occurs only at high peak intensities while the latter was observed even at lowest intensities that just suffice to obtain an observable UV-PL emission. The fast damage process is thought to occur due to excessively high CB population, leading to severe and abrupt material damage. The slow process, in contrast, is likely due to accumulation of defects over time. Most often type 1 damage causes an initial exponent of >3 to drop to ≈2, which is accompanied by a significant reduction of the UV-PL emission. For very weak input fields, there were rare observations of nanorod sites where m PL ≈ 2. These sites were converted to a m PL > 3 sites upon exposure to higher pump power. It appears, though, that the weak incident field is insufficient for notable 3PA to occur. The very weak remnant short wavelength tail of the pump causes 2PA on these scarce sites where exceptionally efficient spatial configuration of nanorods allows a detectable amount of 2PA to occur. For higher pump powers 3PA becomes efficient, and m PL > 3 are observed, likely accompanied by the destruction of the 2PA site. Further investigation is warranted to explore the causes underlying these changes in the power dependence of UV-PL in ZnO nanorods.

MPA simulation
We model the MPA process by numerically solving the differential equation for the CB electron population N(t) [32]: The purpose here is to estimate the stregth of the UV-PL emission-assumed to be proportional to the CB population-as a function of the incident intensity. In the above equation, the first term on the rhs describes avalanche ionization (AI) with the coefficient α. The second term accounts for MPA of integer orders m ranging from k to l with MPA coefficients β m , and weighting factors µ m . Finally, the third term includes relaxation effects appearing at a time constant τ. ξ is the ratio of the maximum internal and the incident intensity. This structural factor accounts for interference effects in the medium and adjusts the local field intensity. Mero et al used ξ = 0.66 for their thin film samples, while factors in the range ξ ∈ [0, 6] have been reported for ZnO nanorods [30]. For the lack of a better estimate for the impact ionization coefficient α for ZnO nanorods, we took the value for TiO 2 from [32], as the band gap of TiO 2 is almost equal to that of ZnO. The MPA coefficients β m for the mth order absorption process are given by Keldysh theory of photoionization [31,32]: Here, n 0 ≈ 1.96 is the linear refractive index of ZnO [48], and ω 0 is the carrier frequency, taken as the maximum of our separately measured spectrum. We have allowed the band gap energy E g to depend on the CB population N in order to simulate band gap narrowing due to BGR. For this purpose, we use values computed by Bányai and Koch [49], as reported by Versteegh et al [50]. The carrier reduced mass is taken to be m r = 0.28 m e [51], where m e is the electron rest mass. The carrier recombination rate for ZnO is typically described by an exponential decay with two time constants: a slow constant in the order of tens of ps accompanied by a fast initial decay of ≈1 ps [50,[52][53][54]. As the interaction between the sample and both of our pulses takes place within a <200 fs time frame, the fast time constant defines the decay of CB population generated by the first pulse, subsequently seen by the second pulse. We therefore refrained from more sophisticated modeling of the relaxation process in equation (3) and only employ a single τ = 0.7 ps time constant for the decay [50]. For the simulations, a sech 2 pulse intensity profile with the measured full width at half maximum length of 39 fs was used.
The weighting factors µ m account for the overlap of different harmonics of the fundamental spectrum according to the density-of-states (DoS) such that the availability of CB states affects the probability of a given order of MPA, e.g. the minuscule overlap between the SHG spectrum and the DoS leads to a large reduction of 2PA contributions. The µ m are computed as follows. Harmonics of the separately measured laser spectrum S m (ν) (after filtering) are integrated over the approximate DoS D(hν, N) for the CB of ZnO [55] according to The denominator causes normalization of the weighting factor µ 3 to unity. The DoS D(hν, N) is a function of N as BGR can shift the band gap, thus changing the overlap of the harmonics and the DoS, and subsequently the MPA probabilities. For the unperturbed band gap, we set µ m (N) → µ m (0) and obtain the weighting factors compiled in table 3, ranging from µ 2 (0) = 1.3 × 10 −5 to µ 5 (0) = 5.0. The former weighting factor clearly confirms a strong suppression of second-order effects due to the action of the longpass filter. Moreover, as the µ m increase with growing m, higher-order MPA coefficients, in particular of order 4, contribute to the photoionization process. Due to BGR, an increasing field intensity will shift the band gap by ∆E(N) due to larger induced CB populations. Such high carrier concentrations are, however, not expected to be maintained for periods much beyond the fast initial decay time T due to diffusion of carriers and nonradiative recombination at the substrate interface. As PL emission occurs over tens or hundreds of ps, a transient redshift would not be visible for our slow detector in the iFROG experiments, and not even for the streak camera with a time resolution of ≈40 ps. However, even though the emission spectrum is not affected by BGR, the effect may possibly increase the likelihood of 2PA as we will see further on. A shift of the band gap will also affect the MPA coefficients β m , as can be seen from equation (4), and the weighting factors µ m . When BGR is considered ( E g → E g (N)) these two variables become functions of N, and must be computed for each step of the simulation as N evolves.
For the simulation, we consider MPA orders 2 through 6. By including orders as high as 6 we consider the generation of hot carriers far above the CB edge that will thermalize before eventual UV-PL emission.

Discussion
Let us first employ our model to address the previously discussed damage mechanism of type 1. To this end, we calculate the intensity required to reach the value of N corresponding to a plasma frequency equal to the carrier frequency of the pump. This commonly accepted condition is expected to lead to an exponentially increasing linear absorption of the pump. For our experimental parameters, this limit is 2.1 TW cm −2 in bulk ZnO, corresponding to ≈1% partial ionization. However, for some of the ZnO samples we already observe damage at much lower intensities. We explain these findings by the intensity enhancement γ 2 inside the nanorods, which can amount to values up to about 6.2 (table 1), thus allowing the damage threshold to be exceeded.
We then utilized the model to theoretically estimate the resulting effective exponent of the MPA process m PL (I), see figure 6. Focusing our attention to the data corresponding to the iFROG experiments plotted in red, a dominance of third-order processes is seen at intensities in the range I < 10 11 W cm −2 . Further increasing the intensity, one then observes an increase of m PL towards 4, which is reached at I = 10 12 W cm −2 . For comparison, we also included measured exponents for the three nanorod samples in figure 6. The intensity coordinates for the experimental data were computed by multiplying the experimental peak intensity I 0 = 0.5 TW cm −2 by the relative intensity enhancement γ 2 for each sample (table 1). The incident intensities were further attenuated with the respective ND filters used for each sample. Next, we repeated the numerical simulations using parameters corresponding to the time-resolved PL measurements, the results are plotted in blue in figure 6. For this case, a higher slope for the exponent is obtained from the model, mainly due to an increased contribution from avalanche ionization brought by the longer pulse duration employed (90 fs). Again, the experimental data points are shown, with similar corrections made as for the iFROG experimental data. As can be seen from figure 6, correcting for the local field enhancement results in a good agreement between the measured and the simulated power dependence exponents for both iFROG and TRPL. Consequently, it appears that there is no need to invoke any 2PA enhancing mechanisms to explain the varying exponents in the samples. For example, inclusion of BGR merely increases the 2PA probability slightly for near-breakdown intensities. Taking a closer look at the simulated excitation process corresponding to the iFROG experiments (figure 7), avalanche ionization is found to play a significant role in the ionization of the nanorods for high intensities, as is to be expected for our relatively long pulse duration of ≈40 fs [32,56,57]. At the experimentally determined peak intensity of 0.5 TW cm −2 , our simulations indicate that AI accounts for approximately 30% of the generated CB electrons, while the remaining 70% are due to MPA. An equal amount of CB electrons are created by the two effects at an intensity of ≈1 TW cm −2 . The MPA process always precedes and serves to seed the AI, which reaches its maximum rate of generated CB electrons after a time delay of ≈12 fs. In comparison, for a shorter pulse of 7 fs, our model shows that AI does not exceed MPA before material damage is expected. The magnitude of AI is determined by the avalanche coefficient α of the nanorods, which we unfortunately cannot easily measure or find in the literature. For the case of a DC electric field, an ionization coefficient can be found for bulk ZnO [58]. However, we find that the interaction between the rapidly oscillating, strong electric field of the ultrashort pulses and the nanostructures is better described by the dielectric breakdown model for fs pulses in [32], and chose to use values therein.
The relative contribution of each MPA order for the total modeled excitation of ZnO is also presented in figure 7. If we restrict ourselves to the experimentally relevant intensity range from I = 10 11 to 10 12 W cm −2 , it can be readily seen that 3PA dominates, while 4PA makes an order of magnitude weaker contribution. The 3PA and 4PA contributions reach parity at I ≈ 7 × 10 12 W cm −2 , which is well beyond our experimentally accessible range, yet the fractional exponents in the range from 3 to 4 are clearly explained, especially when AI is considered. 2PA, in contrast, does not contribute to a significant population of carriers in the CB, even when BGR is considered.
Repeating the simulations with a longer pulse length corresponding to the time-resolved PL measurements, the relative contribution of AI in charge carrier generation is further increased. This allows exponents greater than 4 to be reached, even for the experimentally obtainable incident intensities, lending an explanation for the high value of 4.3 measured for sample (ii); see the bottom of table 2. This finding is a further evidence towards the significance of AI in the UV-PL process for ZnO nanorods.

Conclusion
We have investigated the UV-PL emission from ZnO nanorods induced by MPA of NIR fs pulses using a novel experimental approach based on a combination of the ultrashort pulse characterization method of iFROG, simultaneously measuring laser harmonic generation and multiphoton induced luminescence with spectral and temporal resolution, with suitable filtering to ensure a minimum spectral overlap of SHG and UV-PL emissions and suppressed direct 2PA contributions. Moreover, iFROG allowed us to characterize the pulses in situ using the same nanorod samples. The acquired pulse intensity profile was used to better estimate the peak incident field intensity within the samples, a crucial quantity for the UV-PL experiments. The UV-PL emission of three differently grown ZnO nanorod samples showed different intensity dependences, with exponents varying from 3 to above 4. Independently conducted time-resolved PL measurements yielded similar exponent values, and the results were found to agree within the error margins. The varying exponents are explained by a differing field enhancements γ depending on the growth method. This enhancement seems to be most effective in the VLS samples, where Au was used as a catalyst, reaching a 2.5 fold increase. Considerably weaker enhancement of γ = 1.5 is found for the VPT sample grown without a metal catalyst. Plasmonic effects in residual Au nanoparticles at the top of the VLS nanorods or elsewhere are identified as a possible explanation for this trend. The field enhancement in each sample was estimated by the quantum efficiencies computed from SEM and TRPL measurements. A seemingly modest difference of one unit between the factors γ for the electric field enhancement for different samples becomes significant when nonlinear processes are considered. SHG efficiency, for example, is increased by a factor of between 5 and 39, depending on the growth technique, i.e. the growth technique can change the SHG efficiency of the nanorods by a factor of 8. We believe that this difference in conversion efficiency as a result of different growth techniques can be of practical importance in application design. after which AI dominates. The relative contribution of each order of MPA is plotted using solid lines. 3PA is dominant in the intensity region of ≈10 11 W cm −2 corresponding to experimental conditions. Introduction of a band gap reduction due to BGR for high N raises the 2PA contribution for high intensities I > 0.4 TW cm −2 close to the damage threshold, but only to a level between 5PA and 6PA.
The computed damage threshold intensity I crit is also shown (dotted magenta line).
Optical damage to the nanorods was found to affect the power dependence and intensity of the UV-PL emission. Moreover, the catastrophic optical damage threshold could be reached for all the investigated samples, with the Au catalyst VLS samples able to sustain a significantly higher incident intensity than the thicker VPT sample. This is explained through greater heat accumulation in the thicker VPT sample, consistent with an observed bandedge redshift and a lower power dependence than was obtained for the VLS samples. For subsequent experiments, the incident intensity was restricted such that no significant heating, indicated by a redshift, was observed in any of the samples. The catastrophic optical damage is generally explainable by the onset of plasma-induced linear absorption in the samples, but we note that we also found evidence for slower degradation mechanisms. Finally, we observed high-intensity induced noncatastrophic changes of the exponent of the photoemission.
In all the samples investigated, the experimental data and numerical simulations, when taken together, strongly suggest a competition between 3PA, 4PA and avalanche ionization effects, which, in combination, give rise to the fractional exponents observed in the experiments. The effects are explainable within the framework of Keldysh theory, and there appears to be no need to invoke 2PA enhancing mechanisms such as BGR. While the latter certainly increases 2PA, these contrib utions do not appear strongly enough enhanced to significantly influence the measured exponents of the UV-PL emission. Moreover, for the same reasons, we could not find any indications for Rabi oscillations acting to enhance 2PA. Furthermore, we observed significant transverse variation in the measured exponents, which is indicative of local variations of the field enhancement due to spatial variations in the nanostructures and Au deposits. In light of these observations, we conclude against any significant contribution of 2PA to the UV-PL emission, provided that a significant portion of the fundamental field does not correspond to photon energies exceeding half the band gap energy.