Spatio-temporal dynamics of shift current quantum pumping by femtosecond light pulse

Shift current—a photocurrent induced by light irradiating noncentrosymmetric materials in the absence of any bias voltage or built-in electric field—is one of the mechanisms of the so-called bulk photovoltaic effect. It has been traditionally described as a nonlinear optical response of a periodic solid to continuous wave light using a perturbative formula, which is linear in the intensity of light and which involves Berry connection describing the shift in the center of mass position of the Wannier wave function associated with the transition between the valence and conduction bands of the solid. Since shift current is solely due to off-diagonal elements of the nonequilibrium density matrix that encode quantum correlations, its peculiar space–time dynamics in response to femtosecond light pulse employed locally can be expected. To study such response requires to analyze realistic two-terminal devices, instead of traditional periodic solids, for which we choose paradigmatic Rice–Mele model sandwiched between two metallic electrodes and apply to it time-dependent nonequilibrium Green function algorithms scaling linearly in the number of time steps and capable of treating nonperturbative effects in the amplitude of external time-dependent fields. This reveals novel features: superballistic transport, signified by time dependence of the displacement, ∼tν with ν > 1, of the photoexcited charge carriers from the spot where the femtosecond light pulse is applied toward the electrodes; and photocurrent quadratic in light intensity at subgap frequencies of light due to two-photon absorption processes that were missed in previous perturbative analyses. Furthermore, frequency dependence of the DC component of the photocurrent reveals shift current as a realization of nonadiabatic quantum charge pumping enabled by breaking of left–right symmetry of the device structure. This demonstrates that a much wider class of systems, than the usually considered polar noncentrosymmetric bulk materials, can be exploited to generate nonzero DC component of photocurrent in response to unpolarized light and optimize shift-current-based solar cells and optoelectronic devices.


