1. Introduction

Metallic nanoparticles were generated in a vacuum in the mid-1980s, and since then, their structural and optical properties have been extensively investigated [1]. The quantum mechanical effects in metallic nanoparticles and their size dependence have been explored theoretically [2]. Subsequently, metallic nanostructures fabricated on a substrate became a typical platform for plasmonic systems, which enabled a new way to delicately control the light-matter interactions utilized in various applications, such as optoelectronic devices, solar cells, and biosensors [3–6]. Recently, quantum mechanical effects in the optical response of nanostructures have attracted renewed interest and have been demonstrated both theoretically and experimentally. Specifically, when the size of the nanostructures is as small as several nanometers, the optical response is affected by spatial nonlocality [7–13] and electron spill-out [14–20]. Furthermore, state-of-the-art fabrication techniques have enabled nanostructures with angstrom-level fineness separated by subnanometer gaps between plasmonic particles [21–26]. At such an extremely microscopic scale, electrons in the gaps exhibit quantum tunneling [21–23,27–34]. These quantum effects have attracted attention for enhancing optical nonlinearity [8,14,20,32,34], which is important for developing and improving a wide range of applications [35–40]. In this context, there is a growing interest in establishing numerical schemes to analyze the optical nonlinearity caused by quantum effects in metallic nanostructures.

The most straightforward way to accommodate this demand is via a first-principles calculation that employs electronic orbitals, such as time-dependent density functional theory (TDDFT) [41,42]. Although the feasibility of first-principles calculations is limited to small systems, using the TDDFT with a jellium model (JM) that can capture important features of the quantum effects may be applicable to nanosystems comprising up to a few thousand electrons [2]. Therefore, the combination of TDDFT and the JM has been extensively used in the last decade to analyze the optical properties of isolated nanodimers (two metallic nanoparticles) with a subnanometer gap, and this combination could accurately reproduce experimental results that include electron tunneling through the gap [21,22]. However, numerical calculations based on TDDFT with the JM for three-dimensional nanostructures is significantly time-consuming because its computational cost is proportional to the product of the spatial size of the system and the number of orbitals. Therefore, even if the calculation is executed on cutting-edge supercomputers, the nanostructure length that can be calculated within a reasonable computational time is limited to approximately 10 nm. However, the typical size of nanostructures in real measurements is considerably larger—several tens or hundreds of nanometers in most cases [3–6].

A semiclassical approach using kinetic energy (KE) functionals was developed to overcome the computational restrictions of the TDDFT with the JM, and it has recently gained interest in the field of nanoplasmonics. In this approach, the electron dynamics is treated as a fluid, the motion of which is described by local physical quantities, such as the electron density $n({{\textbf r},t} )$ and electric current density ${\textbf J}({{\textbf r},t} )$ or velocity field ${\textbf v}({{\textbf r},t} )$. The semiclassical approach no longer involves electronic orbitals; therefore, it enables calculations that include large nanostructures with an acceptable computational cost. Such an approach with KE functionals has a long history in various fields [43–46]. In the last decade, this approach has attracted attention in the field of nanoplasmonics, taking into account the quantum mechanical effects. Specifically, the spatial nonlocality of nanostructures was first described using the Thomas-Fermi (TF) KE functional, which depends on n [7–12]. Subsequently, the von Weizsäcker (vW) KE functional, which has a $\nabla n$ dependence, has been incorporated to address electron spill-out and quantum tunneling [14–20]. This is referred to as $\; \textrm{TF}\lambda \textrm{vW}$, and it includes a phenomenological parameter $\lambda $, the value of which is adjusted to reproduce the electron spill-out in the TDDFT calculation. Very recently, higher-order KE functionals, including the ${\nabla ^2}n$ contribution, have been introduced to more precisely describe electron dynamics [47]. In previous literature, semiclassical approaches with and beyond $\textrm{TF}\lambda \textrm{vW}$ are referred to as quantum hydrodynamic theory (QHT). Notably, in one study, QHT was successfully derived from TDDFT, which clarifies their mutual correspondence [19].

