Impact ionization and multiple photon absorptions in the two-dimensional photoexcited Hubbard model

We study the non-equilibrium response of a 4x3 Hubbard model at U=8 under the influence of a short electric field pulse, with the main focus on multiple photon excitations and on the change of double occupancy after the pulse. The behavior mainly depends on the driving frequency of the electric field. The largest change of double occupancy occurs during the pulse. For frequencies below the Mott gap, we observe multiphoton excitations at large field intensities. For frequencies beyond the gap energy, there is a region where Auger recombination reduces the double occupancy after the pulse. Impact ionization (Multi Exciton Generation), namely a growing double occupancy after the pulse, occurs for frequencies larger than twice the Mott gap. From the Loschmidt amplitude we compute the eigenstate spectrum of the quantum state after the pulse, observing multiple distinct photon excitation peaks, in line with expectations from a quasiparticle picture. We introduce a technique with which we analyze the time evolution of double occupancy in each peak individually. The long-term behavior of the double occupancy almost only depends on the absorbed energy, and we explore the connection of this property to the Eigenstate Thermalization Hypothesis.


I. INTRODUCTION
One exciting area of research is the influence of photoexcitations on strongly correlated electron systems. In strongly correlated materials, one generally cannot neglect the electron-electron interaction when describing and predicting material properties, because strong localization and Coulomb interaction effects play an important role 1 . In the context of computational material science this means that one of its most prominent approaches, the Density Function Theory (DFT) 2 , fails to reliably predict the correct properties 1 . There are several possible but computationally expensive routes one can take to tackle this problem. One of them is to map the many body Hamiltonian to a Hubbard type model and to use Dynamical Mean Field Theory (DMFT) 3 to solve it. The approach we choose here is to consider small lattice sizes, akin to nano crystals or quantum dots [4][5][6][7] , and to use precise techniques related to exact diagonalization to compute the ground state and the time evolution.
An important effect for, e.g., the solar power industry is impact ionization [8][9][10] , also called Multiple Exciton Generation (MEG) 7,11-13 , Carrier Multiplication (CM) 7 , or Multi Carrier Generation (MCG) 10 , where an excited electron with kinetic energy larger than the band or Mott gap of the material generates additional excitations through electron-electron scattering. Impact ionization has been observed experimentally [11][12][13][14][15][16][17] . It would allow to raise the theoretical efficiency limit of solar cells from approximately 30% 18 to approximately 60% 19 . However, impact ionization cannot be used to significantly increase the efficiency of classical silicon based solar cells 20 , because the electron-phonon scattering process is faster by an order of magnitude than the relevant electron-electron scattering processes 21 . In strongly interacting Mott insulators, on the other hand, the relevant electron-electron scattering processes can be much faster than electronphonon scattering processes due to strong localization and Coulomb interaction mentioned above. It has therefore been proposed to use Mott insulators for photovoltaic purposes 8,9,22 . Impact ionization has been numerically confirmed to exist for the infinite dimensional hypercubic lattice in a DMFT calculation 23 and a quantum Boltzmann approach 24 , and in DMFT calculations for a correlated layer connected to two metallic leads 25 and for models of LaVO 3 and YTiO 3 26 . The drawback of DMFT is that it neglects spatial correlations 3 . In the present paper, we take such correlations into full account. We examine the 4 × 3 Hubbard model, exposed to a short laser pulse 23,[27][28][29][30][31][32] . We use Lánczos-based methods 33 to compute the ground state, and a related scaling method 34 to compute the time evolution. An independent study by Kauch et. al. 32 observed impact ionization in Hubbard models, with the focus on different geometries, disorder, and the quasiparticle spectrum.
When the intensity of the electric field is high, one enters the domain of nonlinear optics, where multiphoton absorption processes 11,35-47 become relevant. With multiphoton absorption it is possible to reach excited states beyond a gap, even when the energy of a single photon is smaller than the gap. Multiphoton absorption is a single physical process described by Quantum Electrodynamics 48 , different from the consecutive absorption of single photons 35 . Experimentally, the necessary field strengths are achieved by focusing a pulsed laser beam onto a very small (order of 10 −9 cm 2 ) area 36 . Relevant applications include multiphoton microscopy 36,37 and high-resolution three dimensional polymerization of photoresists [38][39][40] . High resolution in comparison to single photon absorption is achieved because of the nonlinear dependence of multiphoton absorption on the field intensity 40,41 . The effect can be utilized in a variety of materials including molecules 35,42 , nano-crystals [43][44][45] , cuprates 46,47 , and chalcogenide glasses 41 . We observe multiphoton absorption at high electric field intensity with a driving frequency below the Mott gap size. An interesting phenomenon that has been observed for many quantum systems is that the long time mean of the expectation value of an observable can tend to a value which depends only on the energy of the initial state. This is the topic of the so called Eigenstate Thermalization Hypothesis (ETH) 49,50 . The dependence on energy only can be understood when the initial state is dominated by a single peak in the eigenstate spectrum of the Hamiltonian and the relevant observable varies slowly in eigenenergy. In the present paper, we observe and explain a similar dependence on energy only, for states with support in a very large energy range.
In Sec. II we provide an overview over the model, briefly present numerical methods, and discuss expectations based on a quasiparticle picture. Section III A shows the time dependence of the double occupancy during and after the photo pulse, including impact ionization. We calculate the emerging eigenstate spectrum from the Loschmidt amplitude in Sec. III B, observing a clear peaked structure with distances in multiple of the photon energy. The time evolution of the individual peaks is analyzed in Sec. III C. The long time behavior of the double occupancy is mostly governed by the amount of absorbed energy during the pulse. A connection to the Eigenstate Thermalization Hypothesis is explored in Sec. III D.

