Ultrafast nonlinear dynamics of surface plasmon polaritons in gold nanowires due to the intrinsic nonlinearity of metals

Starting from first principles, we theoretically model the nonlinear temporal dynamics of gold-based plasmonic devices resulting from the heating of their metallic components. At optical frequencies, the gold susceptibility is determined by the interband transitions around the X,L points in the first Brillouin zone, and thermo-modulational effects ensue from Fermi smearing of the electronic energy distribution in the conduction band. As a consequence of light-induced heating of the conduction electrons, the optical susceptibility becomes nonlinear. In this paper we describe, for the first time to our knowledge, the effects of the thermo-modulational nonlinearity of gold on the propagation of surface plasmon polaritons guided on gold nanowires. We introduce a novel nonlinear Schrödinger-like equation to describe pulse propagation in such nanowires, and we predict the appearance of an intense spectral red-shift caused by the delayed thermal response.


Introduction and motivations
The design and development of subwavelength photonic devices with metallic components has become a subject of intense research over the last decade. This trend is justified by the need for compact high-performance optical devices and is mainly driven by the enormous technological improvement in nano-fabrication techniques. These state-of-the-art manufacturing tools for metallic nano-circuits have made it possible to design and engineer the effective optical properties of artificial materials, commonly known as metamaterials. In these synthetic materials, the propagation of light is strongly influenced by the geometric properties of the embedded metallic nano-circuits. In particular, the metallic nano-structures can be tailored in such a way that the effective refractive index becomes negative [1]. Negative index materials are potentially important in superlensing [2] and cloaking applications [3]. Novel physical mechanisms occur in anisotropic metamaterials, which under some circumstances can exhibit hyperbolic dispersion [4]. In the nonlinear regime, metamaterials have been studied for harmonic generation [5], soliton propagation [6] and optical modulation and switching [7]. The unusual properties of metamaterials arise from the fact that the metallic nano-circuits are much smaller than the wavelength of light, resulting in a space-averaged macroscopic dielectric response. Conversely, in the case where the optical wavelength is comparable with the dimensions of the metallic sub-structures, the light feels the geometric details and plasmon polariton modes are excited. Surface plasmon polaritons (SPPs) are electromagnetic waves propagating on metallic surfaces. They constitute the best candidates for manipulating light on the nanoscale and for the development of subwavelength all-optical devices [8]. In particular, plasmonic waveguides have important applications as optical interconnects in highly integrated optoelectronic devices [9]. Other relevant applications of SPPs are found in medicine [10], sensing [11] and nano-lasers [12]. The nonlinear properties of SPPs can be used for second harmonic generation [13], active control [14] and nanofocusing [15]. Nonlinear self-action can be exploited for manipulating transverse spatial diffraction by selffocusing and for the formation of plasmon solitons [16,17]. Subwavelength modulational instability, plasmon oscillon formation and nonlinear switching can be achieved in arrays of metallic nanoparticles [18,19]. Fundamental studies of metamaterials and SPPs are closely related. Indeed, relevant phenomena occurring in metamaterials are observed also in plasmonics, e.g. negative refraction [20], anomalous diffraction [21,22] and electromagnetic cloaking [23]. In both fields, an innovative step is the use of nanostructured metals for manipulating light. In most of the nonlinear studies reported above, the optical response of metals is assumed to be linear, while the nonlinearity originates from the dielectric medium; however, experimentalists know well that the Kerr nonlinearity of metals can be enormous. Experimental results indicate strong third-order nonlinear susceptibilities that vary by several orders of magnitude, with values of χ m 3 that vary between 10 −14 and 10 −18 m 2 V −2 [24,25] and that are much bigger than the third-order susceptibility of bulk silica (χ Si 3 ≈ 10 −22 m 2 V −2 ). Recently, the non-local ponderomotive nonlinearity for a plasma of free electrons was proposed as a possible model for the interpretation of experimental results [26]. The predicted value for the ponderomotive third-order susceptibility at optical frequencies (χ 3 ≈ 10 −20 m 2 V −2 ) is, however, insufficient to explain the experimental findings. In addition, the spectral dependence of the ponderomotive nonlinearity (χ 3 ∝ 1/ω 4 ) does not fit with the enormous spectral variation (by several orders of magnitude) observed in the measurements, suggesting that the basic nonlinear mechanism for metals is resonant. Theoretical and experimental confirmation of this hypothesis is to be found in the results of Rosei, Guerrisi et al on the thermo-modulational reflection spectra of thin films of noble metals [27,28]. In their work, the authors theoretically predict and experimentally observe a strong modulation in the reflection spectrum due to light-induced heating. They demonstrate that the temperature change smears out the energy distribution of the conduction electrons, affecting the resonant interband absorption and hence the dielectric susceptibility. This process is intrinsically nonlinear, since the temperature change modulating the dielectric response depends on the optical power. Subsequent pump-probe experiments in thin films [29] and nanoparticles [30] have confirmed the initial results of Rosei and Guerrisi. Theoretical and experimental investigations on the temporal dynamics of the system clearly indicate that the nonlinear response of metals is characterized by a delayed mechanism [31], as is typical of thermal nonlinearities. Very recently, a complete analysis of the nonlinear optical response of noble metals, leading to the first theoretical derivation of a consistent model for the third-order nonlinear susceptibility of gold, was reported [32]. Although experiments in thin films have been satisfactorily explained, a theoretical description of the thermo-modulational interband nonlinearity for ultrashort optical pulses propagating in plasmonic waveguides is still missing.
In this paper, we describe the effects of the thermo-modulational nonlinear susceptibility on SPPs propagating in a gold nanowire surrounded by silica glass. The paper is organized as follows. In section 2 we briefly describe the optical properties of gold, the interband transitions, their effect on the dielectric susceptibility and its temperature dependence. Then, in section 3, we model the temporal dynamics of the electrons through the two-temperature model (TTM), deriving the characteristic temporal response function. In section 4, collecting the results of the previous sections, we are able to derive an analytical expression for the thermo-modulation interband polarization. Finally, in section 5, we model the propagation of SPPs along a gold nanowire by introducing a novel nonlinear Schrödinger-like equation and predicting for the first time to our knowledge a strong red-shift caused by the thermo-modulational nonlinearity of gold.

