Local photo-excitation of shift current in noncentrosymmetric systems

TL;DR

This paper theoretically demonstrates that local photo-excitation in noncentrosymmetric materials can induce a shift current independent of excitation position or system size, highlighting its potential for efficient photovoltaic devices.

Contribution

It introduces a theoretical model showing local photo-excitation can generate shift current regardless of excitation location or system length, unlike conventional photocurrent.

Findings

Shift current is independent of excitation position and system size.

Local photo-excitation can induce a persistent shift current.

Contrast with conventional photocurrent suppression at local excitation.

Abstract

Photocurrent in solids is an important phenomenon with many applications including the solar cells. In conventional photoconductors, the electrons and holes created by light irradiation are separated by the external electric field, resulting in a current flowing into electrodes. Shift current in noncentrosymmetric systems is distinct from this conventional photocurrent in the sense that no external electric field is needed, and, more remarkably, is driven by the Berry phase inherent to the Bloch wavefunction. It is analogous to the polarization current in the ground state but is a d.c. current continuously supported by the nonequilibrium steady state under the pumping by light. Here we show theoretically, by employing Keldysh-Floquet formalism applied to a simple one-dimensional model, that the local photo excitation can induce the shift current which is independent of the position and width of the excited region and also the length of the system. This feature is in stark contrast to the conventional photocurrent, which is suppressed when the sample is excited locally at the middle and increases towards the electrodes. This finding reveals the unconventional nature of shift current and will pave a way to design a highly efficient photovoltaic effect in solids.

Equations

$$H=H_0+H_A,$$

$$H_0$$

$$H_A$$

$$J_s(x_A)\equiv \frac{J(x_A)+J(1-x_A)}{2}$$

$$J_{as}(x_A)\equiv \frac{J(x_A)-J(1-x_A)}{2},$$

$$\delta n_c=\frac{W_{vc}{\rho }_v}{\Gamma +W_{vc}({\rho }_c+{\rho }_v)}.$$

$${\rho }_c\delta n_c={\rho }_v\delta n_v.$$

$$\delta n_v=\frac{W_{vc}{\rho }_c}{\Gamma +W_{vc}({\rho }_c+{\rho }_v)}.$$

$$\delta n_c=\frac{1}{2\pi N}\Im \left\{\mathrm{Tr}{\int }_0^{\infty }d\omega G_{00}^{<}(\omega )\right\},$$

$${\rho }_c=\frac{1}{2\pi N}\Im \left\{\mathrm{Tr}{\int }_0^{\infty }d\omega \left[G_{00}^A(\omega )-G_{00}^R(\omega )\right]\right\}.$$

$$G_{\kappa {\kappa }^{\prime }}(\omega )$$

$$G_n(\omega )=\frac{1}{T}{\int }_0^{T}d{t}_a\int d{t}_{r}G(t_a+\frac{t_r}{2},t_a-\frac{t_r}{2})e^{\mathrm{i}(n\Omega t_a+\omega t_r)}$$

$$G(t,t^{\prime })=G(t+T,t^{\prime }+T),$$

$$(\omega +\kappa \Omega )$$

$${\mathcal{H}}_{\kappa ,{\kappa }^{\prime }}^0=\frac{1}{T}{\int }_0^{T}dt{H}^0(t)e^{\mathrm{i}(\kappa -{\kappa }^{\prime })\Omega t}$$

$$H(t)$$

$${\left[G_{\kappa ,{\kappa }^{\prime }}\right]}_{k\alpha ,i}(\omega )=\sum _j{\left[g_{\kappa ,\kappa }\right]}_{k\alpha ,k\alpha }(\omega )V_{k\alpha ,j}{\left[G_{\kappa ,{\kappa }^{\prime }}\right]}_{j,i}(\omega ),$$

