Limitations of mean-field approximations in describing shift-current and injection-current in materials

We theoretically investigate bulk photovoltaic effects, with a specific focus on shift-current and injection-current. Initially, we perform a numerical analysis of the direct current (dc) induced by a laser pulse with a one-dimensional model, utilizing mean-field theories such as time-dependent Hartree--Fock and time-dependent Hartree methods. Our numerical results, obtained with mean-field theories, reveal that the dc component of the current exists even after irradiation with linearly polarized light as a second-order nonlinear effect, indicating the generation of injection current. Conversely, when we employ the independent particle approximation, no injection current is generated by linearly polarized light. To develop the microscopic understanding of injection current within the mean-field approximation, we further analyze the dc component of the current with the perturbation theory, employing the mean-field approximations, the independent-particle approximation, and the exact solution of the many-body Schr\"odinger equation. The perturbation analysis clarifies that the injection current induced by linearly polarized light under the mean-field approximations is an artifact caused by population imbalance, created through quantum interference from unphysical self-excitation pathways. Therefore, investigation of many-body effects on the bulk photovoltaic effects have to be carefully conducted in mean-field schemes due to potential contamination by unphysical dc current. Additionally, we perform the first-principles electron dynamics calculation for BaTiO$_3$ based on the time-dependent density functional theory, and we confirm that the above findings from the one-dimensional model calculation and the perturbation analysis apply to realistic systems.


I. INTRODUCTION
The bulk photovoltaic effect, the conversion of light into an electric current, underpins technological applications within modern society.Among the various mechanisms associated with the photovoltaic effect, the bulk photovoltaic effect has been intensively studied from both a fundamental and technological points of view [1][2][3].This effect, formerly referred to as the anomalous photovoltaic effect in inversion symmetry broken materials [4][5][6][7][8][9] and nowadays identified as the shift current [10], is a second-order nonlinear optical effect that generates a direct current (dc) in the presence of light.
Historically, theoretical examinations of the shift current encountered a divergence problem in the nonlinear susceptibility at the low-frequency limit.This issue was addressed by developing an explicit expression for the nonlinear susceptibility, using perturbation theory within the framework of the independent particle approximation [10].This expression, once coupled with first-principles electronic structure calculations, has been used extensively to analyze the shift current in real materials [11][12][13][14][15]. Concurrently, experimental research on the shift current has broadened, offering a variety of insights [16][17][18][19][20][21].
Despite the considerable interest in these phenomena, the majority of theoretical investigations into shift cur-rents in actual materials have been limited to the independent particle approximation.Recently, there have been several attempts to elucidate the role of many-body effects on the microscopic mechanisms of shift current.For example, Chan et al suggested substantial enhancement of the shift current due to excitonic effects, utilizing the GW plus Bethe-Salpeter approach [22].Furthermore, Kaneko et al proposed an enhancement of the photovoltaic effect in excitonic insulators via collective excitations [23].These studies indicate that many-body effects could have a pivotal role in enhancing photovoltaic effects.Nevertheless, the understanding of many-body effects in the microscopic mechanisms of the photovoltaic effects is still incomplete.For further development of detailed understanding of the many-body effects, it is crucial to assess the applicability of different theoretical approximations to accurately interpret these phenomena.
In this study, we conduct a numerical investigation of the bulk photovoltaic effects with a onedimensional model, employing mean-field theories such as the time-dependent Hartree-Fock (TDHF) and timedependent Hartree (TDH) methods, to assess the suitability of mean-field approximations.We also analyze the bulk photovoltaic effects through perturbation theory utilizing the independent particle approximation, mean-field approximation, and the exact many-body Schrödinger equation.Our results indicate that the photo-induced current under the mean-field approximation displays qualitatively different behaviors compared to both the independent particle approximation and the arXiv:2310.08875v1[cond-mat.mtrl-sci]13 Oct 2023 exact Schrödinger equation.Specifically, the mean-field approximation generates the dc current as the secondorder nonlinear effect even after irradiation of linearly polarized light, which neither the independent particle approximation nor the exact Schrodinger equation can describe.This suggests that the mean-field approximation could artificially induce the injection current under linearly polarized light.Consequently, analyses of the bulk photovoltaic effect using a mean-field approximation could significantly overestimate the effect due to this artifact, as the divergence of the susceptibility tensor of the injection current could completely overcome the susceptibility of the intrinsic shift current.Additionally, we perform the electron dynamics calculations for BaTiO 3 based on time-dependent density functional theory (TDDFT) to confirm the findings from the onedimensional model and the perturbation analysis apply to realistic systems.
The paper is organized as follows.In Sec.II, we describe theoretical and numerical methods for investigat-ing the bulk-photovoltaic effect.In Sec.III, we show the numerical results obtained by the method introduced in Sec.II and discuss the qualitative difference among the results obtained with different approaches.In Sec.IV, we analyze the bulk photovoltaic effect based on perturbation theory and explore the microscopic origin of the qualitative difference among different approximations.In Sec.V, we perform the first-principles electron dynamics calculation to analyze the shift-current and confirm the findings from the previous sessions for a realistic material.Finally, our findings are summarized in Sec.VI.