Optical properties of gold
The non-resonant optical properties of gold can be described through the free-electron model, where electrons are considered to be free charges moving in response to an optical field Re[ E 0 e −iωt ] oscillating at angular frequency ω. In this model, the dielectric response of the plasma can be derived directly from the non-relativistic single-particle equation of motion: where ω p = ne 2 / 0 m e is the plasma frequency, 0 is the vacuum permittivity, n is the electron number density, e, m e are the electron charge and mass and γ is a characteristic frequency accounting for electron-electron collisions. This model is justified by the fact that for metals the Fermi energy lies within the conduction band and many accessible states exist for the electrons. From a quantum perspective, free-electron motion only accounts for intraband transitions. For wavelengths in the far infrared, the free-electron model provides very good quantitative agreement with experimental data [33]. For the free-electron calculations we use the parameters ω p = 1.1515 × 10 16 rad s −1 , γ = 8.9890 × 10 13 s −1 , obtained by fitting intra (ω) to the experimental data for long wavelengths in the far infrared. At optical frequencies and in the near infrared, the susceptibility of gold is more involved and the non-resonant model is not sufficient to explain the experimental measurements. Indeed, interband transitions between the d-band and the conduction band become important and cannot be neglected. The presence of interband transitions enriches the variety of physical processes and lies behind the strong temperature dependence of the dielectric susceptibility. At optical frequencies, the gold dielectric susceptibility measured in experiments [33] deviates significantly from the predictions of the free-electron model as a consequence of two intense interband absorption peaks at λ = 300 and 410 nm. Thus, the actual dielectric constant of gold can be expressed as the sum Gold is characterized by a face centered cubic lattice structure, where the Wigner-Seitz primitive cell is a rhombic dodecahedron. The reciprocal lattice is a body centered cubic where the Wigner-Seitz primitive cell is a truncated octahedron. At optical frequencies, the interband absorption is resonant around the points X and L in reciprocal space [28], which correspond, respectively, to the centers of the square and hexagonal facets of the truncated octahedron. Notably, such points are highly symmetric and around them the lattice vector k can be expressed as the sum k = k ⊥ + k , where k ⊥ lies on the square (X ) or hexagonal (L) facets and k is perpendicular to them. In addition, the Fermi surface has cylindrical symmetry around the X, L points and the valence and conduction bands can be approximated in (k ⊥ , k ) space by elliptic and hyperbolic paraboloids [28]: where v, c indicate the valence and conduction bands. The values of the constants E 0v , E 0c , m v⊥ , m c⊥ , m v , m c for the X and L transitions that we use in our calculations are known in the literature [28]. Note that, as a consequence of the cylindrical symmetry of the Fermi surface, the conduction and valence bands around the X, L points do not depend on the direction of k ⊥ , but only on its modulus k ⊥ . The Fermi level E F lies in the conduction band E c ( k) and for the sake of simplicity (and without loss of generality) we perform a constant shift of all energies and assume that E F = 0. Note that the paraboloid approximation made in (1) is accurate only if k ⊥ , k k , where k is the distance between the X, L points and the center of the Brillouin zone. The quantum states of electrons in the valence and conduction bands are Bloch wavefunctions ψ v,c = −1/2 u k,v,c ( r ) exp(i k · r ), where is the primitive cell volume, k is the Brillouin wavevector and u k,v,c ( r ) are the unit cell wavefunctions. By using the Fermi golden rule and following the procedure reported in [27,28,32], it is possible to calculate the imaginary part of the dielectric constant due to the X, L interband transitions where is the dipole matrix element and J c,v (ω, T e ) is the joint density of states (JDOS) In the equation above, T e is the electronic temperature and f (E, T e ) is the Fermi-Dirac occupation number. Note that the conservation of energy and the Pauli exclusion principle for every v → c direct transition are carefully considered in the expression for J c,v . Indeed, the factor f [E v ( k), T e ] accounts for the probability that the k-state in the valence band is occupied, while the factor 1 − f [E c ( k), T e ] accounts for the probability that the k-state in the conduction band is empty. The real part of the dielectric constant can then be obtained directly from (2) using the Kramers-Kronig relation where P represents the principal value of the integral. If the conduction electrons are taken out of equilibrium by light-induced heating, the interband absorption is affected by the so-called Fermi smearing effect [28]. Increasing temperature broadens the electron distribution around the Fermi energy, modifying the effective optical properties of gold. In order to understand the temperature dependence of inter , contained in the X, L JDOSs, one needs to compute numerically the thermo-derivatives [28,34]. The resulting interband dielectric thermo-derivative can be fitted to a series of five Lorentzian functions: where ω p is the plasma frequency of gold. The fit parameters are given in table 1. The calculation of the real part ∂ T e inter through the integration of the Kramers-Kronig relation given by (4) is straightforward: As light impinges on the gold surface, the electrons in the conduction band are heated and the dielectric constant is modified by the amount m (ω) = ∂ T e inter (ω) T e . Note that this thermo-modulational process is intrinsically nonlinear, since the increase of temperature depends on the absorbed optical power T e (P A ). The spectral dependence of the Table 1. Fit parameters for the thermo-derivative ∂ T e inter (ω), given by (5). ω p = 1.1515 × 10 16 rad s −1 is the plasma frequency of gold, which has been calculated by fitting the experimental data in the far infrared with the free-electron model. for several values of temperature variation T e = 20, 50, 100 K (blue, cyan and red curves, respectively). Note that the spectral dependence of both the real and imaginary parts m (ω), m (ω) is non-trivial and they can be either positive or negative; hence the optical absorption can increase or decrease, depending on the wavelength λ.