II. MODEL
We investigate a 4 × 3 Hubbard model with open boundary conditions and nearest neighbor hopping v ij (t), at half-filling, i n i,↑ = i n i,↓ = N 2 , N = 12. The initial state at time t = 0 is the ground state, calculated with the Lánczos method 33 . At time t > 0 the system is irradiated by an electric field E of frequency Ω, in plane with the 2d system, and angled at 45 • to the lattice. We incorporate the field into the Hubbard model by the Peierls phase 51,52 , ∂t . Then all nearest neighbor hopping amplitudes obtain the same time-dependent complex phase factor v ij (t) = v e iφ(t) . ( We employ a Gaussian shaped light pulse 23,29,31,32 , centered at time t i = 5v −1 , with a width of σ = 2v −1 and intensity I 0 , see Fig. 1, so that We investigate the model at U = 8v. In the following, all energies will be specified in units of the hopping amplitude v and all times in units of v −1 .
Our main observable of interest is the total double occupancy of the system: as a measure of excitations of the system. We will use an increasing double occupancy after the incoming photon pulse as a measure for impact ionization 23 . We approximate the time-evolution operator U by We use a modified scaling and squaring method to compute the action of the matrix exponential on the state 34,53 . Convergence of the simulations was verified by using several different time steps, ∆t = 10 −2 · 2 −n , n ∈ {0, 1, 2, 3, 4}. For a detailed description of a similar numerical setup see Ref. 31. The density of states was calculated with H(0) from the Fourier transform 54), averaging over all sites i.