QHT is a promising avenue to analyze the quantum mechanical effects in the optical response of metallic nanostructures; however, numerical calculations based on QHT have thus far been limited to linear optical responses, and they are hardly applied to nonlinear optical phenomena. This is owing to the computational difficulty of the KE functionals, in which $1/n$ and/or the cube root of n are involved. Therefore, even if the density becomes zero or negative ($n \le 0$) at a single spatial point during time evolution calculations, owing to numerical errors in the update of n, the QHT calculation immediately breaks down. It was difficult to develop a numerical scheme that exactly assures $n > 0$ in the whole spatial region and every time step; therefore, almost all previous studies using QHT have focused on the linear response of metallic nanostructures, where a linear equation for the change in the density from the ground state is calculated, whereas the update of the density is unnecessary [15–19]. To the best of our knowledge, a few studies have addressed the nonlinear optical response using QHT with $\textrm{TF}\lambda \textrm{vW}$ (called QHT$- \textrm{TF}\lambda \textrm{vW}$ hereafter), and these studies revealed the mechanism of second harmonic generation originating from the electron spill-out at the interface of a metallic thin film with [14] and without [20] the static approximation. However, in their numerical implementation, a perturbative expansion was adopted to describe the second-order optical nonlinearity. Therefore, it is necessary to develop a new code to analyze higher-order nonlinearities that involve growing complexity as the order of the nonlinearity increases. Furthermore, owing to the perturbative treatment, their calculation was compelled to adopt an undepleted pump approximation, in which the fundamental field is not affected by the nonlinear process; however, it becomes important in some nonlinear optical applications [35–37].

Time-dependent orbital-free density functional theory (TD-OFDFT) is an alternative choice to reduce the computational cost of TDDFT calculations [48–53]. It is derived by mapping a real system of interacting fermions onto a fictitious system of noninteracting bosons. In fact, two methods, QHT and TD-OFDFT, are intimately related with each other [53]. In TD-OFDFT calculations, the numerical issue to keep the density positive-definite during time evolution calculations can be handled easily since the density is defined by absolute square of the orbital. However, a potential of a form of $\Delta \sqrt n /\sqrt n \; $ usually appears in the Hamiltonian. Numerical calculation of such term may involve difficulties at low density region, although highly accurate evaluation is required, in particular, in nonlinear response calculations. TD-OFDFT have been applied mostly to linear response properties. In [48,52], it was applied to nonlinear phenomena, dipole response of a small nanocluster and stopping power of warm dense matter, respectively.

In the present study, a numerical scheme is developed to stably calculate both linear and nonlinear electron dynamics based on QHT$- \textrm{TF}\lambda \textrm{vW}$. The key component of our approach is to rewrite the QHT equation using the effective Schrödinger equation (ESE), which ensures a stable calculation of the KE functionals without any numerical instabilities. This ESE involves only one orbital; therefore, the computational cost is considerably lower than that of the TDDFT with the JM. The potential in the ESE does not include a term of a form of $\Delta \sqrt n /\sqrt n $. The time evolution of the orbital can be solved using, for example, the real-time and -space method, which has been well tested in previous studies. We couple the ESE with Maxwell’s equations including effects of magnetic field that is important for nanoobjects of large size. A method to introduce damping is proposed using a fictitious lossy medium. As a demonstration, we present the calculations of the linear and third-order nonlinear optical responses of a metallic nanosphere, and we compare the results to those of the TDDFT with the JM. We show that the two calculations are in reasonable agreement with each other.

2. Theoretical framework

2.1 Summary of TDDFT and QHT formalisms

In this subsection, we present the basic formulas of TDDFT and QHT for which the numerical methods are developed. A direct derivation of QHT from TDDFT is presented based on a description given in the literature [19].

We consider an isolated system, in which a metallic nanostructure is subjected to an electromagnetic field. In TDDFT, the electron dynamics is described by the following time-dependent Kohn-Sham (KS) equation for electronic orbitals ${\psi _j}({{\textbf r},t} )$:

