Effect of Strong Electric Fields on Material Responses: The Bloch Oscillation Resonance in High Field Conductivities

In this paper, we investigated the effect of strong electric fields on material responses and the Bloch oscillation resonance in high field conductivities. For this purpose, a high-order accurate explicit modal discontinuous Galerkin (DG) solver is employed for solving the quantum Boltzmann transport equation (BTE) in the context of electron transport at nanoscales under strongly out-of-equilibrium conditions. Here, we study the transient behavior and the convergence of a steady-state response to an external oscillating electric field switched on at time zero. We first benchmark our numerical results with known analytic steady-state responses at low fields. The computational results show that the present DG scheme is in excellent agreement with analytic solutions over the whole range of parameters and to an extremely high precision, allowing us to achieve good agreement even for the fifth-order response at low fields. We then extend the method to strong electric fields and show how the responses are deviated from the low-field ones and the transition to a dampened Bloch oscillation regime. Most importantly, we report the observation of a new regime induced by the resonance between the standard low-field response and Bloch oscillations.


Introduction
Transport phenomena in solids are the focus of an extremely wide range of research topics. The transport of quasi-particles (for example, electrons, phonons, or more complex quasiparticles like excitons) in nanotechnology, plasmonics, quantum computer technology, and microelectromechanical systems (MEMS) is of great importance due to their tremendous technological and scientific applications [1][2][3]. Transport phenomena are intrinsically out-of-equilibrium. However, until recently, for most applications, treating near equilibrium configurations has been sufficient. The rise of the field of ultrafast dynamics has brought forward the necessity of going beyond common approximations used close to equilibrium [4][5][6][7][8][9][10][11][12]. However, this importantly increases the complexity of the treatment as a large degrees of freedom value is no longer constrained. The unfreezing of these degrees of freedom becomes critical as qualitatively new phenomena emerge from the ensuing complexity. It therefore becomes essential to develop theoretical approaches that can fully handle the arising complexity. In the early years of the semiconductor industry, macroscopic models such as the drift-diffusion model or the hydrodynamic model have been sufficient for device simulation. Conversely, accurate simulations of modern nanoscale devices within the field of ultrafast dynamics require the use of more precise models.
The quantum Boltzmann equation has been widely proven to have the power and flexibility to describe a very wide range of dynamics of quasi-particle excitations in solids. It allows researchers to describe transport and thermalization even in heterogeneous systems. However, its application has been limited to approximated cases. The most famous is the calculation of conductivities in uniform systems, as a perturbation of the equilibrium state. However, when solved far from equilibrium, the Boltzmann equation becomes quickly numerically challenging. Here, we apply a relatively new numerical approach, discontinuous Galerkin, to the Boltzmann equation. Such method has been shown to have a number of technical and numerical advantages, such as hp-adaptivity, conservative, stable, robust with strong mathematical supports, weak approximation of boundary conditions, and built-in parallelism which permits coarse-grain parallelization [13][14][15][16][17][18]. Here, we first analyze its performance on physically relevant quantities and highlight the deep insight it provides. Critically, we also show that the method allows for the accurate description and prediction of strongly-out-of-equilibrium dynamics without compromising its precision close to equilibrium: we show that it achieves an extremely high precision in constructing high-order low-field responses.
Finally, we use the method to highlight the evolution of the responses as the system moves further from equilibrium. We show how the material response can be analyzed in terms of conductivity and nonlinear conductivities even in the presence of increasing electric fields. We show how the responses are altered further from equilibrium. The key finding is that a clear electric field amplitude-dependent resonance appears in the conductivity and nonlinear conductivities. We identify its origin as arising from the resonance between a more standard low-field driven response and dampened Bloch oscillations.

