Time-dependent spectral properties of a photoexcited one-dimensional ionic Hubbard model: an exact diagonalization study

Motivated by the recent progress in time-resolved nonequilibrium spectroscopy in condensed matter, we study an optically excited one-dimensional ionic Hubbard model by exact diagonalization. The model is relevant to organic crystals, transition metal oxides, or ultracold atoms in optical lattices. We implement numerical pump-probe measurements to calculate time-dependent conductivity and single-particle spectral functions. In general, short optical excitation induces a metallic behavior imprinted as a Drude peak in conductivity or an in-gap density of states. In a Mott insulator, we find that the induced Drude peak oscillates at the pump frequency and its second harmonic. The former comes from the oscillation of currents, and the latter from the interference of single- and three-photon excited states. In a band insulator, the Drude peak oscillates only at the pump frequency, and quantities such as the double occupancy do not oscillate. The absence of the second harmonic oscillation is due to the degeneracy of multi-photon excited states. The in-gap density of states in both insulators correlates with the Drude weight and the energy absorption for weak pumping. Strong pumping leads to saturation of the in-gap density of states and to suppression of the Drude weight in the Mott regime. We have also checked that the above features are robust for insulators in the intermediate parameter range. Our study demonstrates the distinct natures of the multi-photon excited states in two different insulators.

In this work, we study an ionic Hubbard model, which can be considered as an effective model of a Holstein-Hubbard model where frozen phonons are dimerized to give static staggered potential modulation. n n n l l l . We focus on half-filling, where the number of particles = å N n l l is equal to the number of sites L. This model has been investigated to understand organic crystals or transition metal oxides [48][49][50][51]. It also has an intriguing topological property related to the Zak phase [52,53], which can be measured experimentally using ultracold atoms in an optical lattice [42].
When the Hubbard interaction U is dominant (U?J 0 , Δ), the system is a MI with a spin-density wave or antiferromagnetic ordering. When the Hubbard interaction U is negligible, the staggered potential Δ leads to a band dispersion where a gap of 2Δ opens at k=π/2. For a half-filled system, the lower band is completely occupied, and the upper band is empty in the ground state. Therefore, it is a BI.
In this work, we describe optical fields by a spatially uniform time-dependent vector potential A(t), and include it by Peierls' substitution, ; the hopping amplitude becomes complex. To investigate nonequilibrium dynamics under optical fields, we employ pump and probe pulses. For pump pulses, we consider the following Gaussian pulse where A 0 is the pump strength, t 0 the central pump time, and τ the full width at half maximum, and ω p the pump frequency. As a probe pulse (if included in the simulations), we use a step function at time t p , 4 p p probe 0 resulting in an electric field of a delta function E probe =E 0 δ(t−t p ) . This instantaneous form is taken for simplicity instead of narrow Gaussian forms, which are more realistic for experiments. The current density is given by

Method
We use an exact diagonalization method to study the ionic Hubbard model in (1). A ground state is found by the conventional Lanczos method [54], and then used as an initial state for the time-dependent Schrödinger equation. Time evolution is implemented by the Krylov subspace method based on the Lanczos method [36][37][38]. We use τ=5/J 0 and match the pump frequency ω p to the first peak of the linear absorption spectrum unless otherwise explicitly stated. The probe pulse strength E 0 is set to 0.005; we have checked that the results are independent of the probe strength for E 0 0.01. L=N=10 with the periodic boundary condition is used. In the following sections, we investigate time-resolved transient spectral properties of photoexcited states in the ionic Hubbard model by two protocols. One is time-resolved transient conductivity given by twodimensional pump-probe analysis. The other is a transient single-particle spectral function, which can be measured by time-resolved photoemission experiments. In the following, we briefly explain the simulation setups for the two protocols.