The time evolutions of ${\textbf A}$ and $\phi $ are governed by the following inhomogeneous wave equations with n and ${\textbf J}$:

In QHT, we describe the electron dynamics using the equation for the electric current density ${\textbf J}({{\textbf r},t} )$. In the case of QHT$- \textrm{TF}\lambda \textrm{vW}$, which is used in this study, it is given by [19]

Here, the first and second terms on the right-hand side are from the electromagnetic fields, ${\textbf E}$ and ${\textbf B}$, respectively, which are equivalent to the Lorentz force. The third term, which is composed of n and ${\textbf J}$, is the convective term, and the fourth term comes from the energy functionals. ${{\textbf A}_{\textrm{XC}}}({{\textbf r},t} )$ and $\delta {E_{\textrm{XC}}}/\delta n = {V_{\textrm{XC}}}({{\textbf r},t} )$ are the XC vector and scalar potentials, respectively, which have already appeared in the TDKS equation, given in Eq. (1). ${T_{\textrm{TF}}}$ is the TF KE functional, and its functional derivative is given as:

In Eq. (6), an adjustable phenomenological parameter $\lambda $ is included. This value depends on the derivation of Eq. (6), which ranges from $1/9$ to $1$ [46]. We will later determine the value from the condition that the electron spill-out is reasonably described.

The history of the fluid dynamical description of fermionic many-body systems originates from the TF theory in 1927 [43,44]. The vW correction term in Eq. (8) was first derived in 1935 [45], which showed that $\lambda = 1$. Then, an expression with $\lambda = 1/9$ was derived using the density expansion method [46]. Notably, Ciracì recently provided a direct derivation of Eq. (6) from the TDDFT [19]. They have also provided high-order KE functionals systematically [47]; however, these are beyond the scope of the present study.

The linear optical response based on QHT$- \textrm{TF}\lambda \textrm{vW}$ can be executed by linearizing Eqs. (6–8), where it is unnecessary to explicitly calculate the time evolution of electron density. However, to extend the application of this method to a nonlinear regime, an evolution of the electron density $n({{\textbf r},t} )$ is required. This is typically achieved by numerically employing the equation of continuity, $\frac{{\partial n}}{{\partial t}} = \frac{1}{e}\nabla \cdot {\textbf J}.$ However, the time evolution calculation of the electron density is often accompanied by numerical difficulty. Density is a positive definite quantity; however, it is not easy to implement a numerical method that strictly ensures it. Additionally, an accurate calculation of the density is difficult in spatial regions where the density is extremely low. The $\delta {T_\textrm{w}}/\delta n$ term includes $1/n$; therefore, an inaccuracy in the density calculation that causes $n \le 0$ even at a single spatial point immediately breaks down the whole calculation. Similarly, the calculation of $\delta {T_{\textrm{TF}}}/\delta n$ also induces challenges because it includes the cube root of n, which has multiple branches for $n < 0$. Thus, it is considerably difficult to realize a numerical calculation by implementing Eqs. (6–8) with an update of the electron density.

2.2 QHT solved using ESE

In 1926, Madelung proposed that the Schrödinger equation for a single-particle problem can be transformed into the form of hydrodynamic equations [54]. In 1988, an idea to overcome the above difficulty of accurately calculating the density in the fluid dynamical approach was presented in the field of nuclear physics utilizing Madelung’s correspondence, namely, converting the fluid dynamical equation into an equivalent Schrödinger equation [55]. We apply this idea to QHT$- \textrm{TF}\lambda \textrm{vW}$, and Eqs. (6–8) are equivalent to the following ESE with $\lambda = {\xi ^2}$:

The derivation of Eq. (6) from Eqs. (9–11) is written in Supplement 1.

Thus, a positive definite density can be obtained after the wave function $\mathrm{\Psi }({{\textbf r},t} )$ is obtained by solving Eq. (9). Accurate methods have been developed to solve the Schrödinger equation, even in a spatial region where the density is extremely low; therefore, we may overcome the numerical difficulty of the accurate calculation of the density. It is noteworthy that the term involving $\delta {T_\textrm{w}}/\delta n$ is transformed into the KE operator in the ESE in Eq. (9). This means that it is unnecessary to calculate quantities that include the $1/n$ term.