Electron temporal dynamics and the two-temperature model
An optical beam impinging on a metal surface modifies the effective interband susceptibility by heating the electrons in the conduction band. This light-induced electron heating can be described through the TTM, which takes into account the energy balance between the conduction electrons and the lattice. The electrons have a relatively small heat capacity and so thermalize through electron-electron collisions with a characteristic time of order τ th ≈ 300 fs. If one wishes to describe the temporal electron dynamics for ultrashort optical pulses (τ 0 ≈ 100 fs), it is necessary to include also the energy contribution of the non-thermalized electrons in the energy balance. This can be calculated directly from the Boltzmann equation in the relaxation time approximation [31]. A phenomenological description of the electron temporal dynamics can be obtained by separating the electron distribution of energy into thermalized and non-thermalized parts: where P A (t) is the mean absorbed power per unit volume, T e (t), T l (t) are the electronic and lattice temperatures, C e , C l are the electronic and lattice heat capacities per unit volume and N (t) is the energy density stored in the non-thermalized part of the electronic distribution. When an ultrashort optical pulse impinges on the metal, it is absorbed and transfers energy to the non-thermalized electrons. In turn, the non-thermalized electrons release energy density γ e N (t) to the thermalized electrons via electron-electron scattering and energy density γ l N (t) to the lattice via electron-phonon scattering. The non-thermalized electrons achieve thermal equilibrium with a characteristic time delay of τ th = (γ e + γ l ) −1 , where γ e , γ l are the electron and lattice thermalization rates. Once heated by an ultrashort optical pulse, the thermalized electrons gradually release energy to the lattice via electron-phonon scattering, which is accounted for by the coupling coefficient C. Ultimately, for long times, the electrons reach thermal equilibrium with the lattice. The parameters used in the numerical calculation are known in the literature [29]. Note that since the lattice heat capacity is much larger than the electron heat capacity, while the temporal variation of the electronic temperature T e (t) is significant, the lattice temperature T l (t) does not change significantly with time, i.e. T l (t) ≈ const. Note also that the electronic and lattice heat capacities C e , C l , in principle, depend on their respective temperatures T e , T l so that the extended TTM model is nonlinear. However, in the limit T e (t) T eq , it is possible to approximate C e , C l as independent of their respective temperatures [31]. Setting ∂ t N = 0 and removing the equation for N (t) is equivalent to neglecting the thermalization time τ th ≈ 300 fs over which non-thermalized electrons release energy to the thermalized ones. For pulses of duration τ 0 ≈ 100 fs such an approximation is not feasible and all the equations must be retained in order to describe correctly the delayed nonlinearity. The TTM model can be solved straightforwardly in the Fourier domain, leading to the solution where ω is the shift from the carrier angular frequency ω 0 of the ultrashort optical pulse and In the expressions above we have used the parameter τ r = C e C l /[C(C e + C l )] representing the relaxation time of the thermalized electrons with the lattice. In the temporal domain, the temperature variation T e (t) = T e (t) − T l (t) is given by where h is the temporal response function. Note that causality is imposed by means of the Heaviside step function θ(t). Note also that in the limit P A → 0 the conduction electrons and the lattice are in equilibrium and have the same temperature T e = T l = T eq . Since the heat capacity of the lattice C l is much greater than the electronic heat capacity C e , the temperature difference T e (t) = T e (t) − T l (t) can be approximated by T e (t) ≈ T e (t) − T eq . The response functions in the frequency h T ( ω) and temporal h T (t) domains are plotted as functions of ω and t in figures 2(a) and (b). In figure 2(a), the blue and red curves correspond to the real h T ( ω) and imaginary h T ( ω) parts. Note that, following the sign convention chosen for the exponential in the Fourier expansion (e −i ωt ), a positive value of the imaginary part h T ( ω) corresponds to loss, while a negative value corresponds to gain. The temporal thermal response h T (t), depicted in figure 2(b), is mainly characterized by a peak delayed in time by t ≈ 600 fs. As a consequence of the delay, blue-shifted frequency components are suppressed, while red-shifted components are amplified, analogously to what happens in solid-core optical fibers as a result of the Raman effect.

