1. Introduction

Light propagation at oblique incidence is one of the fundamental subjects in macroscopic electromagnetism. In linear optics, the reflection, refraction, and transmission of $s$- and $p$-polarized light at a flat surface are described by Fresnel’s equations and this is used as a basic method to measure optical constants. An incidence at Brewster’s angle is used to avoid the reflection in many optical systems [1]. A Kretschmann configuration is used for the generation of evanescent fields and surface plasmon polaritons [2]. In nonlinear optics, the control of incident angle is essential to achieve phase matching [3]. Recently, measurements of highly nonlinear and ultrafast electronic dynamics in solids have been attracting interest [4–7]. In such experiments, measurements are often carried out at oblique incidence. For example, high-harmonic generation from solids at oblique incidence provide useful information on the mechanism and opportunity to control the spectra [8].

Although the theoretical description of light propagation at oblique incidence on a flat surface is well documented in textbooks of macroscopic electromagnetism [1], it becomes a research topic when the light-matter interaction accompanies nonlinear and/or nonlocal effects. For example, for oblique-incidence irradiation on periodic nano-structures such as metamaterials or metasurfaces, several theoretical methods have been developed recently including a field transformation method [9], a split-field method [10–12] and an iterative method [13]. For nonlinear light propagation, constitutive relations with nonlinear optical susceptibilities of second and third orders may be used when the nonlinear effect is sufficiently weak and can be treated perturbatively [14]. However, for a strong field that does not allow perturbative description, it is necessary to solve coupled equations of macroscopic light propagation and microscopic dynamics of electrons and ions.

Recently, multiscale methods have been developed for nonlinear light propagation that cannot be described using constitutive relations. In one such method, macroscopic Maxwell equations for light propagation and microscopic quantum-mechanical equations for electronic dynamics are solved simultaneously. For example, the microscopic electronic dynamics are described by the quantum Liouville equation for the density matrix in the Maxwell-Bloch theory [15], and those by the time-dependent Kohn-Sham (TDKS) equation for the electronic orbitals in the Maxwell-TDDFT (time-dependent density functional theory) [16]. Since these methods do not use any constitutive relations, they are applicable to a propagation with extreme nonlinearity such as high-harmonic generation from thin films [17], the ultrafast modulation of optical properties by intense pulsed light [5], and initial stages of laser processing [18]. However, multiscale calculations have mostly been carried out for the one-dimensional propagation of a plane electromagnetic field at normal incidence where the light propagation can be described by a one-dimensional wave equation. Since microscopic calculations of electron dynamics are computationally expensive, it has been difficult to apply the method to problems with multi-dimensional light propagation.

This paper is intended as a report of a computational method to describe light propagation at oblique incidence in the multiscale Maxwell-TDDFT method. In this method, the equation to describe the light propagation at oblique incidence is transformed into a one-dimensional wave equation. Therefore, the computational cost is similar to that in the case of the light propagation at normal incidence. Because of the discontinuity of the vector potential at the surface, however, careful implementation is required to obtain an accurate solution. We note that a method employing the one-dimensional wave-propagation equation has been developed for a propagation in linear media at oblique incidence [19]. In the field of plasma physics where the microscopic dynamics are described by classical theory such as the Boltzmann equation, a method to describe nonlinear light propagation at oblique incidence has also been developed [20].

The paper is organized as follows. In Sec. 2, we derive a one-dimensional equation that will be used to calculate the light propagation at oblique incidence. Then a numerical implementation of the multiscale Maxwell-TDDFT method is explained. In Sec. 3, we show some results of oblique-incidence calculations using the multiscale Maxwell-TDDFT method, taking the irradiations of a flat silicon surface as an example. Finally, a summary is given in Sec. 4.

2. Theory

2.1 Multiscale Maxwell-TDDFT method

We first explain the formalism of the multiscale Maxwell-TDDFT method [16] for the case of general three-dimensional light propagation. In the method, electromagnetic fields that describe light propagation are expressed using a vector potential, ${\boldsymbol A}({\boldsymbol R},t)$, where we denote the coordinate as ${\boldsymbol R}$ for macroscopic quantities. The Maxwell equation is given by

To relate ${\boldsymbol J}$ with ${\boldsymbol A}$, a local and linear constitutive relation is assumed in ordinary macroscopic electromagnetism,