II. METHODS
In this study, the one-dimensional TDHF method is mainly utilized to simulate light-induced electron dynamics in solids.Each electronic orbital consisting a Slater determinant obeys the following TDHF equation, where b is the band index, k x is the Bloch wavenumber, and u b,kx (x, t) represents the periodic part of the Bloch wavefunction which satisfies u b,kx (x + a, t) = u b,kx (x, t), with the lattice constant, a.Here, A x (t) is a homogeneous vector potential related to an external electric field as E x (t) = −dA x (t)/dt.The ionic potential is denoted as v ion (x), and the spatial periodicity is imposed as v ion (x+a) = v ion (x).In this study, the ionic potential, v ion (x), is defined as where ρ ion (x) is the ionic charge density.Here, w(x) denotes the one-dimensional soft-Coulomb interaction and is defined as where β is a dimensionless adjustable parameter, and a 0 is the Bohr radius.Therefore, for this study, the softening parameter of the soft-Coulomb potential is set to a 0 .To introduce the standard Hartree potential v H (x, t) and the Fock operator vF (t), the one-body reduced density matrix and the one-body density are defined as where N k is the number of k-points in the calculation.
Using the one-body density, ρ(x, t), the Hartree potential is defined as Furthermore, the Fock operator is defined using the one-body density matrix as follows: By solving the TDHF equation, Eq. ( 1), we can eluci-date the electron dynamics induced by the vector poten-tial A x (t).Furthermore, by utilizing the time-dependent electron orbitals, u b,kx (x, t), we can evaluate physical quantities within the time domain.For instance, the induced electric current can be calculated as follows: where the velocity operator vx (t) is defined as: To model a solid-state system with broken inversion symmetry, we employ the following ionic charge distribution, denoted as ρ ion (x): In this study, the lattice constant a is set to 4.0352 Å, analogous to that of BaTiO 3 [24], a material typically utilized for investigating shift-current.To impose charge neutrality, we incorporate four electrons within the unit cell for this study.Furthermore, we set the dimensionless parameter β of the soft-Coulomb interaction to 0.2915 a.u.so as to reproduce the band gap of BaTiO 3 with the model calculation.
Figure 1 shows the band structure obtained by solving the static Hartree-Fock equation under the conditions stated above.The red solid line denotes the two filled bands (valence bands), while the blue solid line illustrates the empty band (conduction bands).The calculated band-gap at the Γ-point (k x = 0) is 3.2 eV.Computed bandstructure of the one-dimensional model of a bulk with broken inversion symmetry with the static Hartree-Fock method.The two spin-degenerate occupied bands are described as the red solid line, whereas the conduction bands are described as the blue solid line.
For the subsequent calculations in Section III, the realspace coordinate x within the unit cell (0 ≤ x ≤ a) is dis-cretized into 16 grid points.Similarly, the First Brillouin zone is divided into 1025 k-points.