Time-resolved transient conductivity σ(t, ω)
In order to obtain time-resolved transient conductivity, we follow the formulation of [55] (see also [43][44][45][56][57][58][59].) A general conductivity response function ( ) , under a weak probe field E probe is defined by [60] Here δj (t) is a current solely induced by the probe field. When currents are nonzero even without a probe field (e.g. due to a pump field), such contributions are subtracted to obtain δj(t). Due to causality, , only depends on the time difference -¢ t t because of the temporal translation symmetry.
We consider a weak impulsive probe to the system at t p as (4). In real pump-probe experiments, a probe field has a finite duration. If its duration is short enough, the results are essentially the same as the one with the impulsive probe [45]. Repeating the pump-probe simulations with different time delay t p , we obtain the twodimensional current responses δj(t, t p ) in time domain Changing one of the time variables to s≡t−t p , we introducẽ , , This leads to˜( Finally Fourier transforming over s, we find In practice, we use a small damping factor η=1/L in the Fourier transformation to remove artifacts of the finite window.

Time-dependent transient single-particle spectral function A(t, ω)
In photoemission experiments, the central quantity of interest is the probability to detect an electron kicked out by photons. Assuming a uniform strength of probing laser pulses for duration of 2T, the probability of detecting an ejected electron after the probe pulse is expressed by lesser and greater Green's functions as [46,47,[61][62][63][64] , , where we have assumed that the coupling between incoming photons and electrons does not depend on their states. The first term corresponds to the photoemission spectroscopy (PES), and the second one to the inverse photoemission spectroscopy (IPES). When the Hamiltonian is time-independent, the expression is reduced to the conventional formula in the limit  ¥ T [61,62]. Here, instead, the time integral is limited by width 2T, which naturally gives a broadening of a spectral peak ∼π/T. We use T=10 in the following calculations.

Mott insulators
We first look at a MI, U=6 and Δ=0. In figure 1(a), we plot the energy expectation value = á ñ  H after photoexcitation with pump amplitude A 0 and pump frequency ω p . The energy is measured with the Hamiltonian without any vector potential. We see that around ω p ≈4.3 and 6.9, the system absorbs energies. This is consistent with the linear absorption spectrum α(ω) in figure 1 where |y ñ 0 is the ground state,  0 the ground state energy, and j the current operator in (5). For a fixed pump frequency ω p =4.3 (figure 1(c)), the energy absorption becomes nonlinear for A 0 0.3, even with a local minimum around A 0 ≈0.5. Such nonlinear dependence on A 0 is also observed in a charge-ordered state in an extended Hubbard model [67]. In figure 1(d), we plot the energy as a function of pulse duration τ. We see a large oscillation, even almost returning to the ground state around τ≈62; this is related to the coherent nature of excited states as discussed in section 5. In the next subsection, we see that such an oscillation is more apparent in a BI. The double occupancy, d º å á ñ n L l n do ,2 l also follows the same trend (not shown) as the energy absorption, which indicates that the energy absorption occurs via the doublon-holon excitation.
The real part of the time-dependent transient conductivity σ(t, ω) for pump strength A 0 =0.3 and A 0 =0.55 is plotted in figure 2(a). We see that the system becomes metallic due to the optical excitation, which is clearly seen as a Drude peak at low frequencies. Furthermore, a sub gap peak around ω≈1.9 appears, and the conductivity becomes negative slightly above the Drude peak. The sub-gap feature resembles the case of an extended Hubbard model [44], which is attributed to a transition from an odd-parity state to an even-parity state. However, various other features in σ(ω) depend non-monotonically on the pump strength A 0 as the energy absorption does. In particular, visible changes occur for the three peaks at ω≈0 (the Drude peak), ω≈4.3 (the first absorption peak), and ω≈6.9 (the second absorption peak), which are plotted in figures 2(b1) and (b2) as functions of A 0 . For small pump amplitude A 0 0.6, the Drude peak shows a similar structure as the energy absorption, e.g. the local minimum around A 0 ≈0.5. However, for stronger pump pulses A 0 0.6, too much energy is given to the system, and the Drude peak diminishes. The second peak at ω≈6.9 follows a similar trend. On the other hand, due to the resonant condition, the first excited peak at ω≈4.3 is depleted even to negative values when the energy absorption is small [68], and for larger A 0 0.7, the peak disappears.
We also observe that the conductivity oscillates in time, which is most evident for the Drude peak. In figures 3(a1) and (a2), the time evolution of the Drude peak D(t)=σ(t,0) and its Fourier spectrum, D(ω), are   plotted. This shows that the oscillation has exactly the pump frequency ω p and its second harmonic 2ω p . The fundamental mode comes from the oscillation of the current at ω p . On the other hand, the double occupancy n do oscillates only at 2ω p (figures 3(b1) and (b2)). Similarly, other quantities such as the kinetic energy and antiferromagnetic order oscillate at the second harmonic only. As we show in the next section using full diagonalization of a smaller system, such second harmonic oscillations come from the interference of the multiphoton excited states. The peaks in the Fourier spectra of the Drude weight and of the double occupancy are fitted by a simple Gaussian, and their A 0 dependence is plotted in figures 3(a3) and (b3). The second harmonic components, D(2ω p ) and n do (2ω p ), both have maximum around A 0 ≈0.6, and follow the same trend; a part of the Drude peak is governed by the doublon-holon excitation. On the other hand, the fundamental component D(ω p ) decreases significantly once the pump strength A 0 gets strong. Now we turn to the single-particle spectral functions A (t, ω) for a single-particle state α=(l=1, ). We focus on a final state at t f =t 0 +50, which is long after the pump pulse at t 0 . In appendix A, we show some other transient cases | | -< t t 50 0 . First, figure 4(a1) shows the PES signal A PES (ω) for various A 0 . The spectral weight is transferred from the lower Hubbard band to the upper Hubbard band, and the transferred weight as a function of A 0 is plotted in figure 4(a2). It shows a similar A 0 dependence as the energy absorption in figure 1. Interestingly, around A 0 ≈0.3, the weight is widely depleted in −3.5ω1.5 including the original Hubbard gap, in contrast to the case of A 0 =1.2 where only the slightly filled Hubbard gap remains. In figure 4(b1), we plot A(t f , ω) for various A 0 . We see that the gap is gradually filled with increasing A 0 , while the upper and lower Hubbard bands do not completely merge. In figure 4(b2), we show the integrated in-gap with Λ=3.0 and A(t f , 0)Λ as functions of A 0 . This again shows a similar behavior as the energy absorption. Therefore, we conclude that the in-gap density of states represents photoexcited carriers, which do not necessarily contribute to the coherent Drude peak, especially at a large pump strength A 0 .

Band insulators
As a BI, we take Δ=3, and U=0. Figure 5(a) shows the energy as a function of pump strength A 0 and frequency ω p . The optical excitation occurs at ω≈6.1 and 6.9 as shown in figure 5(b). As in the MI, the energy absorption depends non-monotonically on A 0 at a fixed ω p . We also find that the antiferromagnetic order is enhanced near the first excitation frequency, which is shown in appendix B. We expect that the spin part does not affect the charge dynamics that we are interested in due to charge-spin separation in one dimension, and thus we do not elaborate on it here. Figures 5(c) and (d) shows the A 0 and τ dependence of the energy. We see nearly perfect oscillations both in terms of A 0 and τ. Since the system has two bands separated by 2Δ, which resembles a two-level system, we consider these oscillations as Rabi-like oscillations. In the case of the MI, correlation effects hinder the simple two-level analogy (see figure 10), and deviations from a perfect oscillation occur especially at strong pump strength.
The time-dependent transient conductivity σ(t, ω) for various pump strength A 0 is plotted in figure 6(a). Contrary to the MI, there is no sub gap peak below the first excited peak ω≈6.1, nor the negative conductivity slightly above the Drude peak. Therefore, these features appearing in the MI are supposed to originate from the correlation effects [44]. The non-monotonic dependence of the conductivity peaks on A 0 is shown in figure 6(b). Up to the pump strength we studied A 0 1.2, we see that the Drude peak correlates with the energy absorption in figure 5. This is different from the MI, where the Drude peak does not follow the energy absorption for large A 0 . We expect that in the case of the BI the photodoped carriers independently contribute to the Drude peak, since there is no interaction. However, in the case of the MI, the photodoped carriers need to form coherent quasi particles to contribute to the Drude peak, which becomes difficult when the system absorbs too much energy. The peak at the pump frequency ω p changes signs due to the resonant condition, while the peak amplitude does not decay as increasing A 0 .
The oscillating features in σ(t, ω) in the BI are also different from the ones of the MI. In figure 7(a), we plot the time evolution and the Fourier spectrum of the Drude peak, D. We find that the Drude peak oscillates at the pump frequency ω p =6.1, which comes from the current oscillation. However, there is no oscillation at the second harmonic as in the MI. Similarly, the double occupancy plotted in figure 7(b) shows no oscillation. The lack of the second harmonic oscillations in the BI is an important difference compared to the MI. As we discuss in more detail in section 5, the absence is due to the structure of the excitation spectrum of the BI. In the BI, many multi-photon excited states exist near multiples of the energy gap 2Δ. Such degeneracy and the incoherent phases destroy the constructive interference.
Before we continue for the full-diagonalization results, we briefly comment on the single-particle spectral function A(t f , ω) for the BI. Here due to the doubling of the unit cell, we take the average of A α with α=( ) =  l 1, and ( ) =  l 2, . Furthermore, in the noninteracting limit, the sum of the photoemission and inverse photoemission in (11) is independent of the initial states (see appendix C). Therefore, we consider only the contribution from the photoemission part a A PES . Figure 8(a) shows the PES signals for different pump strength A 0 . We see that the upper band is more populated as the pump strength becomes larger. The integrated population in the upper band is plotted in figure 8(b). Here the trend follows the energy absorption closely as expected.

Intermediate cases
In order to investigate how the oscillations appearing in the Drude peak depend on the Hubbard interaction U and the staggered potential Δ, we consider intermediate cases here. For each parameter set, we set the pump frequency ω p again at the first absorption peak and keep the pump strength A 0 =0.1 with width τ=5 as before. We focus on the line U+2Δ=6. In the thermodynamic limit, there exists a quantum phase transition between a MI and a BI around U;2Δ?J 0 [69][70][71]. In figure 9(a), we plot linear absorption spectra α(ω) for different values of U. When U=0 (the BI case), there is only one peak at ω≈6.1. This is due to the fact that the BI has two bands separated by 2Δ. As U increases, the peak shifts towards lower energy. Near U≈3.7, the peak structure abruptly changes; the peak at ω≈0.9 jumps to ω≈2.8. This is related to a quantum phase transition, and due to the small system size, the critical value of U seems shifted from the ideal one U c ≈3 to 3.7. As U further increases, additional peaks appear at higher energies. At U=6 (the MI case), peaks are located at ω≈4.3, 6.9, and 9.1. When we pump at ω p =4.3, we excite the states at ω p and its higher harmonics.   Figure 9(b) shows the Fourier spectra of the Drude peaks, D(ω), for various interaction strength U and Δ. For 2.2U4.8, due to the closing of the energy gap, the system has initially a nonzero Drude peak, and the optical excitation reduces it. Since this is not the case that we are interested in, we exclude the data from the plot. We see that in the BI regime, there is no second harmonic peak, while in the MI regime, the second harmonic peak persists even with nonzero Δ. This finding confirms that the second harmonic oscillation of the Drude peak is the signature of the correlation in the MI, and it evidently distinguishes between MI and BI regimes.

Discussions
In this section, we elaborate on the properties of the photoexcited states by full-diagonalization of L=8 sites. All the energy eigenstates | ñ  n for n=1, L, 4900are calculated. We calculate excited states, |y ñ ex , by applying photoexcitation of A 0 =0.4 and τ=5 to the ground state. After the photoexcitation, the Hamiltonian becomes time-independent, and we can decompose the excited states by the eigenstates of the Hamiltonian | ñ  n . In figure 10, we plot the overlap of the photoexcited states with eigenstates, | | | y á ñ  n ex , in the MI and the BI. In both cases, there are non-negligible weights on the states whose energies are integer multiples of the pump frequency ω≈k ω p (k=0, 1, 2, L.). Here k=1 corresponds to the single photon absorption, which can be seen in the linear absorption spectrum. The states for k2 represent multi-photon excited (MPE) processes. Now we first consider the MI case. In the energy spectrum, there are many states between multiples of photon energies, kω p , due to the correlation U. However, the weight is highly concentrated in the non-  In both panels, lines are shifted by offsets depending on the corresponding parameters. The top line corresponds to U=0, Δ=3, and the bottom one to U=6, Δ=3. In between, the parameters change monotonically. In panel (b), data for 2.2U4.8are excluded, since there the system is initially metallic, and the optical excitation reduces the Drude peak. degenerate MPE states; for each ω;k ω p , there is only one dominant state. We label these states as | ñ k as in figure 10(a). The discreetness of the excited energy spectrum has been discussed as a result of charge-spin separation, which is a characteristic of one-dimensional systems [72][73][74]. The contribution from the five states ( | | | ) y å á ñ = k k 0 4 ex 2 accounts for more than 95% of the total weight. Therefore, the excited states can be well described by Then expectation value of an operatorÔ is approximated as The time-dependent oscillation comes from the off-diagonal elements, ¹ ¢ k k . In figure 11(a), the matrix elements of | | ¢ S kk for the double occupancy is plotted. There are non-negligible off-diagonal elements between k=1 and k=3, whose energy difference, 2ω p , gives the second harmonic oscillations. Similar structures exist for the antiferromagnetic order parameter and the kinetic energy.
Next we turn to the BI. In the energy spectrum, the MPE states are separated by 2Δ, and there are no states in between. The MPE states are quasi-degenerate in contrast to the MI, where there is only one state for each ω;k ω p . Thus, the distribution of the weight is not as concentrated as in the MI case. The largest six states give about 80% of the total weight, and we need thirty six states to have more than 95%. We label a state in the manifold of k-photon excited states as | ñ m k . The transition elements ¢ ¢ S m m , k k for the thirty six states are calculated. We find that there exist small amounts of off-diagonal elements between different manifolds, e.g. k=2 and k=3. However, we still need to sum over m k in (16). Due to the near degeneracy of MPE states in each manifold k, we can approximate (16) as   shows the matrix elements of |¯| ¢ S kk . In contrast to the MI case, the off-diagonal elements are ignorable. This is because the degenerate MPE states with different phases are summed up, and the constructive interference between different manifolds are destroyed. Thus, the second harmonic oscillation is absent.

Conclusion
In this work, a photoexcited one-dimensional ionic Hubbard model is studied by an exact diagonalization method. We find that the energy absorption in Mott and BIs exhibits nonlinear dependence on the pump amplitude. A nonzero Drude peak in conductivity or an in-gap density of states appear after short optical excitation indicating metallic behaviors, which can be measured by time-resolved experiments. In a Mott insulating regime, the Drude peak oscillates at the pump frequency and its second harmonic. The former originates from the current oscillation, and the latter from an interference between single-and three-photon excited states. We expect that as far as strong correlation exists, such coherent oscillations (may not be only the second harmonic) appear in other one-dimensional systems. Coherent oscillations also appear in the double occupancy and the kinetic energy, which might be easier to detect experimentally. In the band insulating regime, the Drude peak oscillates only at the pump frequency. This is because the superposition of the degenerate multiphoton excited states with different phases destroy the constructive interference. The relationship between the Drude weight and the in-gap density of states is also discussed. Including dynamical phonons is the next move to account for realistic experimental results [75,76]. Time-dependent spectroscopies on systems with modulated hopping [42,77,78] or multi-orbitals [79,80] are also interesting open problems.
In the end, the sum of the two contributions is reduced to which is independent of the initial distribution {p m }. This is similar to the case in equilibrium, where ( ) w a A is independent of temperatures for noninteracting systems [82].