In practical calculations, we introduce a Bloch orbital $u_{{\boldsymbol R}, b{\boldsymbol k}}({\boldsymbol r},t)$ to describe the microscopic electronic motion at ${\boldsymbol R}$. Here, $b$ is the band index and ${\boldsymbol k}$ is the Bloch wavenumber. A coordinate ${\boldsymbol r}$ is introduced to describe microscopic electronic motion. The Bloch orbital satisfies the following TDKS equation;

From the Bloch orbitals, the electric current density at ${\boldsymbol R}$ is obtained as

Equations (3) and (4) without the coordinate index ${\boldsymbol R}$ are known as the basic equations of real-time TDDFT in crystalline solids and have been utilized for many purposes such as dielectric functions [23], nonlinear susceptibilities [24], high-harmonic generation [25], and energy transfer from a light pulse to electrons [26].

Solving Eqs. (3), (4), and (1) simultaneously, we can describe macroscopic light propagation as well as microscopic electron dynamics. In practical calculations, we introduce a grid system in the macroscopic coordinate ${\boldsymbol R}$ to solve Eq. (1). The Bloch orbital $u_{{\boldsymbol R}, b{\boldsymbol k}}({\boldsymbol r},t)$ is prepared and the TDKS Eq. (3) is solved at each grid point ${\boldsymbol R}$ inside the medium. To start the time-evolution calculation, the Bloch orbital at every grid point is set to the ground state in the static density functional theory, while an incident pulse is set in the vector potential in the vacuum area. We note that, when the electric field is sufficiently weak, the electric current given by Eq. (4) is related to the vector potential that is used in Eq. (3) by the linear constitutive relation of Eq. (2). Therefore, the multiscale Maxwell-TDDFT method results into an ordinary macroscopic Maxwell equation with the linear constitutive relation given by the linear response TDDFT [23].

Since it is necessary to solve the TDKS Eq. (3) at all grid points of ${\boldsymbol R}$ in the medium simultaneously, it becomes a large-scale calculation even for one-dimensional light propagation. Because of the high computational cost, applications of the multiscale Maxwell-TDDFT method have so far been limited to cases of pulsed light irradiated normally on flat surfaces or thin films for which the Maxwell Eq. (3) is a one-dimensional wave equation [5,6,18,27,28].

2.2 Maxwell equations at oblique incidence

We next show that the Maxwell equations for the light propagation at oblique incidence can be expressed as a one-dimensional wave equation. We consider the irradiation of a flat surface by a plane wave of pulsed light with the incident angle of $\theta$, as shown in Fig. 1(a). We express the macroscopic coordinate as ${\boldsymbol R} = (X,Y,Z)$. Uniform matter is in the $Z>0$ region with the surface at the $Z=0$ plane. For the incident plane wave, we choose the $XZ$-plane as the incident plane. The electromagnetic fields are uniform in the $Y$-direction so that the vector potential is independent of the $Y$-coordinate.

 

Fig. 1. (a) Setup of the system and coordinates for the incident angle $\theta$. (b) The incident electric field is shown in the $XZ$-plane. (c) Electric field for the cross section along the $Z$-axis.

Download Full Size | PDF

From the translational symmetry, we find that the vector potential and electric current density have the following dependence on the coordinates $X$, $Z$, and time $t$:

Using ${\boldsymbol a}$ and ${\boldsymbol j}$, the Maxwell Eq. (1) is written for each Cartesian components as follows. For $a_Y$, we have

For $a_X$ and $a_Z$, we have the following coupled equations:

In general, Eqs. (7), (8), and (9) are coupled and must be solved simultaneously. For an isotropic linear medium, Eq. (7) describes the propagation of $s$-polarized light, while Eqs. (8) and (9) describe the propagation of $p$-polarized light.

For the coupled equations of (8) and (9), we further rewrite them as following. First, we integrate Eq. (9) over time to obtain

Differentiating Eq. (10) with respect to $Z$ and combining it with Eq. (8), we have

We will use this equation to calculate the time evolution of $a_X$. Once $a_X$ is obtained, $a_Z$ may be obtained by integrating Eq. (10) over time. Equations. (7), (10), and (12), which are partial differential equations of $Z$ and $t$, are the basic equations to be solved numerically to calculate the light propagation at oblique incidence.

2.3 Boundary conditions

To carry out numerical calculations, it is important to clarify the behavior of the vector potential at the interface, $Z=0$. In this subsection, we will examine the boundary conditions for each component of the vector potential. To confirm that the boundary conditions are consistent with the ordinary results, we examine the reflection and transmission through the surface for an isotropic and linear medium.