Previously, the ESE was used to calculate the ground state in the QHT [18]. However, to the best of our knowledge, the ESE has never been applied to the time-dependent electronic dynamics described by QHT$- \textrm{TF}\lambda \textrm{vW}$. We also note that the previous application of the ESE in nuclear physics was limited to dynamics described by a rotation-free velocity field [55]. However, our implementation includes velocity fields with rotation because we treat the Lorentz force caused by a magnetic field.

2.3 Optical system interacting with the ESE

In our practical implementation, we solve the following four equations for the electromagnetic fields in the Lorentz gauge instead of Eqs. (4) and (5):

To solve Eqs. (12) and (13), we can make full use of numerical methods that are well developed in computational electromagnetic analysis. Therefore, the proposed QHT scheme simultaneously solves Eqs. (9–15). In the next section, we investigate the problem displayed in Fig. 1(a), in which pulsed light with a linear polarization irradiates a nanosphere described by the JM. For this problem, the calculation proceeds as follows. The propagation of electromagnetic fields ${\textbf E}$ and ${\textbf B}$, including incident pulsed light, is described by Maxwell’s Eqs. (12) and (13). The electron dynamics described by the effective wave function $\mathrm{\Psi }$ is updated using Eq. (9) with the electromagnetic potentials ${\textbf A}$ and $\; \phi $, which are given by Eqs. (14) and (15). When the incident pulsed light that is emitted from the electric current source (see Supplement 1) reaches the nanosphere and excites the electrons, the excited electrons generate an electric current density ${\textbf J}$, as defined by Eq. (11). This, in turn, goes into Eq. (12) of Maxwell’s equations as a source term. The feedback from the excited electrons to the optical system generates a new electromagnetic field that interacts with the electronic system again and acts as the origin of the optical near field and plasmon resonance.

The numerical formulas to solve Eqs. (9–15) in discretized real-time and -space grids are described in detail in Supplement 1.

 

Fig. 1. (a) Schematic of the studied Na nanosphere with the diameter a. ${{\textbf E}^{(i )}}$ and ${{\textbf B}^{(i )}}$ are the electromagnetic fields of the incident light propagating toward the negative z direction. The center of the nanosphere is located at the origin. Electron density distribution in the ground states calculated using the (b) DFT and (c) QHT on the $xy$ plane including the origin. Their amplitudes are normalized by ${n_s}$. (d) Comparison of the normalized electron densities ${n_0}/{n_s}$ along the $x$-axis. The blue solid line with circles (red broken line with triangles) shows the density of the DFT (QHT), and the black dotted line with squares denotes the background positive charge density of the JM, $n_{\textrm{JM}}^{(+ )}$. The inset shows a comparison in a log-scale. The phenomenological parameter $\lambda $ in QHT is set to $1/2$.

Download Full Size | PDF

2.4 Numerical aspects

We compare numerical aspects of TDDFT of Eq. (1), original QHT of Eq. (6), and ESE of Eq. (9). In Table 1, we summarize scaling behavior of computational costs of these methods in terms of the numbers of orbitals ${N_o}$, spatial grids ${N_r}$, and time steps ${N_t}$. We note that a common computational cost of electromagnetic analysis given by Eqs. (12–15), which scales ${N_r}{N_t}$, should be added in each theory. Denoting the system size as N (e.g., the number of electrons in the system), the number of orbitals in the TDDFT calculation is proportional to the system size, ${N_o} \propto N$. The number of spatial grids also usually proportional to the system size, ${N_r} \propto N$. Therefore, assuming that a common and size-independent time step is used, the computational cost scales $O({{N^2}} )$ in TDDFT, while $O(N )$ in the original QHT of Eq. (6) and the ESE of Eq. (9). We thus find that the computational cost of ESE is much less than TDDFT, and is comparable to the original QHT.

Table 1. Comparison of computational costs