A. Quasiparticle picture
The density of states is shown in Fig. 2. There is a gap of E gap ≈ 4.7 and the width of the Hubbard bands is E bw ≈ 5.8. The reaction of the system to the pulse will depend on the value of Ω. The excitations expected to happen in the single particle picture are illustrated in Fig. 3. When Ω < E gap , then the system should not react to the pulse, because an incoming photon does not carry enough energy to excite electrons. When E gap < Ω < E gap + 2E bw , a photon can excite an electron to the upper Hubbard band, leaving behind a hole, and for Ω > E gap +2E bw , no final state is available, so the photon should not be absorbed.
Electron scattering can later modify the number of excited electrons. When E gap < Ω < E bw + Egap 2 (third subfigure of Fig. 3), Auger recombination 8 , can reduce the number of excited electrons. Conversely, when 2 E gap < Ω < 2 E bw + E gap (and also E gap < E bw ), then an excited electron can subsequently transfer enough energy by scattering processes to excite another electron into the upper Hubbard band. Thus a single photon can produce two excited electrons in this process of impact ionization 8,23,25,26,32 . Note that this is an idealized view.
In nonequilibrium, the spectral function is time dependent 24,28,29,31,32,42,[55][56][57][58] and there is a photo-induced insulator-metal transition in the Hubbard model in various setups 27,29,31,32,56,57,[59][60][61] . Such an insulator-metal transition has also been observed in experiment in quasi one-and two-dimensional materials [62][63][64][65][66][67] . We also note that the incoming photon pulse has a finite widthσ = 1 Expected excitation processes in the quasiparticle picture, for different ranges of the photon energy Ω. In each subfigure, the horizontal axis represents energy and the semicircles schematically depict the lower and upper Hubbard bands. An electron is excited by a photon (step 1) into the upper Hubbard band, given that the photon energy is not too low (first subfigure) or too high (third subfigure). For the energy range depicted in the fourth subfigure, Auger recombination can reduce the number of excited electrons (step 2). At higher photon energies (last subfigure), a single photon can lead to two excited electrons via impact ionization (step 3). in frequency, due to the finite width of the Gaussian peak in the time domain (see Fig. 1). ble occupancy oscillates slightly during the pulse but then returns to almost exactly its original value. At frequencies Ω between E gap and E gap + 2E bw , where absorption is expected, the double occupancy rises quickly during the pulse, i.e., elecrons are excited across the gap and energy is transferred into the system. Note that after the pulse has decayed, the Hamiltonian H(t) becomes the time-independent bare H(0) again, so that energy is then conserved. At Ω = 9.1, just below the region where impact ionization is expected, the double occupancy stays (almost) constant after the pulse. At Ω = 10.7, the double occupancy shows the expected impact ionization behavior, noticeably rising further after the pulse, which was also observed by Kauch et al. 32 . We note that the corresponding time scale is large compared to the electron hopping time 23,26,32 At Ω = 7.5, the double occupancy goes down after the pulse, compatible with the expected Auger recombination, beginning already after the pulse maximum at time t = 5.
Remarkably, at Ω = 3.2, below E gap , where the quasiparticle picture would forbid excitations, Fig. 4 exhibits a sizeable increase of the double occupancy. We will later show that this is due to multiphoton excitations.
In Fig. 5 we display the double occupancy at t = 20 after the pulse, for different intensities I 0 and different frequencies Ω. There is a strong nonlinear dependence on I 0 . We will show in Sec. III B that it occurs together with multiple photon excitations. Indeed, for Ω below E gap , a sizeable excitation of double occupations does not occur for small intensities, but only at large I 0 , as would be expected for multiphoton absorption.
The energy absorbed by the Hubbard system matches the change in double occupancy in Fig. 5 closely, which will be further explored in Sec. III D.
In Fig. 6 we show the further change of double occupancy after the pulse, from t = 20 to t = 300 where it has converged well for all Ω. In the energy range 2 E gap < Ω < E gap + 2 E bw the double occupancy increases, i.e., there is impact ionization, as suggested by the quasiparticle picture. Impact ionization is larger at higher intensities. We will show below that this is again connected to the absorption of several photons during the pulse. Notably, when the intensity is large, impact ionization even occurs at small Ω < E gap . The double occupancy decreases after the pulse in a range of lower values of Ω. For the small I 0 , this range closely matches the expectation E gap < Ω < E bw + Egap 2 from the single particle picture, while at I 0 = 0.71 and I 0 = 0.97, the range extends to larger energies. We will discuss these processes further in Sec. III C.