Boltzmann Transport Equation for Electrons
The time evolution of quasi-particle excitations in materials is described by a distribution function f (x, k, t), where x ∈ 3 is the position in real space, k ∈ 3 is the momentum vector, and t > 0 is the time. The domain of k is the first Brillouin zone (BZ). This distribution function fulfills the Boltzmann transport equation (BTE) [19,20], where ∂ f ∂t coll represents the Boltzmann collision integral of the interaction between two quasiparticles.
With the prescribed initial condition f (x, k, 0) = f 0 (x, k), Equation (1) must be solved for phase space coordinates (x, k) ∈ Ω x × Ω k = Ω and time t ∈ 0, t f . The quantities v x and v k appearing in Equation (1) represent the velocities in real and momentum spaces (due to external forces) which are derived from single-particle Schrodinger equation and given as where the contributions of anomalous velocity is neglected [21]. Here, e is the absolute value of electron charge,h is the reduced Planck constant, (k) is the energy dispersion relation, E(x, t) is the macroscopic electric field, and B(x, t) is the magnetic field. The right-hand side of Equation (1) is the collision term, or collision operator, representing all the scattering processes of the electrons.
In the present study, we want to (a) benchmark the numerical method for physically relevant observables and (b) observe the behaviour of the system as it is driven progressively further and further from equilibrium. For this purpose, we address only the one-dimensional BTE in the absence of magnetic field, i.e., ∂ f ∂x = 0, and B = 0. Then, In the present work, we are concerned with testing the capability of the numerical technique to produce accurate predictions for transport properties and investigate the effect of strongly out-of-equilibrium configurations on transport. To be able to have analytical benchmarks close to equilibrium, we adopt a simple collisional model: the so-called relaxation time approximation (RTA) for the collisional operator which is defined as Here, f eq (k) is a static distribution function which describes a local equilibrium, while the parameter τ is the relaxation time which characterizes the speed of recovery of the equilibrium state. For electrons, the f eq (k) is the Fermi-Dirac distribution defined as with β = 1/k B T being the inverse temperature, k B the Boltzmann constant, T the absolute temperature, and µ the chemical potential. Concluding, the BTE (1) reduces to In the present study, the electric field E is considered a time-dependent function defined as E(t) = E 0 sin(ωt), where E 0 and ω are the amplitude and frequency of the electric field, respectively. The electric field is switched on at time t = 0, triggering a transient behavior, before the system reaches a steady state.
For the numerical simulations, we consider a simple tight-binding dispersion relation, (k) = cos(2πk), ∀k ∈ Ω k = [0, L] which is shown in Figure 1a. The numerical method has the full flexibility to handle a generic band structure; however, we opted to analyze this simple case, since it is familiar to the vast majority of the readers and yet it allows for the appearance of a number of nontrivial effects compared, for instance, to a parabolic band structure. The associated Fermi-Dirac distribution function is plotted in Figure 1b.

