Nonlinear thermoplasmonics in graphene nanostructures

Nonlinear optical phenomena facilitate the temporal and spectral control of light by light, a crucial resource for emerging integrated photonic and quantum optical technologies [1, 2]. Despite the inherently weak nonlinear response of conventional materials [3], demand for efficient miniaturized photonic devices motivates intensive research on material platforms and schemes to enable nonlinear light-matter interactions on the nanoscale [4, 5, 6]. Graphene is actively explored in this context due to its conical electronic dispersion relation, which endows the atomically thin carbon layer with a proportionately large nonlinear optical response [7]. Measurements of third-harmonic generation (THG), four-wave mixing, and the optical Kerr effect in graphene reveal large third-order susceptibilities, on the order of $10^{-16}-10^{-15}$ m²/V², several orders of magnitude larger than in other highly nonlinear optical materials [8, 9, 10, 11, 12].

In electrically doped graphene, the hybridization of light with plasmons—collective oscillations of conduction electrons—can strongly concentrate electromagnetic near fields to boost nonlinear optical effects. In particular, the combined linear electronic bands and atomic thickness of graphene results in plasmon resonances with high sensitivity to charge carrier doping levels [13], enabling electrical tunability of plasmon resonances in the terahertz (THz) or infrared (IR) spectral bands [14]. Patterning of graphene into ribbons or islands offers the means to couple light with localized plasmon resonances that can be tuned passively through spatial confinement while maintaining active tunability via electrostatic gating [15]. While enhancements in the nonlinear response associated with THz plasmons in graphene nanoribbons (GNRs) have been reported in wave-mixing experiments [16], large optical nonlinearities driven by plasmons at near-IR and visible frequencies in smaller graphene nanostructures, predicted to exhibit an involved dependence on quantum finite-size and edge termination effects [17], have yet to be observed. In this context, narrow GNRs of $\lesssim 10$ nm width constitute an appealing platform for explorations in nonlinear graphene plasmonics, particularly with recent advances in nanofabrication techniques enabling atomically-precise control of GNR width and edge configuration [18, 19, 20, 21].

The linear electronic bands in graphene are not only associated with a large parametric nonlinear optical response, but also a low electronic heat capacity [22, 23], thus presenting opportunities to explore nonlinear dynamics associated with thermo-optical effects [13, 24]. Energy absorbed by pumping graphene with intense ultrashort optical pulses is rapidly converted to electronic heat within $\sim 10$s of femtoseconds [25, 26, 27], such that the two-dimensional electron gas can reach temperatures of several $\sim 1000$ K using pulses of moderate $\lesssim 1$ J/m² optical fluence [28]. Time- and angle-resolved photoemission spectroscopy measurements reveal that the thermalized carrier distribution relaxes on $\sim$picosecond timescales by exciting phonons in the carbon lattice [29, 30]. The collective thermal response of graphene under illumination by THz pulses—for which electronic heating is effectively instantaneous relative to oscillations of the electromagnetic field—results in efficient harmonic generation, with third-harmonic susceptibilities reaching $\sim 10^{-9}$ m²/V² [31]. At IR frequencies, nonequilibrium intra-pulse dynamics in graphene can manifest as spectral blueshifts in third- and fifth-harmonic signals [32], while carrier heating effects have been found to significantly impact wave mixing and harmonic generation [33, 12].

Electron heating in photoexcited graphene facilitates optical modulation of graphene plasmon resonances on ultrafast timescales [34, 35]. For example, experiments on graphene encapsulated in hBN show that optically-induced carrier heating can sustain mid-IR plasmons that are otherwise absent in the equilibrium state [36], while pump-probe measurements of THz plasmons in doped GNRs reveal a transient redshift in localized plasmon resonances [37]. Qualitatively similar spectral shifts and saturation of mid-IR plasmons under resonant excitation are expected for narrower GNRs of a few $10$s of nm width [38], while the distribution of absorbed optical energy among fewer free charge carriers in even smaller graphene nanostructures is predicted to enable more dramatic temperature-induced changes in the optical response [39].