We first examine the boundary condition for $a_Y$. In Eq. (7), there appear $a_Y$ and $j_Y$. The current $j_Y$ is discontinuous at $Z=0$. Accordingly, $(\partial ^{2}/\partial Z^{2})a_Y$ must also be discontinuous at $Z=0$. This fact indicates that both $a_Y$ and $(\partial /\partial Z) a_Y$ are continuous at $Z=0$. We consider the consequence of these boundary conditions for an incident pulse of $s$-polarization on an isotropic and linear medium characterized by the dielectric constant $\epsilon$. The following form is assumed for the solution:

The coefficients $R_s$ and $T_s$ are determined by the boundary condition that both $a_Y$ and $(\partial /\partial Z)a_Y$ be continuous at $Z=0$. Then, we have

From these relations, we obtain the well-known result for $s$-polarization,

We next consider boundary conditions for $a_X$ and $a_Z$ components. In Eq. (12), $(\partial /\partial Z) p_Z$ on the right-hand side is divergent at $Z=0$ since $p_Z$ is discontinuous. We integrate both sides of Eq. (12) over a small region of $Z$ including $Z=0$, $[-\delta, +\delta ]$ with an infinitesimally small $\delta$. Then, we have

From this relation, we find that $(\partial /\partial Z) a_X$ is discontinuous and, then, $a_X$ is continuous at $Z=0$.

To find the boundary condition for $a_Z$, we evaluate Eq. (19) at $Z=\pm \delta$ and combine it with Eq. (10) to obtain

From this relation, we find that $a_Z$ is discontinuous at $Z=0$. We note that Eq. (20) is the ordinary continuity condition for the perpendicular component of the electric displacement, $D_Z = E_Z + 4\pi P_Z$.

We consider the case of an incident pulse of $p$-polarization on an isotropic and linear medium. The solution can be expressed as

To determine the coefficients $R_p$ and $T_p$, we utilize the continuity of $a_X$ and Eq. (20). They provide

From these relations, we obtain the well-known result for $p$-polarization,

2.4 Numerical implementation

To achieve a multiscale Maxwell-TDDFT calculation at oblique incidence, we solve the one-dimensional propagation Eqs. (7), (10), and (12) together with the microscopic Eqs. (3) and (4). To solve the one-dimensional propagation equations, we introduce a uniform spatial grid $\{ Z_i \}$ in the $Z$ coordinate. We solve the propagation equations using finite difference approximation and do not use the boundary conditions at $Z=0$. As discussed in the previous subsection, there appears a function $(\partial / \partial Z) p_Z$ that diverges at $Z=0$ on the right-hand side of Eq. (12), and such divergence is cancelled by a divergence of the $(\partial ^{2} /\partial Z^{2}) a_X$ term. We numerically found that the divergent term induces unphysical oscillations if we naively implement the finite difference approximation.

To avoid this numerical difficulty, we have found it useful to introduce smearing in the vicinity of the boundary, $[Z_{\rm vac}, Z_{\rm mat}]$, as shown in Fig. 2. In the smearing region, we multiply a weight function $w$ to the electric current density,

The weight function $w$ is a smooth function, changing from 0 in the vacuum region to 1 in the matter region. In our implementation, we will use the following function that is defined in $[Z_{\rm vac}, Z_{\rm mat}]$:

 

Fig. 2. Spatial grid systems used in the multiscale Maxwell-TDDFT method. The smearing region is indicated by gradation.

Download Full Size | PDF

To update the vector potential, we adopt the following scheme:

In Eqs. (27),(28), $1/(\cos \theta )^{2}$ becomes divergent as $\theta$ approaches $90^{\circ }$ (grazing incidence), and the method does not work in the limit. This is already evident in the 1D transformed Eqs. (7)-(9), which are singular as $\theta$ approaches $90^{\circ }$.

To verify that the smearing method is effective, we performed a test calculation of pulse propagation at oblique incidence on a bulk surface of a linear medium. We use an incident pulse with the central frequency $\omega _c = 1.55$ eV. A dielectric function of a classical Drude-Lorentz model,

In Fig. 3, we show the reflectance of $p$-polarized light as a function of the incident angle $\theta$. In the calculation, a grid spacing of $\Delta {Z}=0.53$ nm is used. The results of calculations using different numbers of grid points $N_s$ for smearing, $N_s=8$ (orange dotted curve) and $N_s=4$ (blue solid curve), are compared with the results of exact Fresnel formula of reflection (green broken curve) [1]:

