Multiscale numerical tool for studying nonlinear dynamics in solids induced by strong laser pulses

Strong-field phenomena in solids exhibit extreme high-order nonlinear optical effects, which have triggered many theoretical and experimental investigations. However, there is still a lack of highly efficient numerical tools to simulate the relevant phenomena. In this paper, a versatile multiscale numerical tool set is developed for studying high-order nonlinear optical effects in solids, generated by ultrafast strong laser pulses. This tool is based on the tight-binding model approximation of the crystal structure, the related parameters of which are obtained from the density functional theory calculations. And the nonlinear effects are explored by solving the Maxwell equations coupled with the semiconductor Bloch equations. Our numerical tool can provide not only basic electronic structures and optical responses of the crystal, but also the real-time evolution of the macroscopic electromagnetic fields and the current density. The high-performance parallel computing and the interpolation method in our tool make it possible to study the strong-field nonlinear responses and propagation effects on a large spatial and temporal scale. Finally, three theoretical or experimental results published recently are satisfactorily reproduced, showing a good performance of the current toolbox.

Figure 1

Schematic illustration of multiscale calculation with a light propagation through a material thin film. The macroscopic spatial coordinate is expressed by $Z$ and the black dots represent the macroscopic grid for solving the Maxwell equation. Every macroscopic grid point in the material thin film contains a unit cell of crystal, as shown by the blue dashed lines.

Figure 2

The numerical frameworks in our tool and in the salmon package. The blue curved arrows denote the numerical processes in our tool while the green curved arrows denote the processes in salmon. And the gray curved arrows denote the relationships with the investigation of nonlinear optics in strong field. See the text for more details.

Figure 3

The electronic structures of the graphene obtained from our tool. (a) Top view of the monolayer graphene. The parallelogram in magenta is the unit cell and there are two different types of carbon atom (namely, A and B) in the unit cell. $E_0$ is the on-site energy of carbon atom, $E_1$ the nearest-neighbor hopping energy, and $E_2$ the second-nearest-neighbor hopping energy. (b) The reciprocal unit cell of the graphene (the parallelogram in dashed magenta). (c) The energy band along the path $\mathrm{\Gamma }\text{−}\mathrm{K}\text{−}\mathrm{M}\text{−}\mathrm{\Gamma }$. (d) The energy band in the plane of the monolayer.

Figure 4

(a) Comparisons of the ellipticity dependence of the seventh harmonics between the results from our numerical calculation and those from the experiment in Ref. [48], where $x$ denotes the horizontal direction of Fig. 3 while $y$ denotes the vertical direction. The circles and squares are the experimental data in Ref. [48] and solid lines are the numerical data in our calculations. For comparison, the numerical harmonic yield is normalized to the maximum value of the experimental data. (b) Ellipticity dependence of the seventh harmonic from the contributions of the intraband (blue curve) and the interband (green curve) currents. The harmonic yield is normalized to the intraband harmonic yield at the ellipticity $\epsilon =0$.

Figure 5

(a) The face-centered-cubic (fcc) geometry structure of the crystalline silicon in real space (top) and its Wigner-Seitz cell in the reciprocal space (bottom). (b) The energy band structures of silicon along the path $\mathrm{W}\text{−}\mathrm{L}\text{−}\mathrm{\Gamma }\text{−}\mathrm{X}\text{−}\mathrm{W}\text{−}\mathrm{K}$ from the TBM calculation (the blue dashed lines) and the DFT calculation (the red solid lines).

Figure 6

The run time (wall time) for parallel computing under different number of cores in Sec. 3b (see the circles). The red line shows the curve fitted by the inverse proportion function which corresponds to the ideal computing efficiency. The transverse axis of this figure represents the number of cores and the longitudinal axis corresponds to the run time (wall time).

Figure 7

(a) The time evolution of the incident wave, reflected wave, and transmitted wave (electric field) of a silicon thin film. (b) The time evolution of the energy of the electromagnetic field at the back surface macroscopic grid point of a silicon film (see Fig. 1). (c) The HHG spectra from the contributions of the transmitted wave (blue line) and reflected wave (black line). (d) The absorbing energy at the back surface macroscopic grid point of a silicon film (see Fig. 1).

Figure 8

(a) The unit cell of $\alpha$- quartz crystal in real space (top) and the projection of the Brillouin zone along the $z$ direction for a $z$-cut $\alpha$- quartz crystal (bottom). (b) The energy band structures of the silicon along the path $\mathrm{\Gamma }\text{−}\mathrm{M}\text{−}\mathrm{K}\text{−}\mathrm{\Gamma }$ from the TBM calculation (the blue dashed lines) and the DFT calculation (the red solid lines).

Figure 9

The HHG spectra from our numerical simulation and the experimental results in Ref. [55]. The blue line at the top is our theoretical HHG spectrum while the red line is that presented in Ref. [55]. To compare with the experimental result, the first three orders of harmonics are not shown in the figure. The laser parameters in both theoretical and experimental results are consistent; that is, a central carrier wavelength $\lambda =800$ nm, a peak intensity $I_0=3\times {10}^{12} \mathrm{W}/{\mathrm{cm}}^2$, and the pulse duration $\tau =30$ fs. The dephasing time is ignored in our numerical calculations.

Figure 10

The polarization dependence of HHG spectra. (a) The polarization-dependent HHG spectra with the incident wave polarized along the $\mathrm{\Gamma }\text{−}\mathrm{M}$ direction [see the bottom of Fig. 8]. (b) The polarization-dependent HHG spectra with the incident wave polarized along the $\mathrm{\Gamma }\text{−}\mathrm{K}$ direction. In both (a) and (b), $\theta =0$ corresponds to the polarization direction of the incident wave.