B. Eigenstate spectrum
For a better understanding of the excited state after the pulse, we examine the eigenstate spectrum of |ψ(t) . To this end we compute the Fourier transform of the Loschmidt amplitude, with respect to an auxiliary time span τ (Refs. [68][69][70][71]. Equation (6b) can also be viewed as the probability distribution of work done on the system 72 by the electric field. In Eq. (6b), the sum runs over all eigenstates |n of H(0), and E n is the energy of |n . Through Eq. (6b) we obtain the energy spectrum of an arbitrary state, while only using the tool of time-evolution. This is very useful here, because the large Hilbert space dimension prohibits full diagonalization of the Hamiltonian. We note that the eigenstate spectrum is independent of time after the pulse has decayed, because then the Hamiltonian H(t) reverts back to the initial Hamiltonian H(0). We dampened the time evolution in Eq. (6a) with exp −t 2 /(2σ 2 ) , σ = 5/ √ 2, thus widening the spectra byσ = 1/σ ≈ 0.3.
In Fig. 7 we show the resulting spectrum for the case of Ω = 10.7, after the pulse has decayed. The spectrum has a very distinct peaked structure, with the distances between the peaks close to Ω. The excited state |ψ(t) thus consists of groups of eigenstates of H(0) close to multiples of the photon energy. At low intensity, I 0 = 0.26, the peaks are narrow, and excitations at 1Ω and 2Ω dominate. For large intensity I 0 = 0.71, the peaks are wider, shifted slightly to higher energies, and they include higher multiples of Ω. The ground state is then almost depleted, in line with the saturation of double occupancy in Fig. 5. Since 2Ω is larger than E gap + 2E bw here, peaks beyond 1Ω must correspond to the sequential excitation of several electrons from the lower to the upper Hubbard band, by several photons. We will investigate these peaks individually in Sec. III C. Below the gap, Ω < E gap , Fig. 5 showed that the double occupancy is increased by the light pulse when the intensity I 0 is large, even though single photon absorption is energetically forbidden. In Fig. 8 we show corresponding eigenstate spectra for Ω = 3 at different intensities I 0 . At small I 0 = 0.26, where there is very little absorption, Fig. 8 shows that |ψ(t) has almost returned to the ground state (slightly widened in the plot), with an additional small excitation at 2Ω. At larger intensity I 0 = 0.71, however, states with higher multiples of the photon energy have become excited. Thus the increase in double occupancy here is indeed due to multiphoton absorption. The largest amplitude in Fig. 8 is from threephoton excitations at 3Ω = 9, where the available phase space is the largest as indicated by the highest absorption values in Figs. 4 and 5. At excitation energy 4Ω, both 4-photon absorption and two sequential two-photon absorptions can contribute. At 5Ω = 15 and beyond, the end of the bandwidth has been reached, so that the excitations, still sizeable in Fig. 8, must correspond to sequential absorptions. We note that since the density of states develops a small in-gap density after electron excitations 28,29,31,32,56 , additional processes with singlephoton absorption will also be possible after the initial two-photon absorption has taken place.

C. Time evolution of individual photon-absorption peaks
Knowing the structure of the eigenstate spectrum, one can gain considerable additional insight into the development of |ψ(t) . Here we introduce a technique to decompose a state |ψ(t) into states |f i (t) with support around the individual photon peaks. We isolate the contribution of a photon peak to |ψ(t) by applying a filter, We take N (t) such that ϕ(t)|ϕ(t) = 1. This is a Gaussian peak centered around E f . Expanding |ϕ(t) in the eigenbasis of H(0) shows that basis states which are energetically too far away from E f are discarded in |ϕ(t) . Note that the filtering from |ψ(t) to |ϕ(t) is very similar to doing a time-evolution.
Suppose that we have derived a number of distinct states |f i (t) from a state |ψ(t) with the filtering procedure above. To find the best representation |ψ (t) of |ψ(t) that is a linear combination of filter results |f i (t) we use: (8) and choose the coefficients α i (t) such that |ψ(t) − |ψ (t) 2 is minimized. We will look at the expectation value of an operator O with respect to the single states |f i (t) , including the coefficients α i (t): (with the time dependence not denoted). Specifically for the photon peaks we define |f i (t) as follows: In other words, |f i (t) is a filtered state centered around the energy that is i times Ω above the ground state energy, filtered with a width of σ f = Ω 3 . Note that each filtered state still contains many eigenstates, so that observables, like the double occupancy, remain time dependent even after the light pulse.
We find that the set of filtered states |f i (t) provides a good representation of the original state |ψ(t) at all times. The absolute value of the overlap between the filtered states | f i (t)|f j (t) | is about 5·10 −2 when they are Ω apart, 10 −4 when 2Ω apart, and 10 −8 when they are 3Ω apart. The norm squared of the difference between the original time-evolved state and the best approximation |ψ(t) − i α i (t) |f i (t) 2 is of the order of 10 −2 . The relative difference between the expectation value of the double occupancy of the original states and the sum of filtered state contributions Eqs. (9a) and (9b) is of the order of 0.5 %. We examine two cases, Ω = 10.7 where the double occupancy increases, and Ω = 7.5 where it decreases after the pulse. In both cases we find that during the pulse, the initial rise of the photon peaks occurs sequentially, delayed by roughly one hopping time for each additional photon, in line with the expected sequential nature of photon absorptions at these values of Ω.
The eigenstate spectra of the two states after the pulse are shown in Figs. 9 and 10 (center row). At every time t, we separately filter out each peak and calculate the double occupancy according to Eq. (9a), including the coefficient |α i (t)| 2 . The results are shown in Fig. 9 and Fig. 10 (top and bottom rows).
For Ω = 10.7, we saw in Fig. 4 that there is an increase of double occupancy after the pulse, in line with the impact ionization expected in the quasiparticle picture. Figure 9 shows that indeed each photon peak and also the non-diagonal part contribute to the increase. The largest contributions to the double occupancy and to its increase come from the states with two and three absorbed photons, each of which can separately contribute to the impact ionization process sketched in Fig. 3.
At Ω = 7.5, absorption of up to 6 photons is visible, and the remaining contribution of the ground state is tiny. Now |ψ(t) has a decreasing double occupancy as a function of time. The individual peaks in Fig. 10 all contribute to this decrease, with the notable exception of the single photon peak, for which the double occupancy stays almost constant.
This difference in behavior is in line with Auger recombination in the quasiparticle picture, shown in Fig. 3 (fourth subfigure), which is only possible when at least two electrons are excited into the upper Hubbard band. Thus this decay channel is absent in the single photon peak. The upper bound for this process is Ω < E bw + ing Auger-like processes with more initial photons and larger energy range. Furthermore, for more initially excited electrons there are more decay channels, suggesting a stronger decrease of double occupancy, in agreement with the behavior of the large 3Ω peak in Fig. 10, which shows the largest change in double occupancy after the pulse.

D. Eigenstate Thermalization
Another way to learn about the photoexcited system is to look at the double occupancy as a function of ab-sorbed energy, shown in Fig. 11. Interestingly, whereas the double occupancy at time t = 20 depends on intensity and frequency separately, after convergence at long simulation times it depends almost only on the absorbed energy, with an almost linear relation between those two quantities.
In many physical systems it has been observed that long time averages of some expectation values are actually close to a microcanonical average, which has been discussed under the name of Eigenstate Thermalization Hypothesis (ETH) 49,50,73 .
When the Hamiltonian is not time-dependent (long after the light pulse in our case), a quantum state under Under the assumption of no degeneracy, only the diagonal contributions survive for long times: The Eigenstate Thermalization Hypothesis 49,50,73 tries to explain the peculiar observation that this long-time mean of O(t) often coincides with a microcanonical average, where we take the average of the expectation values of the eigenstates in a small energy window ∆E around E ψ = n E n | n|ψ | 2 , For t t 0 : This is rather surprising, because the left hand side of the equation above depends on the specific coefficients of the state in the eigenbasis | n|ψ | 2 , while the right hand side does not. One possible explanation is that for many physical situations the state has a single peak in the eigenstate spectrum and in addition, the difference between two expectation values of two separate eigenstates is small when the difference between their eigenenergies is small: In our case, however, we do not have a single peak in the eigenstate spectrum, but a series of peaks that stretch over almost the whole eigenenergy range, as can be seen in Figs. 9 and 10. We now show that eigenstate thermalization still holds when there is an almost linear relation between the eigenenergies and the expectation values of their respective eigenvalues (Fig. 12, inset) in the energy range(s) where |ψ(t 0 ) has an overlap with the eigenstates: with M the maximum deviation from the linear behavior. Then the long-time average of O(t) has the same linear behavior as a function of E ψ with at most the same maximum deviation M , no matter how peaked the structure of |ψ(t 0 ) is in the eigenstate spectrum. Namely, for large t t 0 , where Eq. (13) holds, we have Similarly, for the difference between the long time average and the microcanonical average Eq. (14) with ∆E and N from Eq. (14) and M the maximum deviation from linear behaviour within the small energy range ∆E. The upper bounds apply when |ψ has support only from eigenstates where the deviation Eq. (16b) is maximal.
We examined the linearity Eq. (16a) of the double occupancy in the Hubbard model by uniformly sampling states from the hypersphere of normed states 74 . For each value of E f in Fig. 12, we took a sample of N r = 9 states |r , projected and normalized each state to a small range around E f by |r f = 1 N e − (H−E f ) 2 2 |r , and plotted the averaged double occupancy d(E f ) = 1 Nr r f r f | i n i↑ n i ↓|r f , which is similar to the microcanonical average Eq. (14), since contributions that are nondiagonal in eigenstates cancel stochastically. The approximately linear behavior in Fig. 12 was verified for a smaller 3 × 2 system where exact diagonalization is still possible. Figure 12 indicates that for the Hubbard model, the double occupancy of eigenstates is indeed close to linear in the eigenstate energies, like Eq. (16a), with the same slope as in Fig. 11. The steps in the figure correspond to individual double occupations.
In our case, the double occupancy after the light pulse converges (up to small fluctuations) at large times t. The converged value, shown in Fig. 11, is then the same as the long time average. Taking into account the actual eigenstate spectra of the excited states, the convergence towards linear behavior in Fig. 11 indeed corresponds to Eq. (17).

IV. SUMMARY
We investigated the non-equilibrium response of a strongly correlated Mott insulator to a short light pulse, using exact-diagonalization based calculations on a 4 × 3 Hubbard model for a large range of light intensities and of photon energies Ω. The pulse excites electrons into the upper Hubbard band, quickly increasing the number of doubly occupied sites. At sufficiently large photon energies, we observed impact ionization (also seen in Ref. 32), namely a further increase in the double occupancy over time after the light pulse had ended. Conversely, at lower photon energies, we observed Auger recombination with a reduction in double occupancy.
We calculated the eigenstate spectra of the nonequilibrium states, i.e. the probability distribution of the work done by the light pulse, as the Fourier transform of the Loschmidt amplitude. The resulting spectra exhibit distinct peaks at distances of about Ω, corresponding to the absorption of multiple photons. The absorption rate is strongly nonlinear in light intensity. Multiphoton absorption was identified for small photon energies below the band gap, which leads to electron excitations when the light intensity is large.
We introduced a technique, using tools similar to time evolution, to isolate photon peaks in the eigenstate spectra. This enabled us to investigate the nonequilibrium evolution of double occupancy in individual photon peaks. We showed for example that, as expected from the quasiparticle picture of Auger recombination, double occupations excited by a single photon do not contribute to the overall reduction over time which occurs at intermediate light frequencies.
We found that at large times, the double occupancy moves towards a function that only depends on the absorbed energy, reminiscent of the Eigenstate Thermalization Hypothesis. Eigenstate Thermalization is usually observed for states with support in a narrow region of eigenenergies. We showed that the observed dependence on energy alone will also happen for states which contain a wide range of eigenenergies, when the relevant observ-able is almost a linear function of energy.
The analysis of eigenstate spectra via the Loschmidt amplitude and the filtering of the relevant energy ranges provided valuable insight into multiple photon absorptions and should prove to be useful tools to investigate strongly correlated systems, when full diagonalization is not possible, but the computation of the time-evolution is accessible.