The numerical results are in good agreement with the exact value. The minimal reflectance at the Brewster angle of $\theta \approx 75^{\circ }$ is accurately described by the numerical calculation with the smearing. A comparing calculations with $N_s=4$ and $8$, the accuracy is improved by decreasing the smearing points $N_s$. This indicates that, as expected, the smearing region $N_s \Delta Z$ should be small to obtain on accurate result. However, we find that a finer grid spacing $\Delta Z$ is required to obtain a stable numerical solution when a small $N_s$ value is adopted.

 

Fig. 3. Reflectance of $p$-polarized light from a surface of a linear medium as a function of the incident angle $\theta$. Numerical results with different smearing points, $N_s=4$ and 8, are compared with the exact value.

Download Full Size | PDF

Let us discuss the limiting behavior of the smearing method when we use a small smearing area and a small grid spacing. For simplicity, we restrict ourselves to the linear regime in which our Maxwell-TDDFT formalism reduces to ordinary macroscopic electromagnetism with a dielectric function in TDDFT. In our smearing method, it is possible to improve the description as much as we want by simultaneously decreasing the length of the smearing area and also the grid spacing. In the limit, it is possible to reproduce exact results of boundary-value problem at an interface in macroscopic electromagnetism. Therefore, using our smearing method with sufficiently small smearing area and grid spacing, it is possible to describe phenomena that can be described in ordinary macroscopic electromagnetism such as the Fresnel reflection, as confirmed above, and the surface plasmon at the metal-dielectric interface.

3. Illustrative calculations

In this section, we show several calculations of light propagations at oblique incidence on a silicon thin film of 50 nm thickness using the multiscale Maxwell-TDDFT method. For the incident pulse, we use a pulse with a cos-square-shaped envelope with the central frequency of $\hbar \omega = 1.55$ eV and the full duration of $T_p = 20$ fs. We set the $Z=0$ plane as the $(001)$ surface of silicon. A linear $p$-polarization is assumed for the incident pulse with the incident plane parallel to $(010)$. Because of the symmetry, there is no electric current density in the $Y$-direction. Therefore, only two propagation Eqs. (10) and (12) are utilized to describe the propagation.

We have implemented the light propagation at oblique incidence in the SALMON code [30] that has been developed by the authors’ group. In the code, a real-space 3D grid is used to express electron orbitals. A norm-conserving pseudopotential is used for the electron-ion interaction [31,32]. Numerical parameters are chosen as follows. For the propagation of the light pulse, we use a spatial grid size of $\Delta Z = 1.0$ nm so that the number of grid points in the Si film is $N=50$. We choose the number of smearing points $N_s$=8 for the front and back surfaces. For the electron dynamics calculation, we use a cubic unit cell containing eight silicon atoms. Each side of the cell is divided into 16 grid points. The $k$-space is sampled by $8 \times 8 \times 8$ Monkhorst-Pack grids. A common time step is used for the light propagation and electron dynamics calculations, $\Delta t = 0.43$ as.

We use the adiabatic local density approximation (ALDA) as the exchange-correlation potential [33]. It is well known that the LDA systematically underestimates the optical gap energy of dielectrics, which results in the overestimation of the index of refraction $n$. In the case of bulk Si, ALDA gives $n=4.0$ at the frequency of $1.55$ eV, while the measured value is $n \simeq 3.7$ [29].

Figure 4(a) shows snapshots of the two components of the electric field, $E_X$ (red) and $E_Z$ (green). The thin film is placed in the $Z=0$ region and appears as a line without thickness at this scale. The incident angle is chosen at $60^{\circ }$, which is smaller than Brewster’s angle. The maximum intensity of the pulse is set at $I=10^{9}$ W/cm$^{2}$. At $t = 0$ fs (top), the incident pulse is placed at the left of the thin film and propagates in the positive-$Z$ direction. At $t=10$ fs (middle), the pulse reaches the film and interacts with it. At $t=20$ fs (bottom), the pulse is separated into reflected and transmitted waves. Figure 4(b) shows expanded views of the electric field in the vicinity of the thin film. It can be seen that $E_X$ is continuous at the surfaces of the film, while $E_Z$ is discontinuous. This behavior is consistent with the boundary condition discussed in Sec. 2.3.

 

Fig. 4. (a) The electric field in the multiscale Maxwell-TDDFT calculation is shown at $t=0$ fs (top), $t=10$ fs (middle), and $t=20$ fs (bottom). $E_X$ is shown by the red line and $E_Z$ by the green line. The silicon thin film is located within $0 < Z < 0.05$ $\mu$m. (b) The expanded view of the electric field in the thin film region is also shown.