$$(\omega +$$

$${\left[{\Sigma }_{\kappa ,\kappa }^{\prime }\right]}_{i,j}(\omega )=\sum _{k\alpha }{V}_{i,k\alpha }{\left[g_{\kappa ,\kappa }\right]}_{k\alpha ,k\alpha }(\omega )V_{k\alpha ,j}.$$

$${\left[{\Sigma }_{\kappa ,\kappa }^{\prime }(\omega )\right]}_{i,j}$$

$${\Gamma }_{\alpha }=2\pi \sum _k|V_{k\alpha ,i_{\alpha }}{|}^2\delta (\omega -{ϵ}_{k\alpha }).$$

$${\Gamma }_{\alpha }(\omega )={\Gamma }_{\alpha }=\Gamma .$$

$$J_{\alpha }$$

$$J_{\alpha }=i{q}_e\int \frac{d\omega }{2\pi }\mathrm{tr}$$

$${\left[{\hat{\Gamma }}_{\alpha }\right]}_{i,j}(\omega )=2\pi \sum _{k\alpha }{V}_{j,k\alpha }{V}_{k\alpha ,i}\delta (\omega -{ϵ}_{k\alpha }),$$

$${\left[{\hat{\Gamma }}_{\alpha }\right]}_{i,j}(\omega )={\delta }_{i,i_{\alpha }}{\delta }_{j,i_{\alpha }}{\Gamma }_{\alpha }(\omega )$$

$$J_{\alpha }=i{q}_e\int \frac{d\omega }{2\pi }$$

$$J$$

$$J_a$$

$$=$$

$$(\omega +\kappa \Omega )\hat{1}-{\mathcal{H}}^0-{\Sigma }^{\prime },$$

Full text

Keywords: photovolatic effect, shift current, nonequilibrium Green’s function method.

1 Introduction

Photo-carrier generation and the associated photocurrent is a subject of intensive studies in condensed matter physics and its applications. For example, the excitation spectrum of photocurrent together with the optical absorption spectrum offers a powerful tool to investigate the electronic states in bulk solids. Namely, the absorption can occur both by the excitons, i.e., the bound states of electron and hole, and electron-hole continuum, while the photocurrent can be carried only by the latter, as the exciton is a charge neutral object and cannot contribute to the current. By increasing the external electric field for the photocurrent, the excitons split into electrons and holes and begins to contribute to the photocurrent. In solar cells, the external electric field is generated by the p-n junction, where the electrons and holes feel the potential gradient in the opposite directions and are separated [1, 2].

In contrast to the conventional view of the photocurrent discussed above, the shift current in noncentrosymmetric systems has been attracting intensive interest [3, 4, 5, 6, 7, 8, 9, 10, 11]. In this phenomenon, since the system itself has a fixed “direction”, the light irradiation alone can produce the d.c. current without the external electric field. One might think that the internal electric field plays a similar role to that of the external one. However, the situation is much more subtle. Due to the periodicity of the crystal potential, there is no macroscopic potential gradient. Microscopically, the band structure has certain symmetry between $\vec{k}$ and $-\vec{k}$ due to the time-reversal symmetry $T$ ( $\vec{k}$ is the crystal momentum ). Namely, $T$-symmetry gives the relation $\varepsilon_{\uparrow}(\vec{k})=\varepsilon_{\downarrow}(-\vec{k})$ with $\uparrow(\downarrow)$ represents the up (down) spin. Therefore, it is highly nontrivial whether the excitation of electrons and holes contributes to the shift current. The lack of inversion symmetry $I$, however, opens the possibility of the Berry phase of Bloch wavefunction to be finite. (More precisely, one can choose the gauge of the Bloch wavefunction where the Berry phase connection is zero in the presence of both $T$ and $I$ symmetries, i.e., no Berry phase contribution to material properties.) It have been shown that, the shift current comes from the Berry connection of the Bloch wavefunction $|\vec{k}\,n\rangle$ for conduction ($n=c$) and valence $(n=v$) bands [12, 13, 11]. An important observation here is that the Berry connection has the meaning of the intracell coordinate in the sense that the position $\vec{r}_{n}$ of the wavepacket made from the Bloch wavefunctions shifts proportionally by the Berry connection. Indeed, it is related to the polarization in ferroelectric materials [14, 15].

This scenario for the ground state finds its generalization to the nonequilibrium steady state, where the continuous pump from the valence to conduction band occurs by the photo-excitation. In this case, the difference of the Berry phase between the conduction and valence bands corresponds to the shift of the electron position induced by the transition, and contribute to the d.c. photocurrent. This is the picture of shift current, and hence is different from the transport current supported by the free carriers. It is more quantum mechanical origin, i.e., cannot be attributed to the motion of individual carriers, and is expected to be nonlocal in nature. Therefore, the optical response of noncentrosymmetric photovoltaic systems to local photo-excitation is nontrivial.

To study the shift current due to the local photo-excitation, in this paper, we consider a noncentrosymmetric non-interacting fermion chain as a basic model that exhibits shift current. As we ignore the Coulomb interaction of electrons, our model has no internal bias due to the ferroelectric polarization, therefore, no photocurrent from the conventional mechanism similar to photoconductors. We study the model by a nonequilibrium Green’s function technique used to study nonequilibrium phenomena in semiconductors [16] and in organic solids [17, 18, 19]; the details of the method are elaborated in appendices. By applying the method, we numerically study the shift current induced by the local photo-excitation using the fermion chain up to $N=1200$ sites. Although there are no potential gradient, we find that a finite photocurrent appears in this model. Furthermore, we find that the locally excited shift current in the fermion chain behaves in a distinct manner from the photocurrent in photoconductors and in p-n junctions, as shown in Fig. 2. In particular, we find that the magnitude of shift current does not depend on the position of the photoexcitation. This is a contrasting feature to the conventional case, in which the current is suppressed for the photoexcitation at the center of the chain. This is an interesting feature that experimentally distinguishes the shift current from the conventional photocurrent. Furthermore, this is an advantageous feature for application to highly efficient optoelectronic devices.

We also find that the magnitude of photocurrent does not change by increasing the width of the bonds excited. As presented in Sec. 3.2, we discuss that this is likely to be a consequence of high carrier density excited by light.

The remainder of the paper is organized as follows. In Sec. 2, we introduce the model we consider in this paper. In Sec. 3, we present the result of theoretical calculations on the light-position dependence of the shift current in a fermion chain. Sections 4 and 5 are devoted to discussions and summary, respectively. The nonequilibrium Green’s function method for the Floquet theory and its numerical implementation are elaborated in A and B, respectively. The detail of the quasiclassical theory for carrier concentration and the shift current is explained in C.

2 Model and Method

To study the effect of local photo-excitation in noncentrosymmetric systems, we consider Rice-Mele model [20] coupled to leads (electrode) at each end and with monochromatic incident light [Fig. 1(a)]; the light is taken into account as the time-dependent electric field. The Hamiltonian for the chain is given by

[TABLE]

where the first term is the Rice-Mele Hamiltonian,

[TABLE]

and the second term is the interaction with the incident light,

[TABLE]

Here, $c_{i}$ ($c_{i}^{\dagger}$) is the annihilation (creation) operator of the electron at the $i$th site; $t$ ($B$) is the uniform (staggered) hopping between the nearest-neighbor sites, and $d$ is the staggered potential. In Eq. (3), $\Omega$ is the frequency of the light and $A$ is the amplitude of vector field. Here, we only consider the leading order in $A$. The sum $\sum^{\prime}$ is over consecutive $l_{A}$ sites centered at $i_{A}$ [see Fig. 1(a)].

In periodically driven systems, there are no eigenstates as in the time-independent Hamiltonians. However, according to the Floquet theory [21, 22], one can define similar eigenstates for the periodically driven Hamiltonian called Floquet states. In the Floquet theory, in addition to the bands that consists of the eigenstate of time-independent Hamiltonian, we consider the side bands of the original band with the energy shift of $-\kappa\Omega$ ($\kappa\in\mathbb{Z}$) [Fig. 1(b)]. When the chain is exposed to light, it gives the mixing term between the side bands with different $\kappa$; the leading order contribution from the light comes in as the mixing between two neighboring subbands (subbands with the indices $\kappa$ and $\kappa+1$). Therefore, to take into account of the leading order contributions from the light, in the subsequent calculations, we focus on the bands with indices $\kappa=0$ and $-1$ bands and ignore all others (B).

To study the nonequilibrium transport phenomena in the Rice-Mele model, the photocurrent from the local excitation is studied using a nonequilibrium Green’s function theory [23, 24, 16]. We generalized the method to the case of periodically driven systems and computed the shift current numerically. The general formalism and further details on the application of the method to the model in Eq. (1) are elaborated in A.

3 Results: Shift Current in Noncentrosymmetric Crystals

3.1 Light-Position Dependence of Shift Current

In this section, we study the shift current induced by the local excitation by light. Figure 2 shows the light-position dependence of the shift current by the local excitation at $x_{A}=i_{A}/N$ for $N=1200$. The three different lines correspond to different $l_{A}$, the number of bonds we shine the light onto [see fig. 1(a)]. The results show that the shift current by the local excitation is almost independent of $x_{A}$. Besides, the shift current is symmetric with respect to $x_{A}$, i.e., $J(x_{A})\simeq J(1-x_{A})$. For later convenience, we here define the “symmetric” and “asymmetric” components of the photocurrent by

[TABLE]

and

[TABLE]

respectively. By calculating $J_{s}$ and $J_{as}$ from the data in Fig. 2, we indeed find that $J_{as}\lesssim 10^{-8}$, which is below the numerical precision. Hence, our results show that the symmetric part of the photocurrent is dominant in the Rice-Mele chain.

This is a contrasting feature from what is expected in the conventional mechanism in p-n junctions, by which the photocurrent is borne by charged quasi-particles. In p-n junctions, the photocurrent is essentially an interfacial phenomenon. At the interface (or surface), a potential gradient occurs due to the nonuniform structure. Due to this potential gradient, when an electron-hole pair is created by light, the electron and hole drift toward opposite directions, resulting in a net charge current. In the setup shown in Fig. 1(a), the junction of the chain and the lead plays the role of the interface. Therefore, in this mechanism, $J$ is expected to increase as we approach the edges. Besides, if the carriers are generated close to the edge, we expect that the excited carriers can escape to the electrodes easily, before being scattered back into the valence band. Hence, the current is expected to be larger at the ends of the chain than in the center.

Another important feature is the symmetry. Although there is no potential gradient in the bulk, in principle, a finite photocurrent can appear from the conventional mechanism when the position of the light is sufficiently close to one of the edges. This is a mechanism similar to p-n junctions where the potential gradient appears due to the interface. If this mechanism is dominant, for a centrosymmetric chain (or if the effect of broken inversion symmetry is negligible), the photocurrent from the local excitation is expected to be antisymmetric. On the contrary, if the photocurrent is induced by the bulk effect, e.g., anomalous photovoltaic effect, then the $x_{A}$ dependence of the photocurrent is expected to be independent of the position of light for a sufficiently long chain, with possibly some effect of the surface near the edges, i.e., they are (almost) symmetric. Therefore, we expect the light-position dependence of the locally excited photocurrent reflects the microscopic mechanism of the photocurrent; $J_{as}$ represents the contribution of the edges, e.g. the photocurrent by the conventional mechanism at the edge, while $J_{s}$ reflects the bulk effect such as the anomalous photovoltaic effects or shift current.

These observations imply that the symmetric feature of the photocurrent with respect to $x_{A}$ reflects the unconventional photocurrent in the noncentrosymmetric chains, namely, the shift current. Therefore, in this paper, we call the symmetric photocurrent “shift current”. Our result in Fig. 2 indicates that shift current is the dominant source of the photocurrent in Rice-Mele chain. To compare the result of Rice-Mele model with a model with no shift current, we also performed the calculation of $x_{A}$ dependence for a symmetric chain with $t=1$, $B=1.2$, and $d=0$ (not shown). In this model, however, we find no observable photocurrent within the numerical precision (only small $J_{as}$ that is below our numerical precision), consistent with the above discussion.

In the last, we discuss the effect of finite size effect. Fig. 3 shows the result of the asymmetric chain for $l_{A}=39$ with $N=600$, $900$, and $1200$. The results for the different $N$ show similar magnitude of the shift current. They are also symmetric with respect to $x_{A}$. From this observation, we conclude that our data for $N=1200$ reflects $N\to\infty$ behavior in the ballistic limit.

3.2 Light Width Dependence of Shift Current

Another point to be noted in Fig. 2 is the $l_{A}$ dependence of the shift current. In Fig. 2, the induced shift current appears to be roughly the same for all $l_{A}$. Naively, the photocurrent is expected to be proportional to the amount of the photon energy absorbed by the chain, therefore, proportional to $l_{A}$. Our results, however, show that there is nearly no dependence of the induced current on $l_{A}$. In our calculation, such behavior was observed for all range of $\Omega$. Figure 4 shows the result of shift current excited by the local excitation at $x_{A}=0.5$. Figure 4(a) is the result for different $l_{A}$ in $N=1200$ noncentrosymmetric chain. A finite photocurrent is observed above $\Omega\sim 2.6$ and form a triangular shape with a peak at $\Omega\sim 3.2$. The lower limit $\Omega\sim 2.6$ roughly corresponds to the band gap $\sqrt{\Delta^{2}+4B^{2}}=\sqrt{5}\sim 2.24$, while the peak corresponds to an $\Omega$ which a subband of the conduction band perfectly overlap with that of the valence band, $\frac{1}{2}(\sqrt{\Delta^{2}+16t^{2}}+\sqrt{\Delta^{2}+4B^{2}})=\frac{1}{2}(\sqrt{5}+\sqrt{17})\sim 3.17$. The results show almost no $l_{A}$ dependence for all range of $\Omega\in[2.5:4.0]$, where the shift current is observed.

Figure 4(b) shows the $N$ dependence of the shift current with $l_{A}=39$ and $A=0.2$ at $x_{A}=0.5$. For most values of $\Omega$, the results of shift current show almost no dependencies to $N$. On the other hand, some finite size effect is seen at $\Omega\sim 3.1$. The $N$ sensitive behavior around $\Omega\sim 3.1$ is likely to come from the sparse density of states in the center of the band. In the Floquet theory, the effect of light primarily appears around the crossing of two bands with different Floquet indices $\kappa$. In our model with $t=B=d=1$, it is around the crossing of the conduction band with $\kappa=0$ and valence band with $\kappa=-1$ (see Fig. 1); the two bands touch at their edge for $\Omega\sim 2.24$, and crosses at the center of the bands for $\Omega\sim 3.17$. As the density of states is relatively sparse at the center of the band than at the edges, the discreteness of the energy levels due to the finite length become more evident at the band center. However, the $N$ dependence at $\Omega\sim 3.1$ also appears to be converged for the chains we used here ($N=1200$).

The result indicates that, for a sufficiently long chain, there seems to be almost no $l_{A}$ dependence. Furthermore, from the $\Omega$ dependence, it seems the $l_{A}$ independent behavior has less to do with the microscopic structure of the wave functions, e.g., Berry connection. According to the Floquet theory, as the effects of light have much smaller energy scale than other parameters, the effect of light on the electrons appears primarily around the band crossing point, where $\varepsilon_{v}(k)+\Omega=\varepsilon_{c}(k)$ [28]. Here, $\varepsilon_{v}(k)$ is the energy eigenstate of the valence band with wavenumber $k$, and the $\varepsilon_{c}(k)$ being that of conduction band [see Fig. 1(b)]. As the crossing point changes by changing $\Omega$, if the behavior reflects the microscopic structure of the electronic states, e.g., Berry phase, then it is natural to expect $\Omega$ dependence. Therefore, it is unlikely that the $l_{A}$-independent behavior is related to the mechanism of the shift current.

From this observation, the $l_{A}$-independent behavior is likely to come from the fact that a large number of carriers are excited in the conduction band. To show this, we here consider a quasi-classical model schematically shown in Fig. 5. The model consists of valence band, conduction band, and the lead(s) as shown in Fig. 5. The carriers are transferred between the valence and conduction bands by transition rate $W_{vc}$ and between the conduction band and the lead by hybridization $\Gamma$. Further details on the model and the calculation are given in C.

From the general formalism of the Green’s function theory, the charge current induced in the model is proportional to the product of $\Gamma$ and the change of electron distribution function from the equilibrium distribution, as shown in Eq. 34. Therefore, we here focus on the carrier density in the conduction (valence) band, $\delta n_{c}$ ($\delta n_{v}$), assuming that the current is proportional to $\delta n_{c}$ and $\delta n_{v}$. In our case, $W_{vc}$ comes from the irradiation of light. A perturbation theory based on Fermi’s golden rule shows that this term is proportional to $l_{A}$.

In the steady state ($\delta\dot{n}_{c,v}=0$), the calculation based on this simple model gives

[TABLE]

Here, we used charge neutrality condition,

[TABLE]

Similarly, we get

[TABLE]

Therefore, for $W_{vc}\ll\Gamma$, the carrier densities reads $\delta n_{c,v}\propto W_{vc}\propto l_{A}$; the number of carriers and the photocurrent is proportional to $l_{A}$. On the other hand, if $W_{vc}\gg\Gamma$, $\delta n_{c}$ and $\delta n_{v}$ become (almost) independent of $W_{vc}$, i.e., it is independent of the width of excitation.

From the numerical calculations, it is difficult to directly give an estimate of $W_{vc}$. However, we can investigate to which limit the system is close to by calculating the density of excited carriers. If the system is in $W_{vc}\gg\Gamma$ limit, there exists a large density of carriers. On the other hand, the carrier density is about $\delta n_{c}\sim W_{vc}\rho_{v}/\Gamma$ for $W_{vc}\ll\Gamma$; the carrier density is much smaller than the density of state. In Fig. 6, we calculated the carrier density in the conduction band with respect to the number of states, $\delta n_{c}/\rho_{c}$. The number of carriers is calculated by

[TABLE]

and the number of the number of states by

[TABLE]

Here, $G^{<}_{00}(\omega)$ is the $N\times N$ matrix of lesser Green’s function and $G^{R,A}_{00}(\omega)$ is the retarded and the advanced Green’s function with Floquet index $\kappa=0$. The $l_{A}$ dependence of $\delta n_{c}$ is plotted in Fig. 6(a). Here, the calculation were done for $l_{A}=39$ and $x_{A}=0.5$ with $A=0.2$. The three curves were for the different length of chain, $N=600,900,1200$. For the noncentrosymmetric chain, the result shows a rapid increase of the carrier density below $l_{A}=20$ and a saturation to $\delta n_{c}/\rho_{c}\sim 0.25$ above that. In coincidence with the saturation of the carrier density, the photocurrent also shows saturation above $l_{A}=20$. These features are consistent with the above mechanism for the $l_{A}$-independent behavior of the shift current.

3.3 Current-Voltage Characteristic

In this section, we study the current-voltage characteristic of the chains with shift currents by applying the external bias to the leads on the both sides. In our calculation, however, we find no current dependence on the external voltage within the numerical precision, if the bias is smaller than the charge gap (not shown). This is a natural consequence as the chain have no density of states at the Fermi level. As the carriers have higher energy than the bias, the effuluence of carriers to the leads is not affected by the change of electron distribution in the leads, which only takes place close to the Fermi level.

4 Discussion

Our results presented in Fig. 2 show that a large shift current can be induced in a ferroelectric chain by local excitation. In the clean limit of a fermion chain, our result implies that the photocurrent in bulk photovoltaic materials can be robustly transported through the material for a long distance without decaying. This sheds light on to the potential realization of a highly efficient photovoltaic device using ferroelectric insulators [4, 5, 6, 7, 8, 9, 10].

An important point related to the robustness of the current is the effect of impurities. In experiments, the impurities often suppress electric current. In one-dimensional systems, the fermions are susceptible to Anderson localization. If the localization occurs, the current is likely to be suppressed at the scale of the localization length, i.e., no d.c. current is transported for the chains much longer than the localization length. Therefore, we expect the shift current to decay by the length scale of localization when a light is shined onto only a part of the chain.

In addition to the localization, another potential cause of the dissipation is the inelastic scattering by phonons. The inelastic scattering may contribute to the decay of shift current, limiting the distance the current is transmitted. Study on the effect of electron-phonon scattering to the shift current is left for future works.

5 Conclusion

In this paper, we theoretically studied the shift currents from local excitation by light. Considering a noncentrosymmetric insulating fermion chain, we studied the dependence of shift current on the position of excitation. By applying a non-equilibrium Green’s function theory in 2, we show that the local excitation in the noncentrosymmetric chain induces photocurrent that is almost independent of the position of excitation and is symmetric to the location of light ($x_{A}$), as shown in Fig. 2. This is in contrast to the conventional mechanism, where the photocurrent is antisymmetric and maximized at the edge of the device, where the carriers can easily escape to the coupled leads. This difference in the position dependence of the photocurrent provides an experimental method to distinguish the unconventional photocurrent from that of the conventional ones, such as in the semiconductor junctions. Furthermore, this feature is advantageous for highly efficient optoelectronic devices.

Acknowledgements

The authors thank M. Kawasaki, M. Nakamura, M. Sotome, Y. Tokura, and M. Ueda for fruitful discussions. This work was supported by JSPS Grant-in-Aid for Scientific Research (No. 16H06717, No. 24224009, No. 26103006, and No. 26287088), from MEXT, Japan, and ImPACT Program of Council for Science, Technology and Innovation (Cabinet office, Government of Japan).

Appendix A Nonequilibrium Green’s Methods

To study the transport phenomena of the model in Eq. 1 coupled to the leads, we use a nonequilibrium Green’s function method for time-independent Hamiltonians [23, 16]. In this section, we summarize the basic formulation we use to study the photocurrent. The details of the numerical implementation are given in B.

A.1 Keldysh-Floquet Theory for Systems Coupled to Leads

For our calculation, we use an extension of non-equilibrium Green’s function approach used widely to study non-equilibrium phenomena [23, 16]; the theory is extended to Floquet space in order to study periodically driven systems [22, 16]. In this section, we elaborate the general formalism and the derivation of Dyson’s equation which we use to calculate the Green’s function. We also elaborate the application to the model we study in the main text.

In the Floquet theory, the Green’s functions have additional indices that corresponds to the Floquet states [21, 27],

[TABLE]

where

[TABLE]

is the Wigner representation of the Green’s function,

[TABLE]

with $T=2\pi/\Omega$ being the periodicity of the Hamiltonian in the time direction. Using this notation, the general form of Dyson’s equation for the free particle systems is generally given by

[TABLE]

Here,

[TABLE]

is the Fourier transform of the single particle Hamiltonian.

As the Hamiltonian, we here consider a free fermion system coupled to free fermion leads,

[TABLE]

Here, $c_{i}$ ($c_{i}^{\dagger}$) is the annihilation (creation) operator for the $i$th state in the system and $d_{k\alpha}$ ($d_{k\alpha}^{\dagger}$) is the annihilation (creation) operator for $k$th eigenstate in the $\alpha$th lead. $\epsilon_{k\alpha}$ is the eigenenergy for the $k$ th state in the $\alpha$th lead, and $V_{k\alpha,i}$ is the hybridization between the $k\alpha$th state and $i$th state in the system.

For this type of models, the leads can be treated as a self-energy. The explicit form of self-energy can be deduced as follows. From Eq. 14, we get the relation

[TABLE]

where $g_{\kappa,\kappa^{\prime}}(\omega)$ is the one-particle Green’s function for $V_{k\alpha,i}=0$. By substituting Eq. 17 to Eq. 14, the Dyson equation for $\left[G_{\kappa,\kappa^{\prime}}\right]_{i,j}(\omega)$ reads

[TABLE]

with

[TABLE]

Not that, in Eq. (LABEL:eq:dyson_ren), the Hilbert space is reduced to the degrees of freedoms in the chain, not the direct sum of chain and the leads.

In this paper, we particularly consider the case in which the leads are coupled only to the ends of the chain. From Eq. (19), the self-energy that comes from the leads, $\Sigma^{\prime}_{\kappa,\kappa}(\omega)$, reads

[TABLE]

where $\alpha=L$ and $\alpha=R$ corresponds to left and right leads, respectively, and

[TABLE]

Here, $f_{R}(\omega)$ [$f_{L}(\omega)$] is the distribution function of electrons in the right (left) lead, and $i_{\alpha}$ is the site the lead is coupled to: $n_{L}=0$ and $n_{R}=N-1$, where $N$ is the number of sites in the fermion chain. In our calculation, we employed the wide-band approximation, namely, assuming $\Gamma_{\alpha}(\omega)$ to be a constant

[TABLE]

A.2 Calculation of the Current

By using a similar procedure we used to derive our Dyson’s equation, the charge current that flows into/out of the system can also be calculated from the Green’s function of the chain [23]. In our setup, the current between the $\alpha$th lead and the system is given by

[TABLE]

where $V_{i,k\alpha}=(V_{k\alpha,i})^{\ast}$. Here, $\kappa\in\mathbb{Z}$ is an arbitrary number. Following a similar procedure in Ref. [23] and by substituting $\left[G^{<}_{\kappa,\kappa}\right]_{i,k\alpha}(\omega)$ by Eq. 17, Eq. (27) reads

[TABLE]

Here, $\hat{\Gamma}_{\alpha}(\omega)$ is a square matrix with indices that corresponds to each state of the system,

[TABLE]

and $\hat{G}^{<}_{\kappa\kappa^{\prime}}(\omega)$ [$\hat{G}^{r}_{\kappa\kappa^{\prime}}(\omega)$, $\hat{G}^{a}_{\kappa\kappa^{\prime}}(\omega)$] is the matrix for lesser (retarded, advanced) Green’s functions with $\kappa$ and $\kappa^{\prime}$ indices for the Floquet states.

In the case of the model in Eq. (1), $\hat{\Gamma}_{\alpha}(\omega)$ reads

[TABLE]

for $\alpha=L,R$. Hence, Eq. (LABEL:eq:current) become

[TABLE]

Using these results, the current that flows through the system is given by [23]

[TABLE]

In the wide band approximation, Eq. (31) reads

[TABLE]

In the second equation, we set $\Gamma_{\alpha}=\Gamma$.

Appendix B Numerical Calculations

For numerical calculation using the method explained in A, we introduce some approximation to implement the theory onto the computer. In addition, to accelerate the calculation, we implemented a matrix decomposition method. In this section, we illustrate the details on the numerical implementation of the nonequilibrium Green’s function formalism.

B.1 Truncation of the Dyson Equation

The Green’s function for our model in Eq. (1) can be calculated from Eq. (LABEL:eq:dyson_ren), by inverting the matrix

[TABLE]

where ${\cal H}^{0}$ is a matrix of ${\cal H}^{0}_{\kappa\kappa^{\prime}}$ given by

[TABLE]

and $\Sigma^{\prime}$ is defined in the same manner with $\Sigma^{\prime}_{\kappa\kappa^{\prime}}$; $\hat{1}$ is the unit matrix. This matrix, however, is an infinite dimension matrix due to the Floquet indices. In our calculation, as we are interested in the leading order effect by the external light, we truncate the Floquet indices and leave the states that are close to the Fermi level, i.e., $\kappa=0$ and $-1$. Physically, this corresponds to considering only the leading order of scattering between the electrons in the valence band and that in the conduction band induced by light.

One point to be noted is that, due to this approximation, the periodicity of the Green’s function with respect to $\omega$ is violated. To reduce the error that arises from the above approximation, for the calculation of shift current using Eq. (34), we used $G_{00}(\omega)$ for $\omega\geq 0$ and $G_{-1-1}(\omega)$ for $\omega<0$.

B.2 Eigendecomposition of $G^{-1}(\omega)$

For the calculation of the Green’s functions in long chains, an efficient algorithm based on the eigendecomposition of matrices is used. Here, we briefly review the algorithm we used in this calculation. Suppose we have a $\omega$-independent Hamiltonian $H$ and self-energy $\Sigma$. Then the inverse of the Green’s function, $\omega\hat{1}-H-\Sigma$, is decomposed as

[TABLE]

where $R=(\vec{\epsilon}_{1},\vec{\epsilon}_{2},\cdots)$ and $L=(\vec{e}_{1},\vec{e}_{2},\cdots)$ are the $4N\times 4N$ matrices of eigenvectors ($\epsilon_{i}$ and $e_{i}$ are right and left eigenvectors), and $\Lambda$ is a $4N\times 4N$ diagonal matrix with eigenvalues as its diagonal components. Using this decomposition, the Green’s function is given by

[TABLE]

Therefore, if the Hamiltonian and the self-energies are $\omega$ independent, then the arbitrary element of the Green’s function can be calculated efficiently.

In our formalism, in general, the self-energies are $\omega$ dependent. However, within the wide-band approximation and at zero temperature, the $\omega$ dependence of the self-energy appears as a step function, as shown in Eq. (25). Therefore, for a time-independent model, the above method is applicable by splitting the $\omega$ into two sectors, above and below the chemical potential. In the Floquet formalism, the number of sectors is given by $N_{s}+1$, where $N_{s}$ is the number of subbands. For the numerical integration in the calculation of the shift current using Eq. (34), we find that the cutoff of $\Lambda_{\omega}^{\pm}=\pm 4$ and discreteness of $\delta\omega=10^{-3}$ gives sufficient convergence.

In Fig. 7, we show an example of the calculation using the approximation introduced above. The figure shows the results of the light-position dependence of shift current for $l_{A}=39$ with different length of chain (the results for larger sizes are given in Fig. 3). The result for $N=100$ shows deviation from the results for larger $N$ with an arch-like structure. The results for $N\geq 200$ show almost the same magnitude of the photocurrent and it is nearly independent of $x_{A}$. The oscillation of the result around $x_{A}=0.5$, on the other hand, remains robust for larger sizes. Presumably, this is related to the competition between the typical distance between the states in a band and the broadening of the pole of Green’s function due to the leads. In general, the typical distance between the states is proportional to $1/N$; the numerical calculation converges when the typical distance become smaller than that of the broadening of the poles. In our setup, on the other hand, we only couple the edge sites to a lead. Hence, the poles of the Green’s function become sharper with increasing size deep inside the bulk. Due to the competition between the two energy scales, the finite size effect around $x_{A}$ persists up to larger sizes. In the main text, we mostly use $N=1200$ for calculation.

Appendix C Quasiclassical Analysis on the Carrier Density

This section elaborates the quasiclassical method discussed in Sec. 3.2. In the Green’s function theory developed in A, the current given in Eq. (34) is proportional to the change of the charge density from its “equilibrium” distribution, $f_{\alpha}(\omega)\rho_{i_{\alpha}}(\omega)$, where $\rho_{i_{\alpha}}(\omega)$ is the local density of states at site $i_{\alpha}$. Therefore, for the analysis on the $l_{A}$-independence, we focus on the carrier density in the conduction band. We consider a simple model that consists of a valence band, a conduction band, and the lead(s) as shown in Fig. 5. The carriers are transferred between the valence and conduction bands by $W_{vc}$ and between the conduction band and the lead by $\Gamma$. The carrier concentration of the conduction band is given by solving a differential equation:

[TABLE]

Here, $\rho_{c}$ ($\rho_{v}$) is the density of states for conduction (valence) band, $\delta n_{c}$ ($\delta n_{v}$) is the concentration of the electrons (holes) in the conduction (valence) band. The first and the second term represents the transition of electrons between the conduction and valence bands (we discuss the detail in the following paragraphs), and the last term represents the effluence of carriers to the lead(s). The form of the last term can be guessed intuitively, or formally derived from the Keldysh theory, which the current that flows out of the chain is given by the product of hybridization and the modulation of electron density from the equilibrium distribution.

To give a qualitative understanding, we estimate $W_{vc}$ using a perturbative approach. From the Fermi’s golden rule, for a perturbative Hamiltonian $H^{\prime}(t)=H^{\prime}(0)\sin(\Omega t)$, the transition amplitude from initial state $i$ to the final state $j$ reads

[TABLE]

where $\omega_{fi}\equiv\varepsilon_{f}-\varepsilon_{i}$ with $\varepsilon_{i}$ being the eigenenergy of $i$ th state. Here, we assumed the initial state $c_{i}(t)\sim 1$ and $c_{f}(t)\sim 0$ ($f\neq i$). In the $t\to\infty$ limit, the inside of the braces in this equation has two peaks in its absolute value, at $\omega=\pm\Omega$. We take the perturbative Hamiltonian to be the irradiation by light for electrons between sites $l_{L}$ and $l_{R}$ ($l_{A}=l_{R}-l_{L}$):

[TABLE]

By the Fourier transformation of $c_{n}$’s, $H_{A}$ reads

[TABLE]

with

[TABLE]

Here, $c_{Ak}$ ($c_{Bk}$) are the annihilation operators for electrons on the even (odd) sites with wave number $k$. For simplicity, we assumed $l_{L}$ to be an even number site and $l_{R}$ to be odd. For a given $l_{L}$ and $l_{R}$, $I_{n}(dk)$ is given by a periodic function with periodicity $\pi$ (which corresponds to the size of the Brillouin zone), and with a peak of hight $(n+1)/N$ and the full width at half maximum $\sim\frac{\pi}{(n+1)}$. As the eigenstates for the bulk limit with wave number $k$ is given by a superposition of $c_{Bk}$ and $c_{Ak}$, the change in the wave number while exciting the electrons from the valence state to the conduction state is limited by the width of $I_{n}(dk)$.

In the $t\to\infty$ limit, the transition from the valence band electron with momentum $k$ to the conduction band electron with $k^{\prime}$ occurs only between the $k$ and $k^{\prime}$ with energies precisely separated by $\Omega$. The transition amplitude for a state in the valence band with $k$ to that in a conduction band with $k^{\prime}$ is given by $\sim A^{2}|I_{n}(k-k^{\prime})|^{2}\sim(l_{A}+1)^{2}/4N^{2}$ ($l_{A}=l_{R}-l_{L}$). Due to the constraint from energy conservation, the transition occurs only between the state close to the crossing between $\varepsilon_{v}+\Omega$ and $\varepsilon_{c}$ [see Fig. 1(b)]; the width of the wave number the transition takes place is given by $\sim\pi/n+1=2\pi/l_{A}+1$. Finally, the density of electron states per a width of $k$ is given by $N/2\pi$. From these facts, $W_{vc}$ is estimated to be

[TABLE]

linearly proportional to $l_{A}$ and $A^{2}$.

In the steady state ($\delta\dot{n}_{c}=0$), the solution of Eq. (42) is given by

[TABLE]

Here, we used charge neutrality condition,

[TABLE]

Similarly, we get

[TABLE]