View Table

3. Illustrative calculations

In the previous section, we explained the proposed numerical scheme used to stably execute the time evolution calculations of electric current density and electromagnetic fields based on QHT$- \textrm{TF}\lambda \textrm{vW}$ using the ESE. In this section, we present illustrative calculations of the linear and third-order nonlinear optical responses of a metallic nanosphere using this method. The results calculated based on QHT$- \textrm{TF}\lambda \textrm{vW}$ using the ESE are compared with those obtained using TDDFT with the JM, which reveals the reliability and limitations of the QHT.

3.1 Setup for the calculation of a nanosphere

Figure 1(a) displays the studied system, in which a metallic nanosphere with a diameter a is subjected to incident light that is specified by the electromagnetic fields ${{\textbf E}^{(i )}}$ and ${{\textbf B}^{(i )}}$. The nanosphere is described by the JM, which replaces an ionic structure with a positive background density $n_{\textrm{JM}}^{(+ )}({\textbf r} )$ of a spherical shape with the following shape boundary [2]:

3.2 Ground state

Before exploring the optical responses, we calculate the initial electron density ${n_0}({\textbf r} )$ in the ground state, and we use the ESE for the ground state calculation. The time-independent ESE for the effective wave function ${\mathrm{\Psi }_0}({\textbf r} )$ in QHT$- \textrm{TF}\lambda \textrm{vW}$ is given by:

To solve Eqs. (17) and (18) numerically, we discretize a spatial area using a three-dimensional Cartesian grid. The spatial grid spacings $\mathrm{\Delta }x$, $\mathrm{\Delta }y$, and $\mathrm{\Delta }z$ are set to 0.05 nm. This value is chosen so that the electron spill-out at the subnanometer scale can be accurately described, as will be demonstrated later. The computational domain is set to a 14.6 nm cube, which is substantially larger than the diameter of the nanosphere a. A finite difference formula is used for the spatial derivatives, and an imaginary time propagation method [57] is used to calculate the ground state self-consistently. These methods are explained in Supplement 1. The boundary value of the static scalar potential ${\phi _0}$ is determined by the multipolar expansion.

Figure 1(b) and 1(c) display the spatial distribution of the ground-state electron density ${n_0}$ calculated using DFT and QHT, respectively. We set the value of the phenomenological parameter $\lambda $ as $\lambda = 1/2$ that well reproduces the electron spill-out in the DFT calculation, as discussed below. Figure 1(c) exhibits a lack of ring-like patterns. These patterns originate from the Friedel oscillation, which is characterized by the Fermi wavenumber, and the quantum mechanical shell effect formed by the electronic orbitals ${\psi _j}$ interfered with each other. Therefore, QHT, which is categorized as a semiclassical approach without the orbitals, cannot reproduce the density oscillation. For a clearer comparison, we show their electron density distributions along the $x$-axis in Fig. 1(d), where the inset shows a comparison in a log-scale. The blue broken line with circles and the red broken line with triangles display the electron density of the DFT and QHT, respectively. As a reference, we also plotted the background positive density of the JM, $n_{\textrm{JM}}^{(+ )}$, using a black dotted line with squares. Although QHT does not reproduce the oscillation caused by the quantum effect that appears at $x \le 2$ nm, the density tail at $x\; >\;2$ nm, including the electron spill-out, exhibits a good agreement between QHT and DFT. This agreement is due to an appropriate choice of the $\lambda $ value, which critically affects the density tail and electron spill-out. In the following subsections, we use $\lambda = 1/2$ in the real-time calculations for linear and nonlinear optical responses.

3.3 Linear response

In this subsection, we numerically solve the QHT using the ESE, given in Eqs. (9–11), coupled with the electromagnetic fields described by Eqs. (12–15) employing real-time and -space grids. The linear response of the metallic nanosphere shown in Fig. 1 is calculated using the dipole approximation.