Thermo-modulational interband nonlinear susceptibility
In this paragraph we sum up the results obtained so far calculating the nonlinear susceptibility due to Fermi smearing of the conduction electrons. The instantaneous power per unit volume absorbed by a metal is W (t) = E(t) · ∂ t D(t), where D(t) is the displacement vector. In the continuous wave (CW) monochromatic case, the electric and displacement fields are where ω 0 is the angular frequency. The mean absorbed power per unit volume can be calculated by averaging the instantaneous power over the fast oscillations e −iω 0 t : Hence, the nonlinear polarization is given by P CW is the thermo-modulational interband nonlinear susceptibility: and γ T (ω) = τ r τ th (γ e /C e − γ l /C l ) ∂ T e inter (ω). The real and imaginary parts of χ (3) T are plotted as functions of optical wavelength λ in figure 3. Note that, as a consequence of the resonant interband transitions, the nonlinear susceptibility is strongly dispersive at optical frequencies and can be much greater (≈7 orders of magnitude) than the Kerr susceptibility of bulk silica (χ Si 3 ≈ 10 −22 m 2 V −2 ). The strong frequency dispersion of gold dramatically changes its optical properties, as well as the signs of Re χ (3) T , Im χ (3) T . Note that, for wavelengths λ 750 nm,  can be approximately considered to be a constant is also indicated.
the thermo-modulational nonlinear susceptibility can be approximated as χ (3) T ≈ const. For ultrashort optical pulses, the calculation of the nonlinear dielectric polarization is more involved.
In the slowly varying envelope approximation (SVEA) the electric field can be expressed as E(t) = ψ(t)e −iω 0 tn , where ω 0 is the carrier angular frequency,n is the polarization unit vector and ψ(t) is the envelope amplitude, which is slowly varying compared to the fast oscillations e −iω 0 t . The expression for the mean absorbed power per unit volume P A (t) in this non-monochromatic case includes also the contributions of first-order dispersion and is explicitly given by Note that, in the limiting case of infinite pulse duration (ψ(t) = ψ 0 ), the expression for P A converges to the CW result, which is the first term in the expression above (∝ m (ω 0 )|ψ| 2 ). The second term (∝ ∂ t |ψ| 2 ) does not contribute to the total density of energy absorbed U A = +∞ −∞ P A (t)dt, since |ψ| 2 → 0 for t → ±∞. The last term (∝ i(ψ * ∂ t ψ − ψ∂ t ψ * )) depends on the spectral momentum and basically represents the dispersive correction of the first term. In other words, it accounts for the modification of the absorbed power per unit volume P A within the frequency spectrum of the pulse. In conclusion, the nonlinear polarization created by an ultrafast optical pulse can be expressed in terms of a double convolution integral In this expression, h T (t) is the temporal response function of the TTM and γ T (t) is the interband response function (measured in the units of m 3 W −1 s −1 ), which is given by the inverse Fourier transform of γ T (ω) (measured in the units of m 3 W −1 ).

