4.1 Classical Model
The results are first analyzed with the Fresnel equations. The reflectivity R with Fresnel loss is calculated as follows[19]:
$$R={\left\frac{N{cos}\left(\theta \right)+\sqrt{1\frac{{{sin}}^{2}\left(\theta \right)}{{N}^{2}}}}{N{cos}\left(\theta \right)+\sqrt{1\frac{{{sin}}^{2}\left(\theta \right)}{{N}^{2}}}}\right}^{2}$$
1
where θ is the angel of incidence, N is the complex refractive index \(N=n+\text{i}\kappa\), \(n\) is the refractive index and \(\kappa\) is the extinction coefficient. In this experiment, the reflectivity is measured at θ = 45º.
The imaginary part of the complex refractive index is calculated from the absorption at normal incidence[19]:
$$\kappa =\frac{{\lambda }_{0}\text{log}\left(T/(1R)\right)}{4\pi d}$$
2
where d is the thickness of the sample, \({\lambda }_{0}\) is the wavelength of the laser in vacuum, R is the reflectivity and T is the transmissivity.
The complex refractive index is numerically calculated from the solid fitted curves in Fig. 2 with equations 1 and 2. The real and imaginary parts of the refractive index are shown in Fig. 3. The calculated imaginary refractive index at low energies is 4.9, which agrees the results reported by Johnson et al. in 1972[20].
Due to the scarce amount of data at 0–10 kJ/m2, different fittings of the transmissivity were tested in Fig. 4 (a) shown in different colors, and the corresponding real part is shown in Fig. 4 (b). The results show that Re(n) is robust to different estimations of the transmissivity below 10 kJ/m2 as the curves overlap at intensities below this value, while the maximum of the calculated Re(n) could vary from 1 to 5 if a lower transmissivity was chosen, and its position could shift from 10 to 30 kJ/m2.
The maximum peak around 10 ~ 30 kJ/m2 from classical models suggests that there are higher orders of nonlinearity at higher energies and needs further data, but the overlapping regions below 10 kJ/m2 allow us to calculate the second order nonlinearity below the damage threshold, and the nonlinear refractive index is then estimated to be:
$${n}_{2}=\frac{0.7}{3.1\times {10}^{16}W/{m}^{2}}=2.25\times {10}^{17}{m}^{2}/W$$
3
This is consistent throughout different estimations below about 1kJ/m2 in Fig. 4 (a). This value is about one order of magnitude smaller than the upper bound of \(1.2\times {10}^{16}{m}^{2}/W\) provided by Goetz et. al., measured with a 130 fs, 800 nm laser below the damage threshold[13].
The real part of the nonlinear susceptibility could be calculated from the nonlinear refractive index:
$${\chi }_{Re}^{\left(3\right)}=2{ϵ}_{0}{n}_{0}^{2}c{n}_{2}=3.87\times {10}^{21}{m}^{2}/{V}^{2}$$
4
And similarly, the nonlinear absorption, taking the region below 1kJ/m2 is calculated to be
$$\beta =\frac{\left(5.37.82\right)\times {10}^{7}}{3.1\times {10}^{16}}=8.15\times {10}^{10}m/W$$
5
And the corresponding imaginary part of the nonlinear susceptibility is
$${\chi }_{Im}^{\left(3\right)}=\frac{{ϵ}_{0}{n}_{0}^{2}c\beta \lambda }{2{\pi }}=6.64\times {10}^{21}{m}^{2}/{V}^{2}$$
6
Which is about two orders of magnitude smaller than previously reported results, as Rotenberg et al.[21] reported to be \((76.8+4.3\text{i}) \times {10}^{19}{m}^{2}/{V}^{2}\) using a zscan method and Renger et al. [8] measured a real part of \(2.0\times {10}^{19}{m}^{2}/{V}^{2}\) with a fourwave mixing method, while Smith et al[7] suggested that the thirdorder susceptibility could be as high as \((3.7+5\text{i}) \times {10}^{14}{m}^{2}/{V}^{2}\). Despite the descrepancy, this result agrees with the previously mentioned decreasing trend[5]. Boyd also concluded that the real part of the nonllinear suscptibility spans a wide range and even changes sign due to the low transmission of gold films. Although the thick 100 nm gold film exhibits the same low transmission, the real nonlinear index is positive at 10 kJ/m2 level, as shown by the positive slopes in Fig. 4 (b).
Although Re(n) is insensitive to transmission values below 10 kJ/m2, this behavior changes above the damage threshold where transmissivity can reach several times its original value depending on the estimated transmission. However, it cannot be arbitrarily large because as the transmissivity gets lower, a classical solution for Re(n) will not exist.
This increase of the real refractive index is supported by the morphology of laser irradiated samples. Wang et al. [22] reported that certain surface structures would appear when samples are irradiated around the damage threshold, and the size of the structures is related to the interference or diffraction of waves inside the metal within the skin depth[22]. A reduced size resembles a reduced light wavelength inside the sample, and thus a real part refractive index greater than one. Similar experiments on gold were performed, and the results are shown in Fig. 5 where a laser slightly more intense than the damage threshold was focused onto a gold surface. Ripples are found at ridges and edges with a size about 1/4 of the original wavelength, which indicates a real refractive index of about 4 in this case and supports a positive second order nonlinearity at this energy level.
4.2 PlasmonPhoton Exchange (PPE) Model
A classical model can be used to estimate the indices below the damaging threshold, but it is unable to explain the optical constants at higher fluence. Although higher orders of permittivity could be introduced, a model studying the process is needed to help guide future research. For gold films with a few tens of nanometers thick, when irradiated by a nearinfrared ultrafast laser, the optical behavior is dominated by the plasmonic intraband absorption, which is caused by free carriers in gold’s s and p bands[23]. The plasmonic properties of gold are closely connected to its optical nonlinearities,[5] and a plasmonic model was chosen to explain the process. Despite the sample being damaged by the laser pulse, it takes picoseconds for the hot electrons to transfer the energy to the lattice to cause the damage depending on laser intensity[14], while the plasmonic response time is about 6 fs for gold[24]. This leaves us a window to study the electron excitations before the lattice is involved.
To set up a model, the energy flow of the system was considered. As shown in Fig. 6, energy transferred to induced plasmons and absorbed photons with a ratio \(p\), plasmons dephase into heat through damping \({\Gamma }\), photons exit the system and are received by the detector through an emission term \(E\). Plasmons and photons could also turn into each other, \({T}_{ph}\) is the transition rate for a plasmon to decay into a photon, and \({T}_{pl}\) is the transition rate from photons to plasmons.
The energy flow allows us to evaluate the plasmon number \({n}_{pl}\) and photon number \({n}_{ph}\) as a function of time, by solving the two following differential equations:
\(\frac{d{n}_{pl}}{dt}=pI\left(t\right){\Gamma }{n}_{pl}{T}_{ph} \left({n}_{ph}+1\right){n}_{pl}+{T}_{pl}\left({n}_{pl}+1\right){n}_{ph}\)