As stated in previous studies [15–20,47], QHT cannot describe the damping effect that originates from electron-electron interactions among the electronic orbitals ${\psi _j}$. Therefore, a friction term to model the damping effect is added to Eq. (6) with the phenomenological damping rate parameter $\gamma $. However, in QHT using the ESE, it is difficult to introduce such a friction term in the Schrödinger equation. Thus, instead of the friction term, we introduce a background lossy medium that adds the following conductive current density ${{\textbf J}_C}({{\textbf r},t} )$ in Eq. (12):

We adopt a dipole approximation in the linear response calculation because the diameter of the sphere is sufficiently smaller than the light wavelength. As a perturbation, we use the following temporarily impulsive and spatially uniform field with polarization in the $x$-direction:

The linear optical absorption cross-section is given by

We compare the results of the QHT and TDDFT. The TDDFT calculation method is described below. We employ the length gauge in the dipole approximation. Namely, while Eq. (2) remains unchanged, Eqs. (1,3,4) are modified as follows:

Figure 2(a) shows the temporal profile of the total dipole moment ${d_x}$ (the $x$-component of ${\textbf d}$) for the nanosphere shown in Fig. 1, with a diameter of $a = 4.3$ nm and an electron number of ${N_e} = 1074$, subjected to the impulsive incident field ${{\textbf E}^{(i )}}$ of Eq. (22). In the TDDFT calculation exhibited by the blue solid line, the damping effect that originates from electron-electron interactions among the electronic orbitals ${\psi _j}$ appears. The dipole moment by the TDDFT induces damping by $t \approx 30$ fs, after which a small oscillation continues. The red broken line is the result of the QHT using the ESE with $\sigma = 0$. As described above, QHT with no loss does not include the damping effect, which results in an undamped oscillation. To reproduce damping in the TDDFT calculation, we choose the value of $\sigma $ to be ${\sigma _{\textrm{TDDFT}}} = 5.05 \times {10^3}$ S/m ($1.10 \times {10^{ - 3}}$ a.u.). The green dotted line in Fig. 2(a) shows the result obtained using QHT with this value of ${\sigma _{\textrm{TDDFT}}}$. Although this value of $\sigma $ is determined using a comparison with TDDFT in the present study, it may also be determined by referring to the experimentally observed absorption spectrum of nanoparticles.

 

Fig. 2. (a) Temporal profile of the total dipole moment ${d_x}$ (the $x$-component of ${\textbf d}$) for a nanosphere with a diameter of $a = 4.3$ nm subjected to the impulsive incident field of Eq. (22). The blue solid line shows the calculation using TDDFT, and the red broken and green dotted lines show the QHT using the ESE with $\sigma = 0$ and ${\sigma _{\textrm{TDDFT}}}$, respectively. (b) Spectral distribution of the linear optical absorption cross-section ${S_{\textrm{abs}}}$ for the same nanosphere. The blue solid line shows the TDDFT, and the other three lines (red, green, and black) show the QHT with $\sigma = {\sigma _{\textrm{TDDFT}}}/2$, ${\sigma _{\textrm{TDDFT}}}$, and $2{\sigma _{\textrm{TDDFT}}}$, respectively. ${S_{\textrm{abs}}}$ of two different nanospheres with (c) $a = 3.1$ nm and (d) $a = 5.3$ nm with $\; \sigma = $ ${\sigma _{\textrm{TDDFT}}}$. The blue and red lines show the results using the TDDFT and QHT, respectively.

Download Full Size | PDF