Thermo-modulational nonlinear dynamics in plasmonic devices
In the previous sections we have calculated the thermo-modulational interband nonlinear susceptibility of gold starting from the basic properties of its band structure. In this section we study the optical propagation of SPPs guided along gold nanowires surrounded by silica glass, including the novel nonlinear effects originating from the heating of gold. A common theoretical approach to modeling optical propagation in optical fibers and plasmonic waveguides uses the nonlinear Schrödinger equation for the slowly varying amplitude of a guided pulse, perturbatively derived on the assumptions of low loss and nonlinearity [35,36]. We start the analysis from the time-dependent Maxwell equations for the optical electric ( E) and magnetic ( H ) fields, coupled with the nonlinear polarization P NL ( r , t). As a consequence of the cylindrical symmetry of the gold nanowire around the z-axis, the dielectric profile depends solely on the modulus of the position vector ρ = | r |: L (ρ, ω) = m (ω)θ(r − ρ) + d (ω)θ(ρ − r ), where θ (x) is the Heaviside step function, ω is the angular frequency, r is the radius of the nanowire and m (ω), d (ω) are the linear dielectric constants of gold and silica.
In the calculations below we use the Sellmeier expansion for the dielectric constant of silica d (ω) and a Lorentzian fit to the experimental data [33] for the linear dielectric constant of gold m (ω). A sketch of the plasmonic structure discussed in this section is depicted in figure 4.