Numerical Method Based on Modal Discontinuous Galerkin Approach
In order to discretize Equation (6), the domain Ω k = [0, L] is partitioned into N uniform and nonoverlapping elements as follows: The grid is defined by I n = [k n−1/2 , k n+1/2 ], ∀n ∈ 1, 2, · · · N, where k n = 1 2 (k n−1/2 + k n+1/2 ) is the center point and ∆k = k n+1/2 − k n−1/2 denotes the mesh size of the element I n . For the domain Ω k , we introduce the piecewise polynomial space of the functions v : Ω k −→ as where L 2 (Ω k ) is the Lebesgue space of square integrable functions over the domain Ω k and P l (I n ) is the space of polynomial functions of degree at most l in element I n . The numerical solution of Equation (6) is approximated by f h ∈ V l h in the local element I n , wheref i h is the local degree of freedom and ϕ i (k) is the basis function. The number of basis functions N p depends on the order of approximation l and is given by N p = l + 1. In the present study, the orthogonal scaled Legendre (modal) basis function is used for function ϕ i (k) in the reference element I n = [−1, 1].
where P 0,0 n (ξ) is the Legendre polynomial and k j is the centre of element. Subsequently, the BTE (6) is multiplied with the test function, which is taken to be equal to the basis function ϕ n (k), and then integrated by parts over an element I n . This results in the following weak formulation of Equation (6) where The ∂I n denotes the boundaries of the element I n and n is the outward unit normal vector. The flux function n · F appearing in Equation (11) is substituted by a numerical flux function. Here, the local Lax-Friedrichs flux F LLF is considered [22]: where Here, the superscripts (+) and (−) respectively denote the inside and outside of an elemental interface.
Among them, the quadrature theory is widely used in the discontinuous Galerkin (DG) community because of several desirable properties such as positivity of weights and symmetry. Quadrature methods are also optimal for integrating polynomials in one dimension. In the present study, the volume and boundary integrations appearing in Equation (11) are computed by Gauss-Legendre quadrature rule with (2l + 1) quadrature points to ensure accuracy [14,15].
By assembling all of the elemental contributions together, the semi-discrete DG formulation for the BTE system (6) yields an ordinary differential equation (ODE) in time for each element as Here, M −1 is the inverse of orthogonal mass matrix and L( f h ) is the residual function. In our present work, an explicit time scheme of the solution is performed with high-order strong-stability-preserving (SSP) Runge-Kutta methods that preserve the monotonicity of the spatial discretization in any norm or seminorm coupled with first-order forward Euler time stepping [24]. The explicit third-order accurate SSP Runge-Kutta method proposed by Shu and Osher [25] is employed, where L( f n h ) represents a numerical approximation of the solution at time t n , and the time step ∆t is chosen according to the following relation [15]: Here, CFL is the Courant-Friedrichs-Lewy condition (CFL ≤ 1).

Analytical Solution for the Steady-State BTE
In this section, for the sake of completeness, we report the analytical solution for the steady-state BTE response in the low-field regime. We will use this regime to benchmark our numerical results. In case of small electric field, the BTE with the relaxation time approximation (6) can be expanded in orders of the perturbation. Neglecting the electron's transient behavior, the BTE solution at low electric field can be written as the imaginary part of where ω is the electric field frequency and each term at a generic order N can be calculated as When the BTE dynamics are known, important collective quantities, such as the current density, can be calculated as From Equation (17), one can obtain the current at the appropriate frequency and consequently the complex conductivity as well as all the nonlinear corrections to it (in the presence of time reversal symmetry, the even order currents vanish), Thus, the current from Equation (18) can be expressed as

Convergence Study
In order to verify the code and estimate the order of accuracy of the numerical scheme, we first solve a simple problem for which an analytical solution exists. Here, we examine the accuracy of the present numerical scheme with P 1 , P 2 , P 3 , and P 4 elements. The L 2 and L ∞ errors are measured in terms of the discrete norms by wheref (k, t) is the exact solution and f h (k, t) is the approximate solution. The error(N 1 ) and error(N 2 ) are the errors (for L 2 , and L ∞ ) for the numbers of elements N and 2N, respectively. Our convergence study was done for the following initial value problem, with the periodic boundary condition f (0, t) = f (2π, t). The exact solution of this problem is given by f (k, t) = sin (k − t). This problem was calculated up to 4th-order of accuracy at time t f = 1 by varying the number of grid points N(20, 40, 80, 160). The results of the convergence study with the order of accuracy L 2 and L ∞ are illustrated in Table 1. From the results, it was confirmed that the present numerical scheme achieved the desired order of accuracy (l + 1).

