Photoinduced $\eta$-pairing in One-dimensional Mott Insulators

Employing the density-matrix renormalization group technique in the matrix-product-state representation, we investigate the photoexcited superconducting correlations induced by the $\eta$-pairing mechanism in the half-filled Hubbard chain. We estimate the characteristic pair correlation function and verify the accuracy of our numerical results by comparison with exact-diagonalization data for small systems. The optimal parameter set of the pump that most enhances the $\eta$-pair correlations, is calculated in the strong-coupling regime. For such a pump, we explore the possibility of quasi-long-range order.


Introduction
Superconductivity is one of the most fascinating phenomena in condensed-matter physics and has attracted continuous attention since its discovery in 1911 [1,2]. While the basic features of the superconducting (SC) phase are well established for equilibrium systems, the photoinduced superconductivity is subject of intense research because of the recent pumpprobe experiments, e.g., in copper oxides [3][4][5] or in K 3 C 60 [6]. In these materials, the transient optical spectra imply the SC-like properties even at elevated temperatures.
Very recently, it was theoretically demonstrated that unconventional SC correlations can be induced by pulse irradiation in a simple Mott insulator described by the half-filled Hubbard model [7]. This SC state stems from the so-called η-pairing mechanism [8][9][10], characterized by staggered pair-density-wave oscillations in the off-diagonal long-range correlations [11]. However, since these results were obtained by the exact diagonalization (ED) technique, the accessible system size was only L 16 with periodic boundary conditions (PBC) [7]. It is, thus, vitally important to discuss the finite-size effects in the η-pairing state and pair correlations especially at larger distances.
To make progress in this direction, we focus on one-dimensional systems, where the unbiased density-matrix renormalization group (DMRG) technique [13] can be applied. (We note that η-pairing state in the photoexcited state can also be found in two dimensions [7], but is hard to calculate in an approximation-free way). Simulating the real-space pair correlation function and its structure factor, we explore the conditions under which η-pairing is most likely. For large enough system sizes, we demonstrate a peculiar increase of the pair structure factor with the parameter set of maximally enhanced η-pairing state after the photo excitation, which might be the signature of the (quasi-)long-range order of pairs.