Figure 2(b) shows the spectral distribution of the linear optical absorption cross-section ${S_{\textrm{abs}}}$. The blue solid and green dotted lines are the TDDFT and QHT results with ${\sigma _{\textrm{TDDFT}}}$. Their maximum peaks at approximately 3 eV originate from the plasmon resonance, the amplitudes and frequencies of which exhibit good agreement between the TDDFT and QHT. However, TDDFT includes several peaks around the main peak that are not observed in QHT. These additional peaks are created from a number of electronic orbitals ${\psi _j}$. QHT cannot reproduce such peaks, which indicates its limitations. Furthermore, we observe another peak at approximately 3.8 eV only in QHT. This peak is known as the Bennett state, and it is caused by the nonuniformity of the electron density [59]. This state is known to also appear in TDDFT as a small peak or shoulder [31]; however, QHT overestimates its appearance. Recent studies have revealed that this overestimation of the Bennett state can be mitigated using more sophisticated QHTs beyond the $\textrm{TF}\lambda \textrm{vW}$ level, such as dynamical correction or high-order KE functionals [19,47], which are not considered in the present study. The red broken and black chained lines in Fig. 2(b) exhibit spectra with different choices of $\sigma $: ${\sigma _{\textrm{TDDFT}}}/2$ and $2{\sigma _{\textrm{TDDFT}}}$. These lines clearly show the importance of the choice of $\sigma $ because it controls the peak widths of the plasmon resonance and Bennett state while maintaining their resonant energies.

To illustrate the robustness of QHT using the ESE, we calculate nanospheres of different sizes: (i) $a = 3.1$ nm $({{N_e} = 398} )$; (ii) $a = 5.3$ nm $({{N_e} = 2030} )$. Figure 2(c) and (d) display their ${S_{\textrm{abs}}}$ for $a = 3.1$ and $5.3$ nm, respectively, where $\sigma $ is set to ${\sigma _{\textrm{TDDFT}}}$. In Fig. 2(c), for the $a = 3.1$ nm nanosphere, the TDDFT calculation plotted using the blue solid line shows several fragmented peaks around the plasmon resonance frequency of 3 eV, owing to a decrease in the number of occupied orbitals ${\psi _j}$ and the larger energy spacing between them. The QHT plotted using the red broken line cannot reproduce such spectra. In Fig. 2(d), for the $a = 5.3$ nm nanosphere, a good agreement is observed between the TDDFT and QHT results. In particular, the spectral widths of the plasmon resonance coincide well with each other, compared with the case of the $a = 4.3$ nm nanosphere, shown in Fig. 2(b). The plasmon in the TDDFT calculation is almost single-peaked, owing to an increase in the number of orbitals, which makes the energy levels more dense.

3.4 Nonlinear response

We now consider the calculation of the nonlinear optical response of the $a = 4.3$ nm nanosphere displayed in Fig. 1. The following incident light pulse ${{\textbf E}^{(i )}}$ polarized along the $x$-axis is employed:

Figure 3(a) shows the power spectra of the total dipole moment ${d_x}$, where the Fourier transformation is applied from $t = 0$ to ${T_{\textrm{max}}}$ with the window function $w(t )$ defined in Eq. (27). The blue solid and red broken lines are calculated by the TDDFT and QHT using the ESE, respectively. The two vertical gray dashed lines indicate the frequencies of the plasmon resonance ${\omega _{\textrm{Mie}}}$ and Bennett state ${\omega _{\textrm{Bennett}}}$. We show the calculations for three incident pulses with the following fundamental frequencies: $\hbar {\omega _i} = 0.86\; ({ \approx {\omega_{\textrm{Mie}}}/3.5} )$, $1.00\; ({ \approx \hbar {\omega_{\textrm{Mie}}}/3} )$, and $1.22\; ({ \approx \hbar {\omega_{\textrm{Mie}}}/2.5} )$ eV. The top panel of Fig. 3(a) corresponds to the pulse irradiation with $\hbar {\omega _i} = 0.86$ eV. Below 3.5 eV, we observe three peaks: the fundamental and 3^rd harmonic generation at approximately 0.86 and 2.58 eV, respectively, and the plasmon resonance at ${\omega _{\textrm{Mie}}}.$ QHT succeeded in reproducing these peaks. However, in the energy region higher than 3.5 eV, two additional peaks appear in QHT, but they are not found in TDDFT. The first peak at approximately 3.8 eV is owing to the Bennett state, which is overestimated in QHT, as demonstrated in Fig. 2(b). The second peak at approximately 4.3 eV corresponds to the 5^th harmonic generation. The 3^rd harmonic generation of QHT and TDDFT coincide accurately with each other; however, the 5^th harmonic generation is only visible in the QHT. The reason for this is currently unclear. The Bennett state that appears strongly in the QHT may affect this phenomenon. The middle panel of Fig. 3(a) corresponds to the pulse irradiation with $\hbar {\omega _i} = 1.00$ eV. The 3^rd harmonic generation of ${\omega _i}$ overlaps with the plasmon resonance ${\omega _r}$. Therefore, the 3^rd nonlinear signal is plasmonically enhanced here, and the QHT reproduces the enhanced nonlinearity well. In the case of the top panel, $\hbar {\omega _i} = 0.86$ eV, two peaks are formed by the Bennett state and 5^th harmonic generation. Similarly, the 5^th harmonic generation is observed only in the QHT. The bottom panel corresponds to the pulse irradiation with $\hbar {\omega _i} = 1.22$ eV. In this case, the 3^rd harmonic generation overlaps with the Bennett state. A large difference appears in the 3^rd nonlinear signal between the TDDFT and QHT because the Bennett state appears strongly in QHT. The 5^th harmonic generation is slightly visible at approximately 6 eV in the QHT.

 