(7.1)

\(\frac{d{n}_{ph}}{dt}=(1p)I\left(t\right)E{n}_{ph}+{T}_{ph}\left({n}_{ph}+1\right){n}_{pl}{T}_{pl}\left({n}_{pl}+1\right){n}_{ph}\)

(7.2)

where \(I\left(t\right)\) stands for the total quasiparticles induced by the laser and is proportional to the laser intensity. For photon to plasmon transitions, the \(\left({n}_{pl}+1\right)\) factor in the term \({T}_{pl}\left({n}_{pl}+1\right)n\) is proposed by Finazzi et al. based on a secondquantization method [25], and is a result of stimulated emission whereby existing photon occupation increases the transition probability. Similarly, existing plasmons increase the transition probability from plasmons to photons. Although this secondquantization term is derived for nanoparticles, it is valid for 100nm thick gold films as the optical behavior at 800 nm is dominated by the intraband transition while the interband transition wavelength at 520 nm is determined by the band structure of gold[26]. Although the interband transition wavelength of gold can decrease with size[27] because of a quantum confinement when the size is smaller than the Fermi wavelength[28], which is not obvious when the size is larger than a few of nanometers[29].
The simulation with equations 7.1 and 7.2 shows that the intensity dependence of the predicted transmissivity is a direct result of these nonlinear terms, without which the intensity increases plasmon and photon numbers proportionally and does not influence the calculated transmission as explained in the supplementary data.
Once the equation is solved, the transmission of such a system is predicted by dividing the total number of emitted photons by total incoming photons, which can be compared with experimental data.
$$Transmission=\frac{E\int {n}_{ph}\left(t\right)dt}{\int I\left(t\right)dt}$$
8
In the simulation, most of the variables are calculated instead of fitted. The \({\Gamma }\) damping factor is determined from plasmon decoherence time, which is about 0.18 fs− 1[24, 30, 31]. The emission factor E was numerically solved to be 8.43×10− 5 fs− 1 to match the lowenergy transmissivity of about 2.4×10− 3, which may depend on the thickness and wavelength. The simulation results are insensitive to the energy partition parameter \(p\) at the range of 0.1–0.9, and 0.8 was chosen for better convergence.
For the energy exchange terms \({T}_{pl}\) and \({T}_{ph}\), Finazzi et al. proposed a calculation from second quantization [25]:
where (κ, α ) is the wave vector and polarization of the photon, \(w\) is the wave vector of the plasmon, \({\omega }_{w}\) and \({\omega }_{k}\) are the angular frequency of the plasmon and photon, \(V\) is the volume of a box containing the system.
The summation is taken over all possible frequencies, polarizations and wave vectors, and the two expressions should be identical. However, this summation is not valid for an arbitrary photon state when a photon excites a plasmon. The equations assume a dipole approximation and neglect interband transitions at high energies, and thus a cutoff frequency must be introduced when summing over all possible photon states, which leads to the small deficient of \({T}_{pl}\) from \({T}_{ph}\). Therefore, we introduce a new variable \(\delta\) to account for the difference \(\delta ={T}_{ph}{T}_{pl}\). The transition rates themselves are estimated from experiments to be around 0.14 fs− 1 [25]. The laser term \(I\left(t\right)\) is in the form of a delayed gaussian pulse[14]:
$$I\left(t\right)=\frac{{I}_{0}}{\sqrt{\pi }w}{e}^{{\left(\frac{t2w}{w}\right)}^{2}}$$
10
Where \(w\) is the pulse duration of the laser pulse and \({I}_{0}\) is the total number of photons in the laser pulse, each 800 nm photon will contain \(\frac{hc}{\lambda }=2.5\times {10}^{19}J\) of energy. For 6 mJ pulses, 2.4×1016 photons were calculated. Using this value assumes the entire area of irradiation is as a single system, which is a cylinder with 50 \(\mu m\) radius and 100 nm thickness. This radius is smaller than the range of surface plasmon polaritons[32] that could span millimeters. Our simulation shows that if the system is divided into two independent subsystems and each has halved total photons, doubling the difference \(\delta\) would give the same result. Thus, the effective size of the subsystems could not be determined from our model, but it gives a numerical correlation between the plasmonic coherent range in the process, which could also be used to determine the cutoff frequency when summing Eq. 9.
With the above conditions, and \(\delta\) fitted to be 7.09×10− 16 fs− 1, the simulated transmission matches the trend of experimental data, as shown in Fig. 7. The solid line in Fig. 7 presents a loosely linear increase of transmissivity with respect to energy, while the nonlinearity introduced this way is less significant at lower intensities. For simulations in Fig. 7, the entire area irradiated by laser was treated as a single system, and all plasmons and photons contribute to the nonlinear terms; however, this coherency might not span 50 µm and the actual system is much smaller. In this case, the surface is subdivided into independent patches, with the pulse energy evenly divided. The result would lead to a much lower total incident photons, and the estimated \(\delta\) that fits the experimental data is approximately inversely proportional to the total photons. The \(\delta\) is estimated to be ~ 10− 7 fs− 1 if a spot with 5 nm radius was chosen as the size of an independent subsystem.
For the case of 1016photon 35fs laser pulse, number of plasmon and photons as a function of time is plotted in Fig. 8.
The model also suggests that the total time of this event lasts about 500 fs, which agrees with the electron scattering time scale[33, 34]. The maximum photon number is much greater than the plasmon number, which is a result of the different damping rates. Consider an equilibrium state with a constant photon input and no exchange terms, the stable plasmon and photon numbers in Eq. 7 become \({pI{\Gamma }n}_{pl}=0\) and \(\left(1p\right)IE{n}_{ph}=0\) respectively. Since \(E\) is much smaller than \({\Gamma }\), the photon number in this steady state is much greater than the plasmon number. By adding the exchange terms back, they become coupled quadratic equations, but a similar conclusion could be reached (see appendix A3).
Going back to a pulsed input, the transition coefficients in the exchange terms are large when plasmon and photon number are large, they balance the photons and plasmons quickly, and each instant could be viewed as a steady state. Thus, the model suggests that most of the energy is in the photon state.
While photon numbers decrease quasilinearly after its maximum, plasmon numbers tend to remain constant before photon number drops to some similar order of magnitude. This is again the result of the previously discussed equilibrium state. After the laser pulse is fully absorbed, ohmic damping and emission will dissipate the energy together. The emission of photons is much slower than the ohmic damping of plasmons, and the system is in a photonrich state compared to equilibrium, and the exchange terms transform excessive photons into plasmons. The photon and plasmon numbers decrease as the equilibrium state would predict, until the event ends. The emitted light is proportional to the total photon number, and this stretched pulse could potentially be verified by an autocorrelationprocess.
This simulation can also be used to explain the transmissivity at nondamaging intensities. Other research suggests a positive nonlinear absorption coefficient \(\beta\), as summarized by Boyd et al. [5]. This positive \(\beta\) stands for a decreased transmission as intensities increase. Since the laser energies differ by several orders of magnitude, the results could be credited to different regimes, and the enhanced absorption is assumed to be governed by an increased damping factor \({\Gamma }\). A temperature damping could explain the result qualitatively, as shown in the appendix A4. If this temperature damping term is considered in the high energy case, increasing the difference \(\delta\) by ~ 15 times would give a similar result.
However, this simulated result did not show a decreasing trend at lower energies, as a positive nonlinear \(\beta\) would do. This could be a result of pulse duration. If a 500fs pulse was used, as shown in Fig. 7 with dotted line, the transmission would decrease at lower energies before the plasmon interaction term dominates at higher energies. In this model, the low energy behavior is dominated by electron temperatures, which reaches the cutoff temperature rapidly (see appendix A4) at higher fluences and is insensitive to the pulse duration. The plasmon interaction term, however, is sensitive to cumulated plasmons, and an extended input pulse provides enough time for plasmons to exit the system through damping. This could also be used to explain why classical nonlinear coefficients decrease as pulse duration decreases. The decreasing trend from the electron temperature damping term is insensitive to the pulse duration, but the increasing trend from the second quantization term dominates when a shorter pulse introduces more plasmons and photons to the system rapidly.