A. Linear Optical Properties
To develop insights into optical responses of solids, we first revisit linear optical properties of the system with mean-field approximations.For this purpose, we calculate the electron dynamics in the presence of an impulsive distortion given by where Θ(t) is the Heaviside step function, and A 0 is the amplitude of the distortion.In this study, we set A 0 to 10 −4 a.u.so that the induced current is proportionate to the amplitude A 0 , resulting in the linear response.
The TDHF equation, Eq. ( 1), is solved with the vector potential as defined in Eq. (11).The static Hartree-Fock method is utilized to compute the ground state of the system, which is subsequently used as the initial condition for Eq. ( 1).The induced current, J(t), is computed using Eq.(8).
Assuming that the amplitude A 0 is sufficiently small, the optical conductivity σ(ω) of the system can be evaluated as where T sim represents the simulation time period, and W (x) is a window function introduced to suppress numerical noise in the Fourier transform resulting from finite simulation time.For practical calculations to analyze the linear responses, the simulation time, T sim , is set to 20 fs.Furthermore, we employ the following form of the window function in this study: Figure 2 shows the real-part of the conductivity of the one-dimensional solid-state system, computed with the parameters described in Sec.II.The red solid line represents the results obtained by using the TDHF method, revealing a distinctive peak structure below the band gap of 3.2 eV.This peak can be understood as the excitation to an excitonic state.
To obtain insights into the characteristics of the meanfield approximation within the context of the linear optical response, we introduce two distinct approximations to the TDHF method.The first approximation is the independent-particle approximation, denoted as TDHF-IP, which freezes the time-dependence of the Hartree potential and Fock operator in the TDHF method.The second approximation is the time-dependent Hartree approximation, denoted as TDH.Real-part of conductivity σ(ω) of the onedimensional solid-state model.The result obtained using the TDHF method is described as red solid line, the result obtained using the independent particle approximation is described as blue dashed line, and the result obtained using the TDH method is described as green dotted line.
When employing the TDHF-IP approach, the static Hartree-Fock method is employed to prepare wavefunctions and eigenvalues as the initial conditions for the time propagation.In this approximation, we neglect the time dependencies of v H (x, t) and vF (t) by replacing the timedependent one-body density matrix, ρ DM (x, x ′ , t), with ρ DM (x, x ′ , 0) for the computation of v H (x, t) and vF (t).Therefore, within the TDHF-IP approximation, we disregard the time dependencies of the mean-field potential completely.The resulting dynamics can be interpreted as those of an independent particle system under a corresponding external potential.
By contrast, the TDH method removes the Fock operator, vF (t), from the TDHF process during both the preparation of the ground state and the time propagation of the system.Consequently, the mean-field potential consists solely of the Hartree potential, v H (x, t).This approximation provides one of the simplest mean-field approximations for quantum many-body systems.
In Fig. 2, the blue dashed line represents the result obtained from the TDHF-IP method, while the green dotted line denotes the result obtained from the TDH method.The TDHF-IP result indicates the vanishing of the excitonic peak below the gap, with the spectral weight of the excitonic peak absorbed by the above-gap absorption, thus reflecting a nature of independent particle systems -the absence of excitons.The TDH method yields a spectral structure similar to that of TDHF-IP but with a substantial red-shift.Such similar spectral structures imply that the TDH method does not capture the excitonic contribution to the excitation spectrum.The red-shift of the spectrum in the TDH method reflects the band gap reduction caused by the exclusion of the exchange interaction.These insights confirm the widely accepted understanding in condensed matter physics: the static contribution of the Fock operator enlarges the gap, while the dynamical contribution of the Fock operator describes the excitonic response of the lowest optical transition [25].Employing these three methods, we further analyze nonlinear optical responses in the subsequent section.