Linear modes
If one neglects the nonlinear polarization P NL , the time-dependent Maxwell equations can be directly solved using the ansatz: where c is the speed of light in vacuum, ψ is the mode amplitude (measured in W 1/2 ), e, h are the linear guided mode profiles (dimensionless), β is the mode propagation constant, φ is the angle between the vectors r ,x and m is the azimuthal mode order. The factor I 1/2 is a constant, chosen in such a way that |ψ| 2 represents the total optical power carried by a linear mode with angular frequency ω and propagation constant β. For both the m = 0, 1 modes, if λ 500 nm, the real (Re β) and imaginary (Im β) parts of the propagation constant increase as the optical wavelength λ decreases, reaching maxima at λ = λ sp 500 nm, the surface plasmon resonance. For λ 500 nm, the behavior of the dispersion relation is more involved. In the ideal case where the metal is lossless, the imaginary part of the propagation constant vanishes, Im β = 0, whereas the real part Re β diverges as λ → λ sp . Optical confinement depends mainly on the real part of the propagation constant Re β: high values of Re β correspond to tightly confined modes. On the other hand, the attenuation coefficient of the linear modes is directly related to the imaginary part of the propagation constant: α = 2 Im β. Hence, if one wants to achieve tight SPP confinement it is impossible to avoid high losses so that the use of materials with large gain is required in practical applications [12]. Note that the fundamental mode m = 0 is TM polarized and that both the real and imaginary parts of the propagation constant Re β, Im β increase as the radius of the gold nanowire decreases (see figures 5(a) and (b)). A contour plot of the time-averaged Poynting vector S z = (1/2)ẑ · Re( E × H * ) of the TM plasmonic mode   (m = 0) of a gold nanowire with radius r = 50 nm at optical wavelength λ = 800 nm is depicted in figure 6(a). Note that the electromagnetic field is tightly bound to the metal surface and that the power distribution is cylindrically symmetric.
In contrast to the TM fundamental mode, the m = 1 mode is hybrid polarized and less well confined. The dispersion relation does not depend on the sign of the azimuthal mode order m so that the m = ±1 modes (characterized by opposite chirality) are degenerate. A contour plot of the time-averaged Poynting vector S z of the superposition of m = ±1 guided SPP modes is shown in figure 6(b) (for the same parameters of figure 6(a)). Note that the power distribution of such a mode is not azimuthally symmetric and depends on the angle φ. If one calculates the time-dependent Poynting vector of the m = ±1 modes, one finds that these SPPs spiral around the surface of the gold nanowire with opposite chirality [37]. However, the azimuthal spiralling is averaged out in time and as a result there is no net angular flow of time-averaged power for the single m = ±1 modes and for their superposition. In figure 7, Im β for the m = 1 mode is plotted as a function of r . Red, green and blue curves correspond to λ = 700, 800 and 900 nm, respectively. Both Re β and Im β depend on r in a manner more complicated than for the fundamental TM mode (see figures 5 and 7). In particular, for fixed optical wavelength λ, Im β is maximum at a characteristic wire radius r 0 , decreasing significantly when r < r 0 . Hence, if the wire radius is much smaller than the optical wavelength r λ, long-range SPPs can be excited [37]. The field penetration within the gold nanowire is limited for these modes, and hence the attenuation is significantly reduced. In the following nonlinear analysis, we focus on the m = 0, 1 modes. Higher order modes are less well confined and cut off for wavelengths greater than a particular value λ co .