Download Full Size | PDF

Figure 5 shows two-dimensional views of the electric field intensity, $\vert {\boldsymbol E}(X,Z,t) \vert ^{2}$, at a certain time $t$ where the electric field is generated by ${\boldsymbol E}(X,Z,t) = -(1/c) (\partial /\partial t) {\boldsymbol a}(Z, t-X\sin \theta /c)$. For the cases of Fig. 5(a) is the case of a weak pulse ($I=10^{9}$ W/cm$^{2}$) at the incident angle of 60$^{\circ }$, (b) a strong pulse ($I=10^{13}$ W/cm$^{2}$) at the incident angle of 60$^{\circ }$, (c) a weak pulse ($I=10^{9}$ W/cm$^{2}$) at the incident angle of 76$^{\circ }$, and (d) a strong pulse ($I=10^{13}$ W/cm$^{2}$) at the incident angle of 76$^{\circ }$. We note that the incident angle of 76$^{\circ }$ corresponds to the Brewster’s angle of Si with the index of refraction $n=4.0$ in the ALDA calculation.

 

Fig. 5. Electric field intensity, $|E(X,Z,t)|^{2}$ at a fixed time $t$ calculated by the multiscale Maxwell-TDDFT method. Two incident angles of $\theta = 60^{\circ }$ and $76^{\circ }$ (Brewster’s angle) , and two incident intensities of $I=10^{9}$ W/cm$^{2}$ (weak) and $10^{13}$ W/cm$^{2}$ (strong) are shown in panels (a)-(d). The magnitude of the field intensity is normalized by the maximum value of the incident pulse.

Download Full Size | PDF

In Figs. 5(a) and (b), there appear both transmitted and reflected waves. We find that the transmitted intensity is much lower in (b) than in (a). Since we set the central frequency of the incident pulse at $\hbar \omega = 1.55$ eV, which is below the direct bandgap of Si, no absorption takes place at the weak intensity shown in Fig. 5(a). At the higher intensity shown in Fig. 5(b), on the other hand, substantial absorption takes place at the film by multiphoton absorption processes, making the transmission weaker at the high intensity.

In Fig. 5(c), the reflected wave cannot be seen and the transmitted wave has almost the same magnitude as the incident wave, as expected for the incidence of the weak pulse without absorption at Brewster’s angle. Figure 5(d) also shows for incidence at the Brewster’s angle but with an increased pulse intensity, $I = 10^{13}$ W/cm$^{2}$. Here, the magnitude of the transmitted wave is smaller than that of the incident pulse owing to the multiphoton absorption, as in the case shown in Fig. 5(b). Looking carefully at Fig. 5(d), an interference pattern between the incident and reflected waves can be seen in the area surrounded by the dashed-line circle. It indicates that that is slight reflection at Brewster’s angle owing to nonlinear light-matter interaction for the strong incident pulse.

To gain a deeper understanding of the propagation, waveforms of the reflected and transmitted pulses are shown in Fig. 6 for the strong incident pulse of $I = 10^{13}$ W/cm$^{2}$ at the incident angles of 60$^{\circ }$ (blue) and 76$^{\circ }$ (orange), as well as waveforms for the weak pulse of $I=10^{9}$ W/cm$^{2}$ (black-dashed curve). The latter is multiplied by a factor of 100 to show it in the same scale as the others. At the angle of 60$^{\circ }$, the transmitted pulse looks similar between the strong and weak cases. On the other hand, the reflected pulse in the strong case gradually becomes weaker with time. This indicates that the optical property of the medium gradually changes with the continued irradiation of the strong pulse. At Brewster’s angle of 76$^{\circ }$, the reflected wave does not appear in the weak case, as expected. However, a substantial reflection appears in the strong case. The waveform of the reflected wave also shows angled structures, indicating that it includes substantial nonlinear components. For the transmitted wave, on the other hand, the amount of modulation looks rather small. We thus find that nonlinear effects appear more clearly in the pulse shape of the reflected wave than in that of the transmitted wave.

 

Fig. 6. Reflected (left) and transmitted (right) pulse shapes in the multiscale Maxwell-TDDFT method. Blue and orange curves correspond to the incident angles of $\theta =60^{\circ }$ and 76$^{\circ }$, respectively, for a strong incident pulse of $I=10^{13}$ W/cm$^{2}$. Black-dashed lines are for a weak pulse of $I=10^{9}$ W/cm$^{2}$ and are multiplied by a factor of 100.