B. Second-Order Nonlinear Optical Responses: Shift Current and Injection Current
Having revisited the mean-field characteristics within linear optical responses, we extend our investigation to nonlinear optical phenomena, with a specific focus on the second-order nonlinear optical effects, namely shift current and injection current [10].To obtain insights into these second-order nonlinear optical phenomena, we first review the nonlinear susceptibility in the frequency domain, and the corresponding current dynamics in the time domain.
For a deeper understanding of these nonlinear optical phenomena, a time-domain behavior of induced responses complements their divergent behavior of the susceptibilities in the frequency domain.We conduct this exploration by analyzing the dynamics induced by a laser pulse represented as: where, f (t) is the envelope function of the laser pulse, and ω 0 is the average frequency.For this analysis, we assume that the envelope function f (t) exhibits slow temporal variation, and that the Fourier transform of Eq. ( 17) is concentrated predominantly around ω = ω 0 .We initiate our analysis with the polarization, P rec (t), associated with optical rectification.The polarization induced by the pulse of Eq. ( 17) can be evaluated as follows: Here, we adopted the assumption that the Fourier transform of E(ω) is predominantly localized around ω = ω 0 .By isolating the low frequency component of Eq. ( 18), we extract the dc-like component of P The results demonstrate that the polarization associated with optical rectification is directly proportional to the square of the envelope function, f (t).Similar analysis can be performed for the current J (2) sft,dc (t) associated with the shift current, and the acceleration K (2) inj,dc (t) associated with the injection current.This investigation yields the following relation: Figure 3 shows a graphical representation of the timedomain behaviors of optical rectification, shift current, and injection current as described by Eqs.(19)(20)(21).The time profile of a sample laser pulse is depicted in Fig. 3 (a), with the envelope function of the laser pulse represented by a black dashed line.The time profile of the polarization for optical rectification and shift current is shown in Fig. 3 (b).The polarization associated with optical rectification shifts only during laser irradiation, whereas the polarization associated with shift current remains finite even after the laser field ends.Figure 3 (c) shows the time profile of the currents associated with shift and injection currents.The shift current is only induced during laser irradiation, while the injection current remains finite even after the laser field ends.This behavior of the injection current can be understood by the fact that the acceleration, K (2) inj,dc (t), associated with the injection current shifts only during laser irradiation, as displayed in Fig. 3 (d).
Guided by the time-domain behavior of the secondorder nonlinear optical responses as shown in Fig. 3, we then numerically examine the impact of mean-field approximations on these nonlinear responses by employing the TDHF method.
For practical calculations, we compute the electron dynamics induced by the following vector potential: in the domain −T pulse /2 ≤ t ≤ T pulse /2, and zero outside this domain.
Here, E 0 is the peak field strength, ω 0 is the mean frequency, ϕ CEP is the carrier envelope phase (CEP), and T pulse is the pulse duration.In this work, we set E 0 to 2 × 10 −4 a.u., ω to 3.3 eV/ℏ, and T pulse to 40 fs.We determined the field strength E 0 such that the resulting dc-like current is dominated by second-order nonlinear optical responses.Furthermore, we choose the photon energy ℏω 0 to exceed the band gap, fulfilling the condition to induce the shift current [10].We treat ϕ CEP as a tunable parameter to extract the dc-like component of the induced current.We denote the current induced by the laser field of Eq. ( 22) as J x (x, ϕ CEP ), explicitly noting ϕ CEP dependence.To extract dc-like component of the induced current, we consider the CEP average of the induced current as By utilizing this average, we effectively eliminate highfrequency components in the induced current, thereby extracting the dc-like component.For practical analysis, we calculate the integral of Eq. ( 23) as the mean of four values of ϕ CEP : ϕ CEP = 0, π/2, π, 3π/2.
Figure 4 shows the time-evolution of the dc-like component of the induced current, J x,dc (t), calculated using Eq.(23).As shown in Fig. 4, there are distinct differences in the behaviors of the results obtained by the TDHF, TDH, and TDHF-IP methods.Both TDHF and TDH methods show a finite current even after the laser field ends.By contrast, the TDHF-IP method does not present any residual current.These observations suggest that mean-field theories, such as TDHF and TDH, produce qualitatively distinct dc-like second-order nonlinear current when compared to the independent-particle approximation.
Importantly, the qualitative discrepancy in the optical response behavior between the mean-field theories and the independent-particle approximation is not due to excitonic contributions, since the TDH method does not incorporate excitonic effects.This discrepancy emerges solely from the time-dependence of mean fields.The slow component of the second order induced current, J x,dc (t) is shown.The result obtained using the TDHF method is shown as the red solid line, and that obtained using the independent particle approximation is shown as the blue dashed line, and that obtained using the TDH method is shown as the green dotted line.As a guide, the envelope function of the applied electric field of Eq. ( 22) is shown as black dash-dot line.
Based on the classification of nonlinear responses outlined in Eqs.(19)(20)(21), a finite dc-like current after laser fields end implies the presence of the injection current.However, it has been extensively discussed that the injection current originates from the breaking of time-reversal symmetry [10,27], such as in the irradiation of circularly polarized light.The residual current depicted in Fig. 4, resultant from the irradiation with linearly polarized light, is expected to be an artifact, connected to the unphysical divergence of the response function.
A few decades ago, there was considerable discussion regarding the presence of unphysical divergences in second-order nonlinear optical responses in materials, and it has been suggested that these divergences are caused by numerical errors and could be removed by employing a sum-rule suitable for semiconducting systems within the independent particle approximation [28].Another theoretical method involves the evaluation of second-order nonlinear optical responses utilizing an explicit susceptibility expression derived from perturbation theory within the independent particle approximation [10].This perturbative approach successfully avoids the unphysical divergence as it analytically manages the susceptibility singularity.
Even though second-order nonlinear optical phenomena hold considerable importance, studies exploring the divergent behavior of the response function for manybody systems beyond the independent particle approximation have been limited [29].In the following section, we examine the effects of many-body and mean-field effects on the second-order nonlinear optical phenomena based on the perturbation theory.

IV. PERTURBATION ANALYSIS
In the previous section we demonstrated numerically that mean-field approximations, including TDHF and TDH methods, may induce unphysical dc-like current as a result of the second-order nonlinear optical process with linearly polarized light.In this section, we investigate the origin of this unphysical current within perturbation theory.We analyze the current induced by a laser pulse and evaluate dc-like component of the induced current after laser irradiation.This perturbation analysis is performed at the exact many-body Schrödinger equation, the independent particle approximation, and the meanfield approximations.