The generalized nonlinear Schrödinger equation
If the effect of loss and nonlinearity on the fast linear oscillations is weak, the nonlinear propagation of an optical pulse in a plasmonic waveguide can be described by the generalized nonlinear Schrödinger equation (GNLSE) for the field amplitude in the SVEA [35,36]. In this approach the ansatz for the electromagnetic field is where ψ(z, t) is the slowly varying envelope amplitude, ω 0 is the carrier angular frequency and δ E, δ H are the residual field corrections due to loss and nonlinearity. We denote β 0 by the linear propagation constant at the carrier frequency calculated by neglecting the nonlinear polarization and the metal loss; e(ρ), h(ρ) are the corresponding unperturbed linear mode profiles (dimensionless) and I is a constant chosen so that |ψ| 2 represents the optical power, as in the previous section. For a gold nanowire surrounded by silica glass, the nonlinear polarization is is the Kerr coefficient of silica glass at λ 0 = 800 nm and P Au NL ( E) is given by (12). Note that in the expression for the nonlinear polarization of silica glass, we have neglected the Raman effect [38], which should, in principle, be retained. However, as we will show, the effective nonlinear coefficient of silica is much smaller than the effective nonlinear coefficient of gold so that neither Kerr nor Raman effects play any significant role. Hence, for the sake of simplicity, we do not consider the Raman term, comparing our results only with the Kerr term. Since we are mainly interested in the thermo-modulational interband nonlinearity, we neglect the dispersive terms in the absorbed power per unit volume P A (t), given by (11). These terms are expected to play only a minor role, since they are small corrections to the carrier term. In what follows, we will focus on the spectral region (λ ≈ 800 nm) where the interband response function is approximately constant (γ T (ω) ≈ γ T (ω 0 ), see figure 3) so that the nonlinear polarization of gold can be approximated by P Au By inserting (15) and (16) into the Maxwell equations, developing a first-order perturbative theory and imposing the Fredholm alternative condition on the residual field corrections δ E, δ H [36], one obtains the GNLSE for the slowly varying amplitude ψ(z, t): where The linear dispersion operatorD(i∂ t ) is complex, accounting as it does for the linear losses of gold. Its action on the envelope amplitude can be calculated in the Fourier domain: where ω = ω − ω 0 and D ω 0 ( ω) = β(ω 0 + ω) − β 0 − dβ dω ω 0 ω. Note that β 0 is the realvalued carrier propagation constant of the linear unperturbed mode, β(ω) is the complex modal wavevector calculated in the previous section and the prime in the equation above indicates the real part (β = Re β). The nonlinear parameters ϒ Si , ϒ Au are measured in m −1 W −1 and account also for the surface nonlinearity [36]. Note that while ϒ Si is a real quantity, ϒ Au is complex and accounts for the nonlinear loss of gold. The nonlinear parameters of silica (ϒ Si ) and gold (ϒ Au ) are plotted as functions of the carrier wavelength λ 0 in figures 8, 9(a) and (b). In both figures, the full and dashed curves represent the m = 0 and 1 modes, while blue and green colors correspond to wire radii r = 50 and 250 nm. Note that the real part of the gold nonlinear parameter is much greater than the Kerr nonlinear parameter of silica in the spectral region considered. Also, if r λ, the nonlinear parameters of the m = 1 mode are much smaller than those for the m = 0 mode since they are much less confined. In this limit, as discussed in the previous section, while the fundamental m = 0 mode is tightly confined to the metal surface and propagates only for a few wavelengths, the m = 1 mode is much less localized and can  propagate for longer distances (long-range guided SPP mode). This reduction in loss is at the cost of a weaker effective nonlinearity.
The formulation of the propagation equation (17) constitutes the main result of this paper. We have numerically solved (17) using the fast Fourier split-step algorithm. As we have already shown, when λ 0 ≈ 800 nm the Kerr nonlinearity of silica does not play a significant role and can be neglected. For the numerical simulations we have considered a hyperbolic secant input pulse: ψ(0, t) = √ P in sech(t/t 0 ), with t 0 = 106 fs. P in is the instantaneous pulse power, which can be directly calculated from the average power of the laser source: P in = C eff P av /(2ν rep t 0 ), where C eff is the launch efficiency of the laser beam into the gold nanowire, ν rep is the repetition rate and 2t 0 is the pulse duration. The instantaneous power is kept well below the damage threshold power of gold, above which melting, ablation and vaporization occur (P dam ≈ 10 6 W) [39]. In figures 10(a) and (b), the numerical propagation along a gold nanowire with radius r = 50 nm surrounded by silica glass is depicted for (a) m = 0, P in = 1 × 10 4 W and (b) m = 1, P in = 5.3 × 10 5 W. In this contour plot the modulus of the Fourier transform of the optical amplitude (|ψ(z, ω)|) is shown. The m = 0 TM mode is highly nonlinear and significant nonlinear dynamics can be observed even for relatively small optical power. However, such a high nonlinearity is paid for by high loss, limiting the effective propagation length to L ≈ 2 µm (see figure 10(a)). The hybrid polarized m = 1 mode supports long-range guided SPP modes so that both the nonlinear and loss coefficients are significantly smaller. In this case, in order to observe a strong nonlinear dynamics, it is necessary to use a considerably higher optical power. For the m = 0, 1 modes a signature red-shift indicates the presence of a thermo-modulational interband nonlinearity, analogously to the Raman effect [38].
This red-shift is the natural consequence of the intrinsic delayed mechanism governing the thermo-modulational interband nonlinear susceptibility of gold. In the time domain the frequency red-shift is accompanied by a small pulse delay of order ≈1 fs. We emphasize that neither the Kerr nor the Raman nonlinearities of silica are large enough to produce the reported red-shift for the propagation lengths considered. The strong red-shift is accompanied by a large time-delayed nonlinear loss, as can be understood from figures 11(a) and (b). In figure 11(a), the transmission spectrum (T = ln|ψ(L , ω)/ψ 0M |, where ψ 0M = max[ψ(0, ω)] and L = 100 µm) of the m = 1 long-range mode of a gold nanowire with radius r = 50 nm is depicted for several input powers: P in = 5.3 × 10 4 W (blue curve), P in = 2.7 × 10 5 W (green curve) and P in = 5.3 × 10 5 W (red curve). The black dashed curve corresponds to the normalized input spectrum on a logarithmic scale (ln|ψ(0, ω)/ψ 0M |). Note that as the input power increases the red-shift increases and the transmission peak decreases accordingly as a consequence of nonlinear loss. Also, as the power increases, the transmission spectrum displays some weak oscillations, which resemble Kerr-related self-phase modulation. Indeed, the dispersion length is much longer than the nonlinear length so that only the linear loss affects the nonlinear dynamics. In figure 11(b), the transmission spectrum (T = ln|ψ(L , ω)/ψ 0M |, where L = 20 µm) of the m = 1 long-range mode is shown for several nanowire radii: r = 50 nm (blue curve), r = 100 nm (green curve) and r = 400 nm (red curve). The black dashed curve represents the input spectrum (ln|ψ(0, ω)/ψ 0M |) and the input power is fixed at P in = 5.3 × 10 4 W. Note that the linear properties of the m = 1 mode are non-trivial (see figure 7) and, as a consequence, the power dependence of the red-shift and the transmission peak is non-monotonic. This means that an optimal radius r = r o exists, where the achievable red-shift is maximum. In figure 12, the redshift of the transmission peak of the m = 1 long-range mode is plotted as a function of wire radius for a fixed input power P in = 5.3 × 10 4 W and for a propagation distance L = 20 µm. The open blue circles represent the results of numerical simulations, whereas the black curve corresponds to an interpolation of the numerical results. The maximum red-shift attainable at this input power is λ ≈ 7 nm for a radius of r o ≈ 110 nm. The dependence of the thermomodulational interband nonlinear coefficient (ϒ Au ) and the absorption coefficient (α = 2β 0 ) on r strongly affects the nonlinear dynamics and as a consequence the red-shift attainable, which reaches a maximum at r = r o . When the input power increases, the maximum red-shift increases accordingly.

Summary
In this paper we have described the thermo-modulational interband nonlinearity of gold starting from its band structure. Electrons in the conduction band are heated by an ultrashort optical pulse and the interband dielectric properties are modulated accordingly. Using a semiclassical approach, we have calculated the imaginary part of the dielectric constant of gold, accounting for interband absorption, which basically depends on the JDOS. In turn, we have been able to describe the temperature dependence of the dielectric susceptibility of gold and have modeled the heating and the relaxation of the conduction electrons using a TTM. We discovered that the metal nonlinearity is basically characterized by a delayed mechanism, similar to the Raman effect in some senses, but with a much longer response time (≈300 fs). Also, in contrast to the Raman effect, the thermo-modulational interband susceptibility is complex valued, providing a delayed nonlinear loss or gain. The optical propagation of SPPs is strongly affected by the metal nonlinearity, which we have found to be several orders of magnitude larger than the Kerr nonlinearity of fused silica. We have derived, for the first time to our knowledge, a GNLSE suitable for modeling the optical propagation of SPPs along a gold nanowire surrounded by silica glass. Solving this equation using a fast Fourier split step algorithm, we found that the signature of the thermo-modulational interband nonlinearity is a red-shift of the optical pulse. This red-shift results from the intrinsic time-delayed nature of the thermo-modulational interband nonlinearity of gold. We believe that this novel nonlinear effect may be important for frequency conversion in plasmonic devices. We have also provided some details of the expected red-shift, its dependence on the wire radius and the optical power necessary to observe it experimentally.