Validation of Numerical Solver
An analytic solution cannot be evaluated when the full BTE with a time-dependent external field and the relaxation time approximation term is used. In order to verify the reliability and accuracy of the present numerical solver, we compute the different orders of electrical conductivity for small electric fields with the low-field analytical expressions for the high-order conductivities obtained at the steady state, as discussed in Section 4. We let the simulations run long enough for the steady state to be reached and then apply the Fourier transform to the current temporal profile over the last period. We then compare the Fourier components. We computed the numerical and analytical solutions for electric field, E 0 = 2 × 10 2 V/cm and at β = 50, τ = 1.0, µ = −0.8, with a wide range of frequency value ω = 0.01 − 20. Figure 2 shows the numerical and analytical conductivity for first-, third-, and fifth-order (higher order responses have been multiplied by a power of the amplitude of the electric field to have the same dimensionality of a conductivity, and then all the responses were normalized to the zero frequency conductivity for a more direct comparison). The computed results for real, imaginary, and magnitude parts of conductivity are found to be in excellent agreement with the analytical solutions, even for the high-order responses. Discrepancy arises when the response becomes too small compared to the numerical error. This shows how this method, that describes the population within the full Brillouin zone, can achieve an extremely high precision, even in the low-field regime, where the absolute modification of the population is extremely small and restricted to very limited momentum regions.
Finally, in Figure 3, we show the dependence of the conductivity on the relaxation time. Consistently with Equation (18), the relaxation time controls both the amplitude of the conductivity and the frequency at which the drop in conductivity is observed.

Numerical Results and Discussion
In this section, we present the numerical results for the BTE dynamics in regimes that range from equilibrium to strongly nonequilibrium. Emphasis is placed on the time evolution of the distribution function, including the effects of physical parameters and steady-state electrical conductivity. In this work, we express the parameters in the following units: volts per centimeter (V/cm) for the electric field; one over nanometer (1/nm) for the momentum; picosecond (= 10 −12 s) for the time; Terahertz (= 10 12 hertz) for the frequency; and picosecond (= 10 −12 s) for the relaxation time. The constant value of the ratio e/h is 1.54 × 10 −4 . In general, the solutions of the BTE system are largely affected by the following flow parameters: the electric fields E, chemical potential µ, inverse temperature β, and relaxation time (τ). These flow parameters are chosen to demonstrate the effects of nonequilibrium situation on the BTE dynamics. In order to investigate these nonequilibrium effects, the electric field values E 0 range from 1∼10 6 V/cm. The chemical potential µ is selected to range from −1.0 to 1.0, while the inverse temperature β ranges from 1 (smooth distribution) to 10 2 (sharp distribution). Three different relaxation time parameters are chosen for the extensive studies, including τ = 0.01, 0.1, and 1.0. A wide range of frequency values (ω) has been chosen from 0.01 − 20. In all the numerical simulations, a polynomial expansion of third-order accuracy is used to approximate solutions in the finite element space. All the computations are carried out on 300 grid points, and the periodic boundary conditions are enforced at the right and left boundaries. The initial condition for the distribution function is considered as Fermi-Dirac function.