Introduction
The conventional photovoltaics is based on semiclassical transport of electrons and holes excited by light and separated by the built-in electric field within a pn-junction where careful control of disorder and flow of light inside a solar cell is required to reduce entropic energy losses [1]. An alternative is distinctly different physical mechanism-the so-called bulk photovoltaic effect (BPVE) [2]-which appears in noncentrosymmetric materials [3], such as ferroelectrics with nonzero electric polarization. The BPVE generates steady-state photocurrent and the above band-gap photovoltage as the signatures of nonlinear optical response. The BPVE can be substantially enhanced in ferroelectric thin films [4], and potentially even more in Weyl semimetals [5]. However, photocurrent generated by BPVE in a closed circuit remains small compared to pn-junctions. This has ignited renewed experimental [4,[6][7][8][9][10][11] and theoretical [12][13][14][15][16][17][18] interest in BPVE with the ultimate goal [15] to optimize its power conversion efficiency.
Aside from potential applications, the nonlinear optical response has also emerged [19,20] as a sensitive probe of local (or geometric) properties of quantum wavefunction in noncentrosymmetric materials, including topological insulators [21] and Weyl semimetals [5,22,23]. The geometric properties of wavefunction are encoded by the Berry-phase-concepts, such as Berry connection and Berry curvature, which govern DC component of injection current or shift current driven by high-energy electric fields of circularly-or linearlypolarized light, respectively. This is somewhat surprising since such fields induce interband transitions in which electronic wavefunctions are not confined to specific bands and adiabatic condition is violated [19,20].
The BPVE can have contributions from different processes [2,9,24,25], where the most intensely studied are: (i) the shift current contribution, traditionally explained as a shift in real space following the carrier interband transition [24][25][26]; and (ii) the ballistic current contribution due to intraband transition of hot photoelectrons with asymmetric momentum distribution during which they lose their energy to descend to the bottom of conduction band while shifting in real space [2]. The shift current has attracted particular attention due to its inherently quantum-mechanical nature involving phase-coherent evolution of electron and hole wave functions. This makes possible rapid propagation of charge carriers toward the electrodes, thereby minimizing carrier recombination (that reduces the magnitude of conventional photocurrent) and energy losses due to carrier-phonon scattering [25]. For example, very recent experiments [10,11] shining a laser spot onto the middle of the sample have detected shift current across the whole sample, independently of the position and the width of the excited region, and extracted it ∼100 μm away from the excited region.
The widely used [12][13][14][15][16][17] 'standard model' [25,26] for computing the magnitude of shift current is based on a perturbative expansion which yields its DC component as a second-order nonlinear optical response to electric field E q (ω) of the incident monochromatic light of frequency ω and polarization in the q-direction. The unpolarized sunlight is modeled by averaging σ pqr (ω) over different polarizations. Thus, equation (1) can only describe photocurrent linear in intensity, as typically encountered also in conventional photovoltaics [27,28]. The third-rank tensor σ pqr (ω) can be expressed as a product of two terms carrying intuitive meaning [12,15,17]-diagonal imaginary part of the dielectric function, which is proportional to the density of states (DOS); and the so-called shift vector (of magnitude [2] 10-100 nm) as the average distance traveled by coherent photoexcited charge carriers during their lifetime. The shift vector depends on the Berry connection [19,20] of the Bloch energy bands of an infinite periodic crystal and it is, therefore, an intrinsic property of the material insensitive to elastic and inelastic scattering processes [25]. However, the 'standard model' tailored for infinite periodic crystals cannot be used to compute open-circuit photovoltage, as one of the key factors determining the power conversion efficiency [14], or examine the effect of the electrodes on the measured photocurrent. Furthermore, it cannot describe photocurrent response to a femtosecond light pulse applied locally to a finite-size open quantum system attached to macroscopic reservoirs, as illustrated in figure 1. On the other hand, such setups have recently become popular in experiments [6-8, 10, 11] aiming to understand: (i) how fast shift current responds to light and how fast it propagates toward the metallic electrodes [10,11]; (ii) how it develops spatially [6,7]; and (iii) how it can deviate [8] from the linear dependence on the intensity of light or, equivalently, quadratic dependence on the electric field as in equation (1).
In this study, we apply time-dependent nonequilibrium Green function (TDNEGF) formalism [29,30] to compute in numerically exact and, therefore, nonperturbative fashion temporal and spatial development of photocurrent induced by light pulse of duration σ light and center frequency Ω irradiating Rice-Mele tightbinding (TB) chain with broken inversion symmetry. To mimic recent experiments [10], we use chain of length N N RM RM light > whose middle N RM light sites are irradiated by light pulse. The Rice-Mele TB chain is attached to two semi-infinite leads (modeled also as TB chains) to form a realistic two-terminal device geometry illustrated in figure 1. This makes it possible to study how fast photocurrent induced in the middle part of such device propagates toward its electrodes, and how much of photoexcited electrons can eventually be extracted into the electrodes.
The paper is organized as follows. In section 2, we introduce the Rice-Mele Hamiltonian and the DOS and DC transport properties of the devices in figure 1. This section also explains TDNEGF algorithms for constructing time-dependent nonequilibrium density matrix of this device viewed as an open quantum system. Section 3.1 discusses two-photon versus one-photon absorption mechanism which generate subgap versus above-the-gap photocurrent in this device. In section 3.2 we demonstrate shift photocurrent as a realization of nonadiabatic quantum charge pumping in two-terminal devices with broken left-right-symmetry. Section 3.3 discusses how fast photoexcited charge carriers propagate from irradiated region of the device toward the electrodes, with additional insight provided as two movies in the supplemental material, available online at stacks.iop.org/JPMATER/2/025004/mmedia 6 . The role of quantum coherence in generation of photoexcited nonequilibrium charge density is rigorously discussed in section 3.4 by analyzing contributions to it arising from off-diagonal elements of time-dependent nonequilibrium density matrix. We conclude in section 4.

Models and methods
whose parameters have been tuned to capture electronic structure of realistic materials like polyacetylene, BaTiO 3 and monochalcogenides [16,17]. Here c n † (c n ) creates (annihilates) electron in s-wave orbital nñ | centered at site n; B 2 g g = -+ + and γ − =−γ−B/2 are the alternating hoppings between the nearest neighbor sites with B parameterizing the dimerization of the chain; and ±D/2 is the staggered on-site potential. The unit cell of Rice-Mele TB chain of size 2a contains two sites, and inversion symmetry is broken by [14] B 0 ¹ and D 0 ¹ . In order to mimic devices used to extract shift current in recent experiments [10,11], we attach Rice-Mele TB chain to two semi-infinite normal metal (NM) leads depicted in figure 1, which are modeled by the same TB Hamiltonian in equation (2) but with B=D=0. We set γ=1 eV in both the NM leads and in the Rice-Mele central region where B=D=1 eV is chosen. The NM leads terminate in the left (L) and right (R) macroscopic reservoirs whose chemical potentials are identical in the absence of DC bias voltage Figure 1. Schematic view of a two-terminal device where photocurrent is generated by a femtosecond light pulse of duration σ light in the absence of any DC bias voltage between its leads. The central region composed of N RM sites is modeled by noncentrosymmetric Rice-Mele TB Hamiltonian in equation (2) with alternating hoppings and staggered on-site potential illustrated. It is attached to twosemi infinite NM leads modeled as TB chains with uniform hopping and zero on-site potential. The incident light pulse is assumed to couple only to N N RM light RM < sites in the middle of Rice-Mele central region. The insets below the device show gapless band structure of NM leads, as well as band structure of an infinite Rice-Mele TB chain with an energy gap E g between the valence and conduction bands. The Fermi energy of electrons in the NM leads is chosen as E F =0 eV, and it falls in the cente of the gap E g , as marked in all three insets. and chosen as μ L =μ R =0 eV, so that electrons in NM leads are at the Fermi energy E F =0 eV (which is band center of the NM leads).
For a sufficiently long Rice-Mele TB chain, this device exhibits an energy gap E g ≈2.28 eV in the DOS (figure 2(a)). We also consider gapless device where the gap of short Rice-Mele TB chain is filled with evanescent wave functions injected by the NM leads ( figure 2(b)). Adding on-site disorder, modeled as a uniform random variable ä[−W/2, W/2], reduces the van Hove singularities in the DOS at the gap edges while keeping finite two-terminal conductance G outside of the gap (figure 2(c)). That is, the scaling of the corresponding twoterminal resistance R=1/G in figure 2(d) with N RM shows that for chosen W=0.3 eV and N RM =100 the device is outside of the Anderson localization regime. This setup makes it possible to quantify how fast nonequilibrium photoexcited charge carriers propagate across diffusive TB chain toward the NM leads. We compute G as a linear response to small DC bias voltage between the NM leads using the zero-temperature ( ) with a Gaussian shaped function for a pulse of duration σ light , centered at time t p and center frequency Ω. Here e x is the unit vector along the direction of TB chain, pointing toward the right NM lead in figure 1, which describes one of the two possible linear polarization of incident light. The corresponding electric field is E=−∂ A/∂t. We neglect the relativistic magnetic field of the laser pulse, so that electronic spin degree of freedom maintains its degeneracy and it is excluded from our analysis. The vector potential couples to an electron via the Peierls substitution in equation (2)  is the dimensionless parameter quantifying maximum amplitude of the pulse.

TDNEGF algorithms
To compute local charges and currents driven by time-dependent terms in the Hamiltonian, we employ TDNEGF formalism which operates with two fundamental quantities [29]-the retarded G t t , Green functions (GFs) describing the density of available quantum states and how electrons occupy those states, respectively. For the device in figure 1 we solve a matrix integro-differential equation [32,33] for the time-evolution of one-particle reduced nonequilibrium density matrix, Here all bold-face quantities denote matrices of size N RM ×N RM in the vector space defined by the central Rice-Mele region of two-terminal device in figure 1. The diagonal elements of the nonequilibrium density matrix yield the local nonequilibrium charge where we subtract local charge in equilibrium. The equilibrium density matrix can be obtained either as the asymptotic limit t neq eq r r  ( ) after transient current dies away (and before the light pulse is applied) in the course of time evolution which couples the NM leads to the central Rice-Mele region, or by evaluating eq r = ( ) ( ) using the retarded GF in equilibrium and the Fermi function f (E) of macroscopic reservoirs [35]. We explicitly confirm that both methods give identical result. The matrix is expressed in terms of the lesser/greater GF and the corresponding lesser/greater self-energies [29] t t , , 2 S a > < ( ) whose numerical construction in order to convert equation (4) into a system of ordinary differential equations can be found in [33]. Equation The computational complexity of TDNEGF calculations stems from the memory effect-the entire history must be stored in order to accurately evolve the GFs. For efficient calculations over long times and for large number of simulated sites, we employ newly developed TDNEGF algorithms [32,33] which scale linearly [30] in the number of time steps.

Results and discussion
3.1. One-versus two-photon absorption mechanisms Although electrons of any system will respond to light pulse by generating a time-dependent photocurrent, photovoltaic applications and the analysis of consequences of broken symmetry are focused on the existence of nonzero DC component of the photocurrent. We define DC component as where σ curr >σ light is the duration of transient photocurrent, and plot it as a function of Ω in figure 3(a) for gapless and in figure 3(c) for gapped device. In gapped device nonzero I t 0 á ñ ¹ a ( ) appears initially at subgap frequency ÿΩ ; Δ/2, which is at first sight surprising since in the 'standard model' analyses [14,17] of an infinite Rice-Mele TB chain DC photocurrent is zero in the gap. However, it is compatible with an electron from the valence band in figure 2(a) absorbing two-photons at the same time to transition to the conduction band. Such two-photon absorption mechanism is confirmed by demonstrating in figure 3 , which is the signature of the usual one-photon absorption mechanism captured also by equation (1). We note that two-photon photovoltaic effect, as the nonlinear analog of conventional one-photon photovoltaic effect, is rarely observed in standard pn-junctions solar cells [27,28].
To include these higher-order processes into the 'standard model' analyses requires to derive an additional expression [37] for the fifth-rank tensor σ pqrst (ω) in order to express DC component of subgap photocurrent, On the other hand, they are naturally taken into account by our nonperturbative TDNEGF approach. Also, the TDNEGF approach applied to two-terminal devices allows one to compute the open-circuit photovoltage simply by using where the two-terminal resistance R (such as the one plotted in figure 2(d)) depends on the disorder and possible domain walls [38] in the central region, as well as properties of the central-region/NM-lead interfaces.  potentials oscillating out-of-phase [40], which leads to I t á ñ µ W a ( ) at low frequencies. This is confirmed in figure 4(c) for Rice-Mele TB chain driven by two on-site potentials oscillating out-of-phase (illustrated as Device 1 in figure 4(e)). In contrast, in the nonadiabatic regime [39][40][41], only one of those two symmetries needs to be broken and this does not have to occur dynamically. The DC component of the pumped current in the nonadiabatic regime is [41,42]  We note that previous theoretical analyses [14] of shift current in bulk materials have concluded that those with broken inversion symmetry but without electric polarization generate DC component of photocurrent only in response to polarized light. Thus, the nonzero polarization vector is required to generate DC component of photocurrent in response to unpolarized light, albeit larger polarization does not always imply a larger photocurrent [14,17].

Shift current as nonadiabatic quantum charge pumping
However, these arguments do not consider realistic devices in two-terminal geometry for which the theory of nonadiabatic quantum charge pumping [39][40][41] predicts how leads made of different materials, or identical leads and static on-site potentials within the central region [42], break the left-right symmetry of the device structure to generate nonzero DC component of photocurrent in response to unpolarized light. We confirm the former possibility in figure 4(d) by using TB chain of finite length with uniform hoppings and zero on-site potential as the central region of Device 3 in figure 4(e). This is irradiated by light pulses of different polarizations (+ sign denotes polarization along the TB chain toward the right NM lead and − sign denotes polarization in the opposite direction), while the left-right symmetry of the device is broken by using two different NM leads. This

Superballistic spreading of photoexcited charge carriers
Motivated by very recent experiments [10,11], exploring position dependence of photocurrent induced by applying CW light across the device or its temporal waveform induced by femtosecond light pulse, we examine spatial profiles of nonequilibrium charge ( figure 5(c)), Q t i neq ( ) in equation (5), and local bond current  To understand how fast the spreading of nonequilibrium charge is toward the NM leads, we analyze spatiotemporal profiles depicted in figure 5(a) which trace those sites of clean or disordered Rice-Mele TB chain where absolute value Q t i neq | ( )| at time t reaches 5% (other cutoffs can be used without changing the conclusion) of the maximum value generated (at t=717.8 fs in panel (c)) within the irradiated region. Thus, upper and lower curve in figure 5(a) can be viewed as the displacement of Q i neq toward the right or left NM lead, respectively. The scaling of the displacement with time, ∼t ν , can be analyzed akin to variance spreading of optical [43,44] or quantum [45] wave packets or classical Brownian particle [46]: ν=1 signifies ballistic propagation in uniform lattices; ν=0.5 or ν=0 signifies diffusion or Anderson localization in disordered lattices, respectively; and subdiffusion (0<ν<0.5) or superdiffusion (0.5<ν<1) are also possible in some quasiperiodic lattices [47]. Surprisingly, we find ν=1.42 (11) in the case of clean Rice-Mele TB chain and ν=1.28 (15) in the case of the diffusive one, which demonstrates superballistic spreading of photoexcited nonequilibrium electrons.
We note that superballistic spreading of optical wave packets, within certain transient time frame, has been observed experimentally in disordered static [43] and temporally fluctuating [44] photonic lattices. The farfrom-equilibrium quantum electron system studied in figure 5 does share some features with the latter case. The superballistic transport of photoexcited nonequilibrium charge carriers unveiled in figure 5, which evolve quantum-coherently and are largely insensitive to scattering off impurities, is remarkably different from photoexcited carriers in conventional pn-junction solar cells where they travel toward the leads via driftdiffusive transport requiring to manipulate their mobilities for efficient DC photocurrent extraction. The insensitivity to impurities of superballistic propagation of quantum-coherent shift current in figure 5 can explain the central result of experiments in [10] observing robust transport of photoexcited carriers without decay over surprisingly long distances ∼100 μm.

Diagonal versus off-diagonal elements of time-dependent nonequilibrium density matrix
Experimentally, the ballistic and shift contributions to BPVE can be distinguished by performing Hall effect measurements because ballistic current is sensitive to external magnetic field and shift current is not [4,9]. Theoretically, the quantum-mechanical nature of shift current is encoded by the off-diagonal elements [24,25] of t neq r ( ), in contrast to ballistic current arising from its diagonal elements [24,25]. However, standard theoretical analysis is performed for infinite periodic crystals, which has lead to confusion [48] about how to identify shift current of carriers which are not illuminated by light [10]. The preferred basis [49] for analyzing quantum coherence encoded in the off-diagonal elements of the density matrix is typically the real-space basis, so akin to 'total site coherence' of [50] we define the following sum of amplitudes of the off-diagonal elements of t signifies processes which involve quantum coherence or superposed states. Figure 6 shows that relative change of the diagonal elements (figures 6(a), (c)) of the density matrix upon transition from equilibrium to nonequilibrium is 1%, in both irradiated and non-irradiated regions of Rice-Mele TB chain. On the other hand, relative change of the off-diagonal elements (figures 6(b), (d)) reaches 100%.
We also perform analysis of matrix elements of t  ¹ belong to the same (valence or conduction) band. The off-diagonal intraband contribution in figure 6(h) is almost exactly canceled by the diagonal contribution (for which E q =E r ) in figure 6(g).Thus, this reveals that photoexcited nonequilibrium charge density far away from irradiated region is governed solely by quantum coherence of superposed electron and whole wavefunctions.

Conclusions
Traditional theoretical analyses of photocurrent in noncentrosymmetric systems have been confined to perturbative treatment of an infinite periodic crystal homogeneously irradiated by continuous wave light. These assumptions are at odds with recent experiments where femtosecond light pulses are applied locally and propagation of photocurrent is studied away from the irradiated region. Furthermore, possibly high intensities of light make it desirable to include higher-order processes that would be naturally accounted for via a nonperturbative approach. Here we make a leap in the theoretical description by moving from perturbative to nonperturbative, and by moving the focus from the bulk material to a realistic device geometry perspective. Specifically, the nonperturbative TDNEGF approach employed in this study makes it possible to analyze spatiotemporal dynamics of photocurrent in realistic devices attached to external leads and exposed to either CW or pulse light of arbitrary intensity. We predict: subgap photocurrent generated by two-photon absorption mechanism; above the gap photocurrent generated by one-photon absorption mechanism; and its superballistic propagation away from irradiated region in both clean and diffusive (for conventional DC charge transport) systems. Importantly, we demonstrate that beyond the traditional target materials (i.e. noncentrosymmetric materials with nonzero polarization vector), a much broader class of systems can be explored to optimize the photocurrent in response to unpolarized light (as relevant for photovoltaic applications). The key requirement defining this broader class of systems is the same as for nonadiabatic quantum charge pumping-broken leftright symmetry of the device structure. For example, it is sufficient to take a two-dimensional material (even non-polar) which interacts with light strongly [51,52] and attach it to two electrodes made of different metals in order to enable quantum-coherent shift current generation with nonzero DC component even when irradiating light is unpolarized.