A. Exact Schrödinger Equation
We first consider the light-induced current within the exact Schrödinger equation for a many-electron system, where Ĥ(t) denotes the many-electron Hamiltonian with a one-body potential v(r), and an interacting potential w(r −r ′ ).Here, the time-dependent vector potential A(t) is included to describe external electric fields as E(t) = − Ȧ(t).For practical analysis with a finite laser pulse, we impose that the vector potential vanishes for t > t f : A(t > t f ) = 0.
The time-dependent Hamiltonian for a many-electron system, Ĥ(t), can be decomposed into the following three components: where N e is the number of electrons in the system, and P is the total momentum of the system given by P = j p j .
For later convenience, we introduce a wavefunction | Ψ(t)⟩ through a unitary transformation as | Ψ(t)⟩ = e i t dt ′ C(t ′ )/ℏ |Ψ(t)⟩.Utilizing this, the many-electron Schrödinger equation, Eq. ( 24), can be rewritten as To investigate second-order nonlinear responses, we begin by introducing a perturbative expansion of the wavefunction | Ψ(t)⟩ up to the second order of perturbation [30]: where |Φ 0 ⟩ refers to the ground state wavefunction of the unperturbed Hamiltonian, Ĥ0 , and E 0 represents its associated ground state energy, fulfilling Ĥ0 The first-and second-order wavefunctions are denoted as |δΨ (1) (t)⟩ and |δΨ (2) (t)⟩, respectively.The correspond-ing first-and second-order dynamical phase factors are determined by the first-order energy shift, E (1) (t), and the second-order energy shifts, E (2) (t), respectively.By substituting Eq. ( 31) into the modified timedependent Schrödinger equation, Eq. ( 30), we derive the following relation for each order: To proceed with the analysis, we introduce the eigenstates of the unperturbed Hamiltonian, Ĥ0 , as follows: If there is a set of degenerate eigenstates with respect to Eq. ( 34), we choose to define the eigenstates such that at least one of the interested Cartesian components of P is diagonalized within the subspace spanned by these degenerate eigenstates.
By utilizing these eigenstates, the perturbative wavefunctions, |δΨ (1) (t)⟩ and |δΨ (2) (t)⟩, can be expanded as |δΨ (1) where Ω a is defined as Ω a = (E a − E 0 )/ℏ.Here, C a (t) and C (2) a (t) represent the expansion coefficients for the first-and second-order wavefunctions, respectively.The expansion excludes the unperturbed ground state |Φ 0 ⟩, given that the energy shifts, E (1) (t) and E (2) (t), are defined as We proceed by substituting Eq. ( 35) and Eq.(36) into Eq.(32) and Eq. ( 33), respectively.Consequently, the derived equations of motion for the expansion coefficients are as follows: From these derived expressions, it is evident that the coefficients C (1) a (t) and C (2) a (t) become time-invariant once the perturbation ends ( V (t) = A(t) = 0 for t > t f ).By using the above perturbative expansion, the second order nonlinear current can be expressed as Given that C (1) a (t) and C (2) a (t) become constant after the perturbation ends ( V (t > t f ) = A(t > t f ) = 0), we can evaluate the dc-component of the current J (2) (t) after the fields end as In order to analyze the dc component of the secondorder current, we turn our attention to the first-order coefficient, C (1) a (t), at time t = t f .This is calculated by integrating Eq. (39) as follows: where Ã(ω) denotes the Fourier transform of A(t).Here, we have exploited the property that A(t) = 0 for t > t f .By employing this explicit expression of C a (t f ), the dc component of the current J (2) dc in Eq. ( 42) can be evaluated as To further analyze the dc-component of the secondorder current, we next explore a characteristic of the eigenstates |Φ a ⟩.The real-space representation of a many-body state |Φ a ⟩ is defined as Φ a (r 1 , • • • , r Ne ) = ⟨r|Φ a ⟩.As a result, the static Schrödinger equation with the unperturbed Hamiltonian can be reexpressed as follows: It can be readily verified that the complex conjugate of an eigenstate, Φ * a (r Moreover, by assuming that the ground state has the time-reversal symmetry as |Φ (−) 0 ⟩ = |Φ 0 ⟩, matrix elements of the total momentum operator hold the following relation: By employing Eq. ( 46) and Eq. ( 47), the dc-component of the second-order nonlinear current in Eq. ( 44) can be further evaluated with the time-reversed states as The final expression for the dc-component of the current, Eq. ( 48), suggests that residual dc-component cannot be induced solely by linearly polarized light as a second-order nonlinear effect.This can be straightforwardly confirmed by evaluating Eq. ( 48) by employing a linearly polarized vector potential in the frequency domain, Ã(ω) = Ã(ω)e p , with a pure-real unit vector along the polarization direction, resulting in J (2) dc = 0.By contrast, if we evaluate Eq. ( 48) with elliptically polarized light represented as Ã(ω) = e x Ãx (ω) + e y Ãy (ω), the dc-component is given by This result indicates that the dc-component of the current may be induced only if the light field is elliptically polarized ( Ãx (Ω a ) Ã * y (Ω a ) − Ã * x (Ω a ) Ãy (Ω a ) ̸ = 0), or namely breaking the time-reversal symmetry.This current represents the injection current induced by a timereversal symmetry broken field in a system lacking inversion symmetry.
Working with the many-body Schrödinger equation, we have confirmed that the second-order dc current cannot exist after a laser pulse ends if the laser field is linearly polarized in a time-reversal symmetry system.In other words, the injection current cannot be induced by a linearly polarized light, but it can be induced only by an elliptically polarized light.This conclusion is consistent with previous works [10,27].

B. Mean-field theory
To extend the perturbation analysis with the exact many-body Schrödinger equation, here we analyze the dc component of the second-order nonlinear current by employing mean-field theories.As a practical mean-field approximation to solid-state systems, we assume that a many-electron system is described by a single Slater determinant, and each electronic orbital is given by a Bloch state, ψ bk (r, t) = e i(k+A(t)/ℏ)•r u bk (r, t), where u bk (r, t) = u bk (r + a j , t) stands for the periodic part of the Bloch function and a j is any lattice vectors.Furthermore, we impose that electronic orbitals, ψ bk (r, t), obey the following mean-field equation of motion: where v(t) denotes a mean-field potential operator.We assume that the Hamiltonian ĥ(t) has spatial periodicity with the lattice vectors, a j .Both the TDHF and TDH methods utilized in the numerical analysis in Section III are described by the same form as Eq. ( 50).Furthermore, it is worth noting that the time-dependent Kohn-Sham equation in the timedependent density functional theory is also described by the same form.
In order to proceed with the perturbation analysis, we introduce the eigenstates of the unperturbed Hamiltonian ĥ0 as follows: where ĥ0 is the unperturbed Hamiltonian, ϕ bk (r) represents the eigenstates of ĥ0 , and v0 is the non-perturbed potential.According to the discussion for the perturbation analysis with the exact many-body Schrödinger equation, if there is a set of degenerate eigenstates, we choose to define the eigenstates such that at least one of the Cartesian components of p is diagonalized within the subspace spanned by these degenerate eigenstates.Next, we expand the time-dependent Hamiltonian ĥ(t) up to the second order in terms of the external field A(t) as follows: ĥ(t) = ĥ0 + ĥ(1) (t) + ĥ(2) (t), ( 53) where δv (1) (t) and δv (2) (t) represent the first-and second-order contributions from the mean-field potential v(t), respectively.

Independent particle approximation
To highlight contributions from the time-dependent mean-field to the dc-component of nonlinear current, we first revisit the results of the independent particle approximation, obtained by setting δv (1) (t) = δv (2) (t) = 0.Under these constraints, the coefficients c (1) a,bk (t) and c (2) a,bk (t) remain time-invariant after the laser irradiation since both ĥ(1) (t) and ĥ(2) (t) become zero.For practical analysis, we assume A(t) = 0 for t > t f .Consequently, the first-order coefficient can be expressed as Next, we calculate the second-order nonlinear current associated with the orbital ψ bk (r, t) as follows: bk (r, t) As discussed in Sec.IV A, we can determine the dccomponent of the second-order nonlinear current after the laser field ends using the following expression: a,bk (t > t f ) Assuming that the unperturbed mean-field Hamiltonian ĥ0 has time-reversal symmetry, the Bloch states at k and −k exhibit the following time-reversal relations: ϕ b,−k (r) = ϕ * bk (r) and ϵ b,−k = ϵ bk .By utilizing these relations, the sum of the dc-component of the current at k and −k can be evaluated as Finally, the dc-component of the total current can be evaluated as where the sum of the index b is taken only for the occupied orbitals (occ), and Ω BZ is the volume of the Brillouin zone.
From the final expression in Eq. (67), it is evident that the dc-component of the current J (2) dc vanishes when linearly polarized light is considered.For instance, assuming a vector potential of the form Ã(ω) = Ã(ω)e p , J dc vanishes due to the integrand of the last line of Eq. (67) becoming zero.Instead, under elliptically polarized light, the dc-component of the current may remain finite.If we consider the vector potential to be of the form Ã(ω) = Ãx (ω)e x + Ãy (ω)e y , the dc-component of the current in Eq. ( 67) is given by J In contrast to the current under linearly polarized light, Eq. ( 68) shows that the dc-component of the current may remain finite under elliptically polarized light, indicating the possibility of inducing an injection current.Moreover, the dc-component of the current in Eq. ( 68) arises from the interference between excited states induced by Ãx (ω)e x and Ãy (ω)e y .This is nothing but the quantum interference among two excitation paths associated with orthogonal components of light fields [10,16] 2. A mean-field approximation with linearly polarized light Here, we investigate the impact of a mean-field contribution on the dc-component of the total induced current.To achieve this, we incorporate the contributions from v(1) (t) and v(2) (t) into the perturbation analysis, employing the independent-particle approximation, as discussed in Sec.IV B 1.
To specifically examine the mean-field contribution only during laser irradiation, we neglect the induced mean-field effect after the laser fields ends.For practical analysis, we impose v(1) (t) = v(2) (t) = 0 as well as A(t) = 0 for t > t f .Under these constraints, we can evaluate the dc-component of the current after the laser fields end as follows: where δn ak represents the population imbalance between k and −k and is defined as: By integrating Eq. ( 61) with the mean-field contribution δv (1) (t), we can evaluate the first-order coefficient as follows: where ṽ(1) (ω) represents the Fourier transform of δv (1) (t).By substituting Eq. (71) into Eq.( 70), the population imbalance can be rewritten using the following decomposition: Each decomposed component is given as follows: In the above expressions, δn A a,bk represents the population imbalance induced by direct excitation from laser fields A(t).The imbalance δn B a,bk arises from interference between quantum states excited by both the laser fields A(t) and the induced mean-field v(1) (t).Finally, δn C a,bk is the imbalance generated solely by the mean-field contribution v(1) (t).
Consistently with the exact many-body Schrödinger equation and the independent particle approximation, the population imbalance δn A a,bk becomes zero if the external field A(t) is linearly polarized light, but it may remain finite only if A(t) is elliptically polarized light.Conversely, δn B a,bk and δn C a,bk may remain finite even under the application of linearly polarized light.The induced mean-field enables additional excitation pathways, leading to quantum interference, which results in the pop-ulation imbalance and subsequent dc current even after laser irradiation ends.This behavior is consistent with the results of our numerical simulations with the meanfield approximations (TDHF and TDH) shown in Fig. 4, which also demonstrate dc current remains finite even after the laser fields end.Notably, the independent particle approximation does not show any dc current after the laser fields end in Fig. 4.
We emphasize that the resultant dc current following the irradiation by linearly polarized light is an artifact of the mean-field approximation.Such a dc current cannot be induced if we consider the exact manybody Schrödinger equation (see discussion in Sec.IV A).This artifact manifests as an unphysical divergence in the nonlinear susceptibility in the low-frequency limit.The unphysical divergence may overcome the intrinsic low-frequency response of the system, inhibiting a proper investigation of the photovoltaic effect, such as shift current, within the scope of mean-field theories unless the artifact is entirely eliminated.
The perturbation analysis suggests that the unphysical dc current originates from excitation paths opened up by the induced mean-field dynamics.Therefore, the artifact is a result of self-excitation via the mean-field that are not present in the fully interacting many-body solutoin (see Sec. IV A).One could potentially improve the mean-field description for the photovoltaic effect by eliminating the unphysical self-excitation effect.This can be achieved by appropriately designing the Hartree-exchange-correlation kernel f Hxc (r, r ′ , ω) in the time-dependent density functional theory.

V. FIRST-PRINCIPLES ANALYSIS
Having established the understanding on the limitation of the mean-field approximation resulting in the unphysical dc current based on the one-dimensional model simulation, we then quantify this limitation for more realistic situations by using the first-principles calculations based on the TDDFT.As an example of realistic systems, we take BaTiO 3 .For the investigation on the nonlinear current, we first compute the ground state of the tetragonal phase of BaTiO 3 with the structural parameters at 300 K [24] by using the local density approximation [31].We note that the polarization direction of the system is the c-axis.
Once the ground state of the tetragonal BaTiO 3 is prepared with the above conditions and parameters, we then compute the light-induced electron dynamics by solving the time-dependent Kohn-Sham equation with the adiabatic local density approximation (ALDA), where the exchange-correlation potential is evaluated by the local density approximation with the instantaneous electron density at each time.As an external field, we employ a vector potential polarized along c-axis with the form of Eq. (22).For practical calculations of BaTiO 3 , we set E 0 to 2.75 MV/cm, ω 0 to 4 eV/ℏ, and T pulse to 20 fs.We repeat the TDDFT calculations with the external fields by using four different CEPs, ϕ CEP = 0, π/2, π, and 3π/2.By averaging the obtained current from these calculations, we extract the dc-like current component of the light-induced current with Eq. (23).
Figure 5 shows the extracted dc-like current in BaTiO 3 as a function of time.As a reference, the envelope function of the applied laser field is also shown as the black dash-dot line.In Fig. 5, the result obtained by using the ALDA is shown as the red solid line, while that obtained by using the independent particle approximation is shown as the blue dashed line.Consistently with the one-dimensional model simulations and the perturbation analysis, the finite dc-like current remains even after the irradiation of linearly polarized light when the mean-field potential, or the Kohn-Sham potential, has the time dependence with the ALDA.By contrast, the unphysical constant current after the field irradiation is removed when the time-dependence of the mean-field is ignored in the independent particle approximation.As seen from Fig. 5, one can confirm that the unphysical dc-like due to the mean-field approximation significantly affects the photovoltaic response even in a realistic material, BaTiO 3 , beyond the one-dimensional model calculations.Therefore, the investigation on the shiftcurrent and injection-current with mean-field theories has to be carefully conducted by properly removing the unphysical dc-current response.

VI. SUMMARY
In this study, we theoretically investigated secondorder nonlinear optical effects in solids, specifically focusing on the shift and the injection currents.We utilized numerical simulations based on mean-field theories, such as TDHF and TDH methods, to explore real-time electron dynamics under laser pulse excitation and the resulting nonlinear current in the time domain.The numerical simulations demonstrated that the dc-component of the second-order nonlinear current remains finite even after irradiation with linearly polarized light, indicating a possible induction of the injection current by the timedependent mean field.However, simulations based on the independent particle approximation did not exhibit such a residual dc-component after linearly polarized light irradiation.
To understand the origin of this residual dc-component observed in the nonlinear current, we performed perturbation analysis using various levels of theories, including the exact many-body Schrödinger equation, mean-field approximations, and the independent particle approximation.The perturbation analysis with the exact manybody Schrödinger equation revealed that linearly polarized light cannot induce such residual dc-component in the second-order nonlinear current when the system has time-reversal symmetry before laser irradiation.Similarly, the perturbation analysis using the independent particle approximation arrived at the same conclusion, indicating that the presence of a residual dc-component in the second-order nonlinear current after laser irradiation may occur only if the applied field is elliptically polarized or breaks time-reversal symmetry.Consequently, the perturbation analysis clarified that the residual dccomponent in the nonlinear current is an artifact of the mean-field approximations.
Further we performed perturbation analysis within the mean-field theories, revealing that the unphysical dccomponent in the nonlinear current arises from population imbalances in the Brillouin zone, specifically at k and −k points.Additionally, the perturbation analysis showed that this unphysical population imbalance is caused by quantum interference between different excitation paths involving self-excitation paths opened via the time-dependent mean field that are not present in the full many-body solution.Furthermore, we performed the first-principles electron dynamics calculations based on TDDFT using the adiabatic local density approximation and confirmed that these findings apply also to realistic materials.The resulting residual current obtained by using the adiabatic approximation indicates a significance of a time-nonlocal memory effect in the TDDFT to capture the proper nonlinear dynamics of interacting many-body systems.
The unphysical dc current, resembling the injection current, induced by mean-field approximations may overcome the intrinsic shift-current contribution due to the higher susceptibility divergence of injection current compared to that of the shift current.
From a different perspective, the induction of unphysical current through self-excitation paths via mean fields suggests opportunities for improving density and currentbased many-body theories to describe light-induced nonlinear phenomena more accurately.An accurate theory should prevent the induction of unphysical dc current by eliminating population imbalances arising from self-excitation paths.This could be achieved by properly designing the Hartree-exchange-correlation kernel, f Hxc (r, r ′ , ω), in time-dependent density functional theory.While this study focused on second-order nonlinear optical responses, it would be crucial to consider the potential significant impact of self-excitation paths in even higher-order nonlinear phenomena as well.Therefore, the limitations of local and semi-local adiabatic approximations and the effects of unphysical self-excitation paths should be carefully evaluated for further investigations on highly-nonlinear optical phenomena.
FIG. 1.Computed bandstructure of the one-dimensional model of a bulk with broken inversion symmetry with the static Hartree-Fock method.The two spin-degenerate occupied bands are described as the red solid line, whereas the conduction bands are described as the blue solid line.

FIG. 3 .
FIG. 3. Schematics of the second-order nonlinear optical responses.(a) The time-profile of an applied example pulse is shown.(b) The polarizations associated with the optical rectification (red) and the shift current (green) are shown as a function of time.(c) The induced currents associated with the shift current (green) and the injection current (blue) are shown.(d) The acceleration associated with the injection current is shown.
FIG.4.The slow component of the second order induced current, J x,dc (t) is shown.The result obtained using the TDHF method is shown as the red solid line, and that obtained using the independent particle approximation is shown as the blue dashed line, and that obtained using the TDH method is shown as the green dotted line.As a guide, the envelope function of the applied electric field of Eq. (22) is shown as black dash-dot line.

FIG. 5 .
FIG.5.The dc-like current in BaTiO3 induced by a linearly polarized laser pulse.The results are computed by using the TDDFT with the ALDA (red solid line) and with the independent particle approximation (blue dashed line).
1 , • • • , r Ne ), satisfies Eq. (45), indicating that Φ * a (r 1 , • • • , r Ne )is also an eigenstate of the Hamiltonian.We introduce a ket vector, |Φ ⟩ is the time-reversed state of |Φ a ⟩, and it may or may not be identical to the original state, |Φ a ⟩. a