Download Full Size | PDF

Figure 7 shows power spectra of the reflected and transmitted waves. We find that the spectra include 3rd and 5th harmonic generations as well as the fundamental wave. In the spectrum of the reflected wave, the fundamental wave component at Brewster’s angle is one order of magnitude smaller than that at $60^{\circ }$. However, the magnitude of the nonlinear components is similar between $60^{\circ }$ and $76^{\circ }$. This finding explains why the reflected wave of the strong incident pulse at the Brewster’s angle shows angled structures. In the spectrum of the transmitted pulse, it is found that the nonlinear component at Brewster’s angle is smaller than that at $60^{\circ }$. This finding indicates that nonlinear effects show complex dependence on the incident angle, and the usefulness of the present approach when investigating them.

 

Fig. 7. Spectral intensities of reflected (left) and transmitted (right) waves for incident angles of $\theta =60^{\circ }$ (blue) and 76$^{\circ }$ (orange).

Download Full Size | PDF

Finally, in Fig. 8(a), we show a reflectance for incident pulses of various intensity and angles in the multiscale Maxwell-TDDFT calculation. For a weak pulse, the reflectance shows a clear minimum with zero reflectance at Brewster’s angle. As the intensities increases, the minimum value of the reflectance gradually increases. The angle at the minimum reflectance also decreases. In Fig. 8(b), we show reflectance using the ordinary Fresnel formula with a complex index of refraction. The value of the index of refraction at each intensity is determined by fitting the result shown in Fig. 8(a). The fitted values of the index of refraction are indicated in the figure. As seen from the figure, the behavior of the reflectance calculated by the multiscale Maxwell-TDDFT method can be reasonably fitted by the complex index of refraction. The imaginary part of the index of refraction determined by the fitting procedure gradually increases as the pulse intensity increases, while the real part decreases.

 

Fig. 8. Angular dependence of the reflectance of a 50-nm-thick Si thin film. (a) Multiscale Maxwell+TDDFT calculation for four incident intensities from $10^{9}$ to $10^{13}$ W/cm$^{2}$. (b) Fresnel reflection using a complex dielectric constant determined to fit the reflectance of the multiscale Maxwell-TDDFT calculation.

Download Full Size | PDF

The change in the index of refraction may be understood in terms of the increase of the carrier density by the multiphoton excitation process. In the TDDFT calculation with a pump-probe set up, it has been shown that the excited carriers behave as free carriers and it adds the Drude-like term $-\sqrt {n_\mathrm {ex}}/(\omega ^{2}+i\gamma \omega )$ to the dielectric function $\epsilon$ [34] where $\gamma$ is a decay constant. As the pulse intensity increases, the carrier density $n_\mathrm {ex}$ increases. Then the real part of the index of refraction $n_r$ ($=\sqrt {\epsilon }$) decreases, and the imaginary part increases. These behaviors of the reflectance are consistent with the results shown in Fig. 8(a).

4. Summary

We presented a theoretical formalism and numerical method to calculate nonlinear light propagation at oblique incidence based on time-dependent density functional theory (TDDFT), a first-principles computational method for electronic motion, coupled with the macroscopic Maxwell equations. In our method, microscopic electronic motion is calculated by solving the time-dependent Kohn-Sham equation, the basic equation of TDDFT, at each grid point that is used to solve the macroscopic Maxwell equations. Since the TDDFT calculation is computationally expensive, the number of grid points should be kept as small as possible. In our implementation, we express the Maxwell equations at oblique incidence as a one-dimensional wave equation for the vector potential and employ one-dimensional grid points to solve it. In this way, the computational cost of the light propagation calculation at oblique incidence is similar to that of the light propagation at normal incidence. Since the normal component of the vector potential is discontinuous at the interface between the medium and the vacuum, we is effective proposed a smearing method to carry out stable calculations and confirmed that it is effective.

As a demonstration of the method, we calculated the light propagation at oblique incidence on a silicon thin film. It was shown that the multiscale Maxwell-TDDFT calculation adequately describes the linear optical response such as the absence of the reflected wave at Brewster’s angle. For a strong incident pulse, the absence of reflectance at Brewster’s angle is lifted by the nonlinear optical response. It was shown that the angle dependence of the reflectance for strong incident pulses can be reasonably described by the nonlinear change in the index of refraction that is caused by multiphoton excitation.