Time Evolution of BTE Dynamics: Effects of Flow Parameters
Flow parameters-the electric field, the inverse temperature, and the chemical potential-play a critical role in ultrafast nonequilibrium electron dynamics. Therefore, we conducted a detailed investigation into the effects of these Flow parameters on BTE dynamics. To demonstrate the effects of the electric field on BTE dynamics, four different electric field values, E 0 = 1, 10 3 , 10 4 , and 5 × 10 4 V/cm, are selected with the same τ = 1.0, β = 50, µ = 0, and ω = 1. Figure 4a illustrates the time evolution of the electronic distribution function for different electric fields.
At a low electric field E 0 = 1 V/cm, the distribution function remains close to the equilibrium state, as shown in Figure 4a. We stress that the shift of the Fermi surface (not shown), although very small, is properly reproduced by the numerics, as showed in the previous section. The steady state is reached after a time of the order of the relaxation time. The transient shows a relatively simple behavior with the distribution progressively approaching the steady-state response.
As the electric field value increases to E 0 = 10 3 V/cm, the steady-state distribution is still near to the equilibrium distribution, as shown in Figure 4a. However, the population is not a rigid shift of the Fermi surface anymore. Even at low temperatures, the Fermi edge, on top of being shifted, is also smoothened. The system is already beyond the close-to-equilibrium regime, and the electric field cannot be considered a small perturbation anymore. The transient time does not look significantly altered by the strength of the electric field.
When the electric field value increases to E 0 = 10 4 V/cm, the system enters into the dampened Bloch oscillation regime. The distribution function undergoes strong accelerations, which push electrons across the border of the Brillouin zone triggering Bloch oscillations. In Figure 4a, one can observe a sizable population traversing the border of the Brillouin zone. Due to the presence of the scattering operator however, the sharp profile of the electrons' distribution is eventually dampened, and the system reaches a far-from-equilibrium steady state. At very high electric field value E 0 = 5 × 10 4 V/cm, the distribution undergoes faster crossing of the Brillouin zone before reaching a steady state, as shown in Figure 4a. The motion is not fully resolved in the figure since only a fraction of the computed time steps have been plotted.
To analyze the effects of the inverse temperature (β = 1/k B T) on the BTE dynamics, four different inverse temperature values, β = 1, 5, 50 and 10 2 , are selected with the same electric field value E 0 = 10 4 V/cm with µ = 0, τ = 1, and ω = 1. We focus on the most interesting regime: the dampened Bloch oscillation regime shown for different inverse temperatures in Figure 4b. At higher temperatures, the electronic distribution is smooth, as shown in Figure 4b. At low temperatures, the distribution displays sharp drops at the Fermi level (in principle, it becomes discontinuous at T = 0). Notice how, in spite of the almost discontinuous solution, the numerical method proved to perform well and remained stable even at very low temperature, as seen in Figure 4b.
Interestingly, one can clearly distinguish the transition between two regimes. At very high temperature, the steady-state distribution shape is strongly smoothened due to the temperature. On the other hand, at very low temperature, the steady-state solution becomes completely dominated by the electric field: as can be observed in Figure 4b, the steady-state distributions are practically independent of the initial temperature.  Further, we show the results for several positions of the chemical potential for a strong electric field. Figure 4c illustrates the effect of the chemical potential µ on the time evolution of the distribution function with E 0 = 10 4 V/cm, τ = 1.0, β = 50, and ω = 1.0. Here, we selected different chemical potential values between −1.0 to 1.0, as shown in Figure 4c. The pictures show the transition from an electron regime to a hole regime.

High Field Current Density and Conductivities: The Bloch Oscillation Resonance
Here, we analyze in more detail the effect of high electric fields onto the material responses, focusing on the effect on the induced current. Figure 5 illustrates the effects of electric fields on the time evolution of oscillating current density with τ = 1.0, β = 50, µ = 0, and ω = 1.0.
We first focus on the transient behavior. At a small electric field E 0 = 1 V/cm, the current density reaches the steady-state response after a transient of the order of few relaxation times, as shown in Figure 5a. Interestingly, during the transient, the current density value first grows above the steady-state response amplitude. This behavior is consistently observed for all the electric fields' amplitudes.