Fig. 3. (a) Power spectra of the total dipole moment ${d_x}$ for the nanosphere with a diameter of $a = 4.3$ nm. The nanosphere is subjected to the pulsed incident light described in Eq. (31). The blue solid (red broken) lines are from the TDDFT (QHT) calculations. The two vertical gray dashed lines indicate the frequencies of the plasmon resonance ${\omega _{\textrm{Mie}}}$ and Bennett state ${\omega _{\textrm{Bennett}}}$. The top, middle, and bottom panels correspond to the pulse irradiation with $\hbar {\omega _i} = 0.86\; ({ \approx \hbar {\omega_{\textrm{Mie}}}/3.5} )$, $1.00\; ({ \approx \hbar {\omega_{\textrm{Mie}}}/3} )$, and $1.22\; ({ \approx \hbar {\omega_{\textrm{Mie}}}/2.5} )$ eV, respectively. (b) ${\omega _i}$ dependence of the third-order optical nonlinearity $d_{\textrm{NL}}^{(3 )}$ defined by Eq. (32). The blue solid line with circles (red broken line with triangles) denotes the results obtained using TDDFT (QHT).

Download Full Size | PDF

Finally, we define the following quantity to explore the difference in third-order nonlinearity between the TDDFT and QHT more quantitatively:

4. Conclusion

This study proposes a numerical scheme that enables stable calculations of linear and nonlinear electron dynamics induced by an optical field described using QHT$- \textrm{TF}\lambda \textrm{vW}$. The key component of our scheme is to rewrite the conventional equations of the QHT$- \textrm{TF}\lambda \textrm{vW}$ to the ESE, and to solve the ESE coupled with Maxwell’s equations using real-time and -space grids. In QHT described by the ESE, all numerical instabilities originating from $1/n$ and the cube root of n terms that appear in KE functionals have been removed, thereby making it easy to expand the computational targets to nonlinear optical phenomena without any approximations. We believe that the proposed numerical method of QHT using the ESE will be useful for exploring nanoplasmonics that include quantum mechanical effects.

The result was compared with that of the calculations for the fully quantum mechanical DFT with the JM for the ground state and the TDDFT calculations for the linear and nonlinear optical responses of a metallic nanosphere. We have shown that an appropriate choice of $\lambda $ in the QHT is essential to reproduce the density tail in the ground state and the plasmonic spectrum in the linear response in the (TD)DFT calculations. In particular, our QHT calculation successfully reproduced the third-order nonlinearity enhanced by plasmon resonance in the TDDFT calculation. However, at an energy region higher than the plasmon resonance, the QHT overestimated the Bennett state in the linear response, which originated from the nonuniformity of the electron density. Owing to this overestimation in the linear response, the nonlinear spectrum calculated using QHT gradually deviated from that of the TDDFT at an energy region higher than the plasmon. To remedy this overestimation, more sophisticated QHTs that rely on dynamical correction or high-order KE functionals may be necessary.