Theoretical approach
Model. Hereinafter, we study the photoinduced η-pairing state in a half-filled Hubbard chain, defined byĤ whereĉ † j,σ (ĉ j,σ ) creates (annihilates) an electron with spin projection σ =↑, ↓ at Wannier site j, andn j,σ =ĉ † j,σĉj,σ . The transfer amplitude t h enables the electrons to hop between neighboring lattice sites, whereas the on-site Coulomb repulsion U tends to localize the electrons, establishing a Mott insulating ground state.
Introducing a time-dependent external field into the hopping term by the Peierls substitution [14], with amplitude A 0 , frequency ω p and pulse width σ p centered at time t 0 (> 0), makes the Hamiltonian time-dependent:Ĥ →Ĥ(t). As a consequence the equilibrium ground state |ψ(0) at t = 0 evolves in time |ψ(t) . To take account of this time evolution, in our numerics, we apply the time-evolving block decimation method [15] in combination with the secondorder Suzuki-Trotter decomposition. In the following, we use t h (t −1 h ) as the unit of energy (time), and the time step δt is set to be δt · t h = 0.01.
Pairing correlations. The photoinduced η-pairing state can be characterized by the real-space pair correlation function with the on-site singlet pair operator∆ j =ĉ j,↑ĉj,↓ . Here, N b = L − r denotes the number of pairs of sites separated by distance r in a system of L sites with open boundary conditions (OBC). For r = 0 the pair correlation is equal to twice the double occupancy, i.e., P (r = as demonstrated in Ref. [7]. Note that it is important to analyze the Fourier transform of P (r, t), i.e., P (q, t) = r e iqr P (r, t), which shows a characteristic enhancement after the pulse irradiation in the half-filled Hubbard model.
Accuracy check. Before we show our main results, we like to discuss the accuracy of the time-dependent DMRG simulations with OBC. During a DMRG simulation one chooses the maximum of the so-called bond dimension χ and makes sure that the truncation error in the singular value decomposition stays sufficiently small. The other way around, the maximum truncation error can also be fixed in the simulation. Figure 1 demonstrates the latter case, simulating the pair structure factor P (q = π, t) by DMRG. The DMRG data are compared with ED results in Fig. 1(a) for L = 12. Keeping the cutoff smaller than 1 × 10 −7 , P (q = π, t) obtained by DMRG is in perfect agreement with the ED results, even up to time t · t h = 80. The deviation only becomes significant when a larger cutoff is used for t t 0 . In each case, the maximum bond dimension χ max increases rapidly around t ∼ t 0 and afterwards stays constant for t · t h 25 [see Fig. 1(b)]. With increasing system size, the memory required by ED calculations increases exponentially: It is already problematic to carry out a timedependent ED simulation for L = 16. Clearly the computational costs of the time-dependent DMRG calculations also increases with the system size L, but the truncation error can be kept below 1×10 −8 in the case of L = 16. To keep the cutoff less than 1×10 −7 up to t·t h = 30, χ max should be about 10000 as in Fig. 1(d). As shown by Fig. 1(c), the difference between the DMRG results for the maximum truncation errors 1 × 10 −7 and 1 × 10 −8 is negligible. Thus, it is sufficient to keep the truncation error less than 1 × 10 −7 in the time-dependent DMRG simulations. This is computationally expensive, however, due to the rapidly increasing bond dimensions. Alternatively, by keeping bond dimensions χ = 800, reasonable agreement can also be obtained up to t · t h ≃ 20 [see the dotted line in Fig. 1(c)].

DMRG results
We now determine the optimal parameter set for the enhancement of P (q = π, t) and examine the behavior of the pair correlations for these parameters. Figure 2 gives the OBC DMRG contour plots of P (π, t) at various times for different values of A 0 and ω p . For t < t 0 , the magnitude of the pair correlation functions is still marginal. When t t 0 , a noticeable enhancement of pair correlations appears in wide parameter regions. Furthermore, the spectral intensity of P (π, t) is concentrated in a single spot after pulse irradiation (t · t h 20). Interestingly, the height of the peak at ω p /t h ≃ 6.8 and A 0 ≃ 0.4 increases at long times t ·t h = 30. Figure 2 also demonstrates the system-size dependence of the pair structure factor. For L = 8 (upper panels of Fig. 2) the stripe structure can be observed after the pulse irradiation as in the case of ED results with L = 14 [7]. As explained in Ref. [7], the peak structure of P (π, t) is essentially the same as the ground-state optical spectrum χ JJ (ω) with some Lorentzian  Contour plot of P (q = π, t) at t · t h = 8, 10, 20, and 30 for U/t h = 8, σ p = 2 and t 0 · t h = 10 in the ω p -A 0 plane. All data obtained by DMRG with bond dimensions χ = 800.
broadening η L , depending on 1/σ p . If the system size is too small for some (small) η L , the spectrum is described by a series of peaks. The stripe structure found in P (π, t) for L = 8 reflects this finite-size effect. These stripes merge and construct a single peak structure for larger system sizes L ≥ 16. With increasing L the enhanced regime is narrowed. Note that the value of ω p ≃ 6.8 at the peak position corresponds to about the size of the Mott gap ∆ c , i.e., the peak position of χ JJ (ω) can be estimated to be ω ∼ 1.461∆ c ∼ 6.84 [16]. This is in accord with the DMRG data for L = 24, although Eq. (8) of Ref. [16] is only valid in the weak-coupling regime.
Let us analyze the pair correlations at the peak position of the contour plot for L = 24 and t · t h = 30, i.e., A 0 = 0.36 and ω p /t h = 6.8. Obviously, after the pulse irradiation, P (π, t) keeps increasing gradually for t · t h 20, as shown in Fig. 3(a), while P (π, t) by ED with PBC saturates at some constant value [7]. Note that the results do not depend strongly on bond dimensions for large times. Figure 3(b) displays the real-space pair correlation function. P (r, t) for t · t h = 30 is clearly enhanced for large distances r L/2 comparing with those for t · t h = 20.
Of course, L = 24 is still too small to examine the behavior of correlation functions and boundary effects are showing up in P (r, t) for large distances. Moreover, the definition of the pair correlation with OBC, Eq. (3), is not equivalent to the usual ones with PBC in Ref. [7]. These data, however, might imply the possibility of a (quasi-)long-range order of η-pairs, i.e., the power-law decay of pair correlations for large times. Further studies for larger system sizes are desirable.

Summary
In conclusion, we demonstrated η-pairing in the one-dimensional Hubbard model at half filling by means of unbiased density-matrix renormalization group simulations. For small system sizes, the time evolution of the corresponding pair correlation functions can be computed with high accuracy-i.e., the maximum truncation error is always less than 1 × 10 −8 -up to times t · t h = 80 (30) for L = 12 (16) with open boundary conditions. Although the numerical accuracy will get worse as L increases, the resonant parameter set can be determined up to L = 24. For these pump parameters the pair structure factor P (q = π, t) is enhanced and magnitude of the pair correlation functions increases for long distances with time after pulse irradiation. This can be taken as strong indication for η-pairing and off-diagonal (quasi) longrange order. Since boundary effects still show up in the pair correlation simulation data, it is highly desirable to prove quasi-long-range order for larger system sizes or, even better, directly in the thermodynamic limit by using, e.g., the infinite time-evolving block decimation approach [17].