Frequency Dependence: The Bloch Oscillation Resonance
We now address the steady-state response, which is reached after the transient response. It can be observed that there are significant differences in time evolution of the current density with increasing electric field value. The first evident effect is the appearance of response at higher harmonics. This is expected even in the small perturbation regime, since the responses of the system at higher harmonics are proportional to higher powers of the field amplitude, and become more and more relevant at higher fields. However, there is a further effect playing an important role: the proportionality constants of the lower harmonics (conductivity) are importantly altered in high fields. Figure 6 shows the frequency-dependence of the first-(conductivity), the third-, and the fifth-order responses at increasing electric fields (all the values are normalized to the low-field conductivity at zero frequency). At low fields (line corresponding to E 0 = 10 3 V/cm) the first order response is simply the conductivity reported in Figure 2a. However as the electric field is increased, the response of the material shows high field effects, that strongly depend on the frequency of the external field. At low frequencies, the conductivity decreases (see Figure 6a) as it enters the dampened Bloch oscillation regime: a large electric field can accelerate electrons across the border of the Brillouin zone, where the electron motion inverts direction, depressing the conductivity. Conversely, the reduction in the conductivity is instead not observed at high frequencies: given the high frequency of the oscillating frequency, even in the presence of a high field, electrons still remain relatively close to the Fermi surface and keep responding in the linear regime.
The most interesting regime is between the abovementioned ones. It is easier to understand the physics behind this regime by first analyzing the dynamics of an electron in a static electric field. As an electron is accelerated across the Brillouin zone due to a static electric field, its velocity changes in time, generating a time-dependent contribution to the current which oscillates with a Bloch oscillation period where ∆K is the size of the Brillouin zone. The reader should note how the Bloch oscillation frequency increases linearly with the electric field's amplitude. However, at low electric fields, such response can hardly be observed since such oscillations are dampened by the scattering term. Therefore, only at sufficiently high electric fields, when the oscillation period is sufficiently shorter than the relaxation time, can these dynamics become observable. In case the external electric field oscillates, a similar dynamics is triggered. However, it is important to realize that, when doing the Fourier transform of the current, it will produce contributions at its resonance frequency and not necessarily at the frequency of the applied field. However, if the applied electric field happens to be at a frequency close or equal to the Bloch oscillation frequency, a resonance appears in the conductivity. This is clearly observed in Figure 2a, where no resonance is observed at low fields; while for increasing fields, the Bloch resonance is observed to both increase in strength and shift to higher frequencies.
The effect appears in the higher-order conductivities as well, but it can be observed to be triggered at lower external field's frequencies. This is indeed understandable since, as mentioned, the Bloch oscillation frequency is dependent only on the amplitude of the external field and always gives a response at that frequency. When the electric field's frequency is one third of the Bloch oscillation frequency, the third-order response in the current happens to be exactly resonant with the Bloch oscillations, giving the peak observed in Figure 2b. In the fifth-order response in Figure 6c, one can observe how the resonance is triggered at even lower frequency, which is one fifth of the Bloch oscillation frequency or the frequency at which the resonance is observed in conductivity.
As a further comment, let us highlight that in the higher-order responses, if the electric field is small, the overall high-frequency component of the current becomes too small and the computed response simply shows the numerical error.

Chemical Potential and Temperature Dependences
Finally, we analyze how the high-field conductivity depends on the other parameters of the system. The dependence on the chemical potential in Figure 7 does not show importantly new effects, as the conductivity remains proportional to the number of free carriers (whether electrons or holes) within the band, even at high fields. The fact that low-frequency responses are more strongly suppressed at higher fields is also clearly observed.

Concluding Remarks
This study focused on an investigation of the strong electric field's effect on material responses and the Bloch oscillation resonance in high-field conductivities. For this purpose, a high-order accurate explicit modal discontinuous Galerkin is applied to study the transport quantum Boltzmann equation far from equilibrium. We first focused on the analysis of performance of the presented method by testing its precision to close-to-equilibrium conditions. We have compared the method to analytic results for the low-field steady-state regime and showed that it achieves extremely high precision, even if the method is devised for strong excitations and not optimized for extremely small deviations from equilibrium. Secondly, we have performed an extensive range of numerical simulations to investigate the effects of strong electric fields, and the dependence of the response on a range of physical parameters. In particular, we have shown how corrections to conductivity arise at high fields. This allowed us to identify a new regime, appearing at high electric fields, triggered by the resonance of typical material responses and dampened Bloch oscillations. Experimental confirmation of the appearance of such resonance requires the choice of a material with a relatively long scattering lifetime (for example, low-temperature, high-purity silicon) in order to allow the effect to be present at not too high frequencies where different contributions to conductivity start playing important roles.