In this letter, we explore possibilities to thermally induce mid- and near-IR plasmons in GNRs with $\lesssim 10$ nm width that can trigger a large nonlinear response. Combined ab-initio modelling and self-consistent perturbation theory reveals that quantum finite-size effects in GNRs significantly impact their supported thermoplasmons, particularly in the nonlinear response associated with third-harmonic generation. Prominent thermoplasmonic features are predicted to emerge in the linear and nonlinear spectral response after doping GNRs with moderate heat energy densities on the order of $\sim 1$ J/m², below the fluences of optical pulses commonly used in ultrafast experiments. Our results indicate that optical pumping constitutes a viable approach to realizing a sizeable transient nonlinear plasmonic response at near-IR frequencies in mesoscopic graphene nanostructures, circumventing the need for cumbersome charge carrier doping schemes.

Electrons in a narrow GNR that are photoexcited by an ultrafast pump pulse with moderate fluence will rapidly thermalize to temperatures on the order of $\sim 1000$s K, effectively modifying the free charge carrier density within $\sim 1$ ps as the electrons transfer their energy to the carbon lattice. As illustrated schematically in Fig. 1a, transient plasmon resonances sustained by the thermalized electron distribution can be optically probed in both linear and nonlinear optical response regimes. Electronic states in armchair (AC) and zigzag (ZZ) edge-terminated GNRs of $\sim 5$ nm width obtained from first-principles density functional theory (DFT) calculations are plotted in Figs. 1b and 1c, respectively, as functions of the electronic wave vector in the direction of translational invariance (see Methods). The population of a state with energy $\hbar\varepsilon_{j,k}$ indexed by band $j$ and wave vector $k$ is determined from Fermi-Dirac statistics as

$$f_{j,k}(\mu,T)=\left[{\rm e}^{(\hbar\varepsilon_{j,k}-\mu)/{k_{\rm B}}T}+1\right]^{-1}$$

(1)

for a chemical potential $\mu$ and electron temperature $T$ satisfying the particle number conservation condition

$$\frac{b}{2\pi}\int_{-\pi/b}^{\pi/b}{\rm d}k\sum_{j}\mathopen{}\left[f_{j,k}(\mu,T)-f_{j,k}({E_{\rm F}},0)\right]\mathclose{}=0,$$

(2)

where the integration over $k$ is defined by the GNR unit cell period $b$ and ${E_{\rm F}}$ is the Fermi energy. The panels plotted alongside the electronic bands in Fig. 1b and c show the corresponding Fermi-Dirac distribution obtained for pristine GNRs, i.e., when ${E_{\rm F}}=0$ (solid curves), as well as for hole-(p-)doped GNRs at doping density $p=10^{14}$ cm^-2 (dashed curves), at the specified electron temperatures and corresponding chemical potentials that satisfy Eq. (2).

To study the optical properties of thermalized GNRs, we apply the ab-initio framework reported in Ref. [40] to construct tight-binding Hamiltonians for carbon p_z orbitals using maximally-localized Wannier functions (MLWFs). As shown in Figs. 1b and 1c, the Wannier tight-binding (WTB) models accurately reproduce the ab-initio bands corresponding to p_z orbitals at low energies around the Fermi level that are most relevant for electron temperatures $T\lesssim 10000$ K, as indicated by the distributions plotted alongside the electronic bands. The obtained WTB Hamiltonians are used to compute the linear and nonlinear optical response of GNRs according to a perturbative solution of the equation of motion governing the single-particle density matrix [41]. Further details on the ab-initio electronic structure calculations and optical response simulations are provided in Methods, while additional calculations for narrower and wider ribbons are provided in Supplementary Material (SM).

In Figs. 1d and 1e, we present the linear absorption cross-section and the nonlinear susceptibility associated with THG computed for the AC and ZZ GNRs considered in Figs. 1b and 1c, respectively, under the same electronic temperature and charge carrier doping conditions. Temperature-induced plasmon resonances appear in the linear response already at $T=2500$ K, which are slightly more pronounced in ZZ GNRs, while the thermoplasmons blueshift and become more prominent as the electron temperature is increased to $10000$ K. Comparison of the pristine and p-doped ribbons reveals relatively minor changes in the linear response associated with thermoplasmons due to charge carrier injection, indicating that the temperature smearing of the electron distribution plays a much more important role than the doping density. Interestingly, an analysis of smaller GNRs presented in the SM reveals that the suggested crossover from charge doping-induced plasmons to thermoplasmons consistently occurs at relatively low temperature. The THG response shown in the lower panels of Figs. 1d and 1e shows strikingly different behavior, e.g., a more noticeable dependence on the charge carrier doping, while the thermoplasmonic THG response of the ZZ GNR tends to decrease with increasing electronic temperature.

The large electronic temperatures required by thermoplasmons in GNRs are associated with thermal energy $Q$ that can be delivered by ultrashort optical pulses. From the electronic band structures computed for the one-dimensional GNRs considered here, we determine the chemical potential $\mu$ and electronic temperature $T$ by simultaneously imposing the particle number conservation condition in Eq. (2) and energy conservation

$$\frac{b}{2\pi}\int_{-\pi/b}^{\pi/b}{\rm d}k\sum_{j}\hbar\varepsilon_{j}\left[f_{j,k}(\mu,T)-f_{j,k}({E_{\rm F}},0)\right]=Q.$$

(3)

The evolution of $\mu$ and $T$ with the heat energy per unit area $Q/{\rm Area}$ is presented in Figs. 1f and 1g for AC and ZZ ribbons, respectively, in pristine (solid curves), electron-doped (density $n=10^{14}$ cm^-2, dotted curves), and hole-doped (density $p=10^{14}$ cm^-2, dashed curves) specimens. Independently of the considered charge doping conditions, electronic temperatures of several $1000$ K are attained for heat energy levels commensurate with pulse fluences used in ultrafast optical experiments [28].

In Fig. 2 we study the linear and nonlinear optical response of the $W\approx 5$ nm AC and ZZ GNRs considered in Fig. 1 when either thermal energy (Fig. 1a) or additional charge carriers (Fig. 1b) are added. Thermoplasmon resonances, appearing as prominent features in the linear absorption cross sections of Fig. 2a that blueshift with increasing heat energy, are relatively well-reproduced by the dashed curve, obtained from the combined intraband conductivity of extended graphene and classical electrodynamic model of plasmon resonances in GNRs as

$$\omega_{\rm p}=\frac{e/\hbar}{\sqrt{-\pi\eta_{1}\epsilon}}\sqrt{\frac{\mu_{\rm D}}{W}},$$

(4)

where $\mu_{\rm D}=\mu+2{k_{\rm B}}T\log\mathopen{}\left(1+{\rm e}^{-\mu/{k_{\rm B}}T}\right)\mathclose{}$ is the temperature-dependent Drude weight, $\epsilon=1$ is the relative permittivity of the environment, and $\eta_{1}=-0.06873$ is the eigenvalue associated with the leading dipolar plasmon resonance supported by the ribbon geometry [42]. Specifically, the dashed curve in Fig. 2a is given from Eq. (4) by setting $\mu={E_{\rm F}}=0$ and using the electron temperature obtained from the heat capacity of extended graphene, $Q/{\rm Area}=\beta(k_{\rm B}T)^{3}/(\hbar v_{\rm F})^{2}$, where $\beta\approx 2.2958$ in the ${E_{\rm F}}\ll{k_{\rm B}}T$ limit considered here [43]. The systematic overestimation of the thermoplasmon resonance energy in the classical model is attributed here to its omission of Landau damping by interband transitions, which, in the case of narrow GNRs, can become more prominent at specific frequencies due to the large energy gaps between electronic bands associated with quantum confinement. In particular, we find that thermoplasmons can strongly hybridize with single-particle transitions at fixed frequencies to create Fano-like spectral features in the absorption cross section. Incidentally, for the same added heat energy, the thermoplasmon energy in ZZ GNRs is reduced compared to their AC counterparts due to the presence of edge states.

Quantum finite-size effects are found to manifest more clearly in the thermoplasmon-driven THG spectra shown in Fig. 2a, which are not only sensitive to the superlinear dependence on plasmon field enhancement at the fundamental frequency that scales as $\sim\sqrt{\mu_{\rm D}}$, but also to single-particle transitions at higher harmonics. The intersection of these effects results in regions of pronounced enhancement and quenching in the THG spectra that are independent of the added heat energy. These findings suggest that, at particular frequencies, photoexcited transient thermoplasmons in narrow GNRs should produce a relatively large THG response over a longer duration as heat energy is dissipated to phonons in the carbon lattice. Qualitatively different behavior is observed in the optical Kerr effect, quantified by the intensity-dependent nonlinear refractive index $n_{2}$ and the two-photon absorption coefficient $\beta$, respectively, which evolve monotonically with the thermoplasmon feature in the linear absorption spectra. Notably, the heat-independent hybridization features seen in the linear response are more clearly resolved in the Kerr nonlinearity, for which the thermoplasmon-enhanced near fields at the driving and generated frequencies coincide. The behavior described above is corroborated by additional simulations of thermoplasmons in narrower and wider GNRs presented in the SM: in narrower GNRs, the thermoplasmonic response exhibits larger spectral “jumps” with varying heat energy, owing to the greater energetic spacing between electronic bands associated with quantum confinement effects; conversely, the reduced energy spacing of electronic bands in wider GNRs leads to a smoother spectral evolution of thermoplasmons with heat energy.

Compared to the thermoplasmons activated with $\lesssim 1$ J/m² heat energy, the plasmonic resonances supported by electrically-doped GNRs at room temperature appear as sharper spectral features in Fig. 2b, however requiring doping densities $p\gtrsim 10^{14}$ cm^-2 that may be difficult to achieve through electrostatic gating. In the linear absorption spectra, plasmon energies are similarly well described by the classical model of Eq. (4) taking $\mu=\hbar v_{\rm F}\sqrt{\pi p}$, where $v_{\rm F}\approx c/300$ is the Fermi velocity in graphene, albeit overestimating the plasmon resonance frequencies in the ZZ GNR owing to the presence of edge states. While the large absorption cross sections at plasmon resonances correspond to higher in-plane electric field intensities that lead to a proportionately larger THG response, the effective Kerr nonlinearity driven by doping charge carrier plasmons is actually reduced compared to that associated with thermoplasmons. We attribute this phenomenon to the interplay between enhanced Ohmic losses at the plasmon resonance and the sensitivity of the intrinsic Kerr nonlinearity to the electronic distribution in the GNRs, both of which are accounted for in the definitions of $n_{2}$ and $\beta$ (see Methods). Interestingly, the thermoplasmon-driven two-photon absorption coefficient is mainly positive, while the doping charge plasmon produces both positive and negative values of $\beta$, corresponding to TPA and saturable absorption, respectively.

We have thus far focused on the thermo-optical response of individual GNRs, while realistic pump-probe experiments should involve ribbon arrays to accumulate detectable far field signals from transient thermoplasmons [13, 37]. The fabrication of narrow GNR arrays, with nanometer-scale control over ribbon width and spacing [44], presents opportunities for engineering tunable plasmonic devices that exploit plasmon hybridization effects [45, 46, 47]. In Fig. 3, we study the linear and nonlinear optical response of $T=5000$ K thermoplasmons in undoped isolated AC and ZZ GNRs of selected widths $W$ and their arrays at a fixed 1 nm edge separation. The enhanced filling factors in arrays of wider ribbons lead to larger plasmon redshifts and absorption cross sections, while also boosting the associated THG and Kerr nonlinearity.

In conclusion, we have shown that narrow GNRs of $\lesssim 10$ nm width can support localized thermoplasmons at near-IR frequencies that enhance parametric nonlinear light-matter interactions. Thermal doping of narrow GNRs can be further exploited to study quantum confinement effects, both in the linear and nonlinear optical response regimes, with third-harmonic generation exhibiting a particularly large sensitivity to electronic structure. Our results indicate that thermoplasmons can be activated with moderate added heat energy, presenting opportunities to optically induce transient near-IR plasmons in small graphene nanostructures, without resorting to complicated electrical doping setups.

The electronic groundstates of one-dimensional AC and ZZ GNRs are computed using density functional theory, specifically using the Perdew-Berke-Ernzhofer functional as implemented in the Quantum Espresso code [48, 49, 50]. Wavefunctions are expanded in plane waves using normconserving pseudo potentials from the PseudoDojo database [51], with a cut-off energy of 86 Ry. We use 12 and 24 k-point grids in the direction of translational invariance for AC and ZZ ribbons, respectively, and a Gamma point sampling in the two directions across ribbons. Each ribbon is edge-passivated by hydrogen atoms, while the interatomic forces and cell stress are minimized. The isolated ribbons are separated by 3 nm in both lateral directions to avoid interactions between neighbouring images. The wannierization of the carbon p_z orbitals and creation of the Wannier tight-binding (WTB) Hamiltonians $\mathcal{H}^{\rm TB}$ are performed using the Wannier90 code [52], from which the energy $\hbar\varepsilon_{j,k}$ of band $j$ at Bloch wave vector $k$ is determined by diagonalizing

$$\sum_{l^{\prime}}\mathcal{H}_{ll^{\prime},k}^{\rm TB}a_{jl^{\prime},k}=\hbar\varepsilon_{j,k}a_{jl,k},$$

(5)

where $a_{jl,k}$ are expansion coefficients that weigh the contribution of orbital $l$ to the ribbon eigenstates.

Following the approach of Refs. [41, 40], the WTB eigenstates are used in the framework of the random phase approximation (RPA) to construct the non-interacting susceptibility

	$$\displaystyle\chi_{s\omega,ll^{\prime}}^{0}=\frac{e^{2}b}{\pi\hbar}\int_{-\pi/b}^{\pi/b}{\rm d}k\sum_{jj^{\prime}}\mathopen{}\left(f_{j^{\prime},k}-f_{j,k}\right)\mathclose{}$$		(6)
	$$\displaystyle\times\frac{a_{jl,k}a_{j^{\prime}l,k}^{*}a_{jl^{\prime},k}^{*}a_{j^{\prime}l^{\prime},k}}{\omega+{\rm i}/2\tau-\mathopen{}\left(\varepsilon_{j,k}-\varepsilon_{j^{\prime},k}\right)\mathclose{}}$$

that determines the induced charge density component

$$\rho_{l}^{ns}=\sum_{l^{\prime}}\chi_{s\omega,ll^{\prime}}^{0}\phi_{l^{\prime}}^{ns}+\rho_{{\rm NL},l}^{ns}$$

(7)

quantifying the response at order $n>0$ oscillating with harmonic $s$ ($|s|<n$) of the external field ${\bf E}^{\rm ext}{\rm e}^{-{\rm i}\omega t}+{\rm c.c.}$ in the perturbation series $\rho_{l}(t)=\sum_{n,s}\rho_{l}^{ns}{\rm e}^{-{\rm i}s\omega t}$, where

$$\phi_{l}^{ns}=\phi^{\rm ext}_{s\omega,l}\delta_{n,1}+\sum_{l^{\prime}}v_{ll^{\prime}}\rho_{l^{\prime}}^{ns}$$

(8)

is the self-consistent potential, accounting for the external potential $\phi^{\rm ext}_{s\omega,l}$ at linear order ($n=1$) and the Coulomb interaction $v_{ll^{\prime}}$ between sites $l$ and $l^{\prime}$, while $\rho_{{\rm NL},l}^{ns}$ is a nonlinear source term (nonvanishing for $n>1$) constructed from lower perturbation orders. In Eq. (6), the dependence on chemical potential and electron temperature enters the occupation factors $f_{j,k}$ according to Eq. (1), while inelastic scattering processes are described by a phenomenological damping parameter that we set to a conservative value of $\hbar\tau^{-1}=50$ meV throughout this work. Combining the above expressions, the density can be expressed as

$$\rho^{ns}=\left(1-\chi_{s\omega}^{0}v\right)^{-1}\cdot\mathopen{}\left(\chi_{s\omega}^{0}\cdot\phi^{\rm ext}_{s\omega}\delta_{n,1}+\rho_{\rm NL}^{ns}\right)\mathclose{}$$

(9)

in matrix notation. For normally-impinging fields polarized along the ribbon confinement direction $\hat{\bf x}$, i.e., $\phi^{\rm ext}_{s\omega,l}=-x_{l}E^{\rm ext}$, we compute the polarizability per unit length $\alpha_{s\omega}^{(n)}/L=\sum_{l}x_{l}\rho_{l}^{ns}/(E^{\rm ext})^{n}b$, where $b$ is the unit cell length along the ribbon. In practice, we use the parametrized Coulomb interaction $v_{ll^{\prime}}$ of carbon $p_{z}$ orbitals described in Ref. [53], which is modified in the case of ribbon arrays to also include the interactions of site $l$ with the corresponding sites $l^{\prime}$ in all other ribbons [45].

Here we define the linear absorption cross section from the polarizability as $\sigma^{\rm abs}=(4\pi\omega/c){\rm Im}\{\alpha_{\omega}^{(1)}\}$, while the nonlinear suceptibilities are calculated from $\chi^{(n)}_{s\omega}=\alpha_{s\omega}^{(n)}/LWd_{\rm gr}$ by estimating the graphene thickness $d_{\rm gr}=0.33$ nm as the interlayer spacing in graphite. In absorbing media, the nonlinear refractive index $n_{2}$ and two-photon absorption coefficient $\beta$ are given as [54]

	$\displaystyle n_{2}$	$\displaystyle=\frac{3}{4\epsilon_{0}c(n_{0}^{2}+k_{0}^{2})}\left({\rm Re}\{\chi^{(3)}_{\omega}\}+\frac{k_{0}}{n_{0}}{\rm Im}\{\chi^{(3)}_{\omega}\}\right),$		(10a)
	$\displaystyle\beta$	$\displaystyle=\frac{3\omega}{2\epsilon_{0}c^{2}(n_{0}^{2}+k_{0}^{2})}\left({\rm Im}\{\chi^{(3)}_{\omega}\}-\frac{k_{0}}{n_{0}}{\rm Re}\{\chi^{(3)}_{\omega}\}\right),$		(10b)

where $n_{0}={\rm Re}\{(1+\chi^{(1)}_{\omega})^{1/2}\}$ and $k_{0}={\rm Im}\{(1+\chi^{(1)}_{\omega})^{1/2}\}$ are the linear refractive and absorptive index, respectively, $\epsilon_{0}$ is the vacuum permittivity, and $c$ is the speed of light.

The authors thank Eduardo J. C. Dias for insightful and enjoyable discussions. J. D. C. is a Sapere Aude research leader supported by Independent Research Fund Denmark (grant no. 0165-00051B). The Center for Polariton-driven Light–Matter Interactions (POLIMA) is funded by the Danish National Research Foundation (Project No. DNRF165). Part of the computation in this project was performed on the DeiC Large Memory HPC system managed by the eScience Center at the University of Southern Denmark.

J. D. C. conceived the project and developed the theory. L. J. refined the idea and performed the numerical calculations. All authors contributed to the analysis of the results and writing of the manuscript.