Photoinduced $\eta$-pairing at finite temperatures

We numerically prove photoinduced $\eta$-pairing in a half-filled fermionic Hubbard chain at both zero and finite temperature. The result, obtained by combining the matrix-product-state based infinite time-evolving block decimation technique and the purification method, applies to the thermodynamic limit. Exciting the Mott insulator by a laser electric field docked on via the Peierls phase, we track the time-evolution of the correlated many-body system and determine the optimal parameter set for which the nonlocal part of the $\eta$-pair correlation function becomes dominant during the laser pump at zero and low temperatures. These correlations vanish at higher temperatures and long times after pulse irradiation. In the high laser frequency strong Coulomb coupling regime we observe a remnant enhancement of the Brillouin-zone boundary pair-correlation function also at high temperatures, if the Hubbard interaction is about a multiple of the laser frequency, which can be attributed to an enhanced double occupancy in the virtual Floquet state.

Introduction.-Optical pumping is not only an excellent tool to investigate complex few-and many-body systems but also makes it possible to create new phases of quantum matter with tunable properties [1][2][3][4]. Inducing superconductivity by light pulses in low-dimensional materials with strong electronic correlations is certainly one of the most fascinating options in this regard [5][6][7]. Thus, it was not surprising that a whole series of theoretical studies has addressed the microscopic modeling and understanding of this nonequilibrium light-matterinteraction phenomenon [8][9][10][11][12].
In this context, the so-called η-pairing, originally proposed by Yang for the Hubbard model [13], has attracted renewed attention [14][15][16][17][18][19][20][21]. Pumping the Mott insulating phase may results in an excited state with enhanced offdiagonal pair-density-wave correlations, which are absent in the ground state [15]. Here the basic mechanism is the creation of η-pairs triggered by the nonlinear optical excitation of the system in conjunction with the selection rules. Interestingly, for low-amplitude pulses, the peak structure of the pair correlation function is essentially the same as that obtained for the optical spectrum in the ground state, implying that the photoinduced state might indeed result from an η-pairing mechanism.
The crucial question is whether these findings will remain valid in the thermodynamic limit and at finite temperature T . Some features, e.g., the stripe structure found in the structure factor of the pair correlations ( Fig. 2 of Ref. [15]) by exact diagonalization (ED) of small systems, have been shown to disappear by increasing the system size [22], exploiting density-matrix renormalization group (DMRG) and time-evolving block decimation (TEBD) methods [23,24]. For sure, determining the temporal evolution of an infinite, driven, strongly correlated electron system at T > 0, is one of the most difficult problems in solid state theory. Since the fermionic Hubbard model [25] can nowadays be realized in optical lattices [26][27][28][29], just as its bosonic counterpart [30], such a theoretical treatment is indispensable, however, for the interpretation of the experimental data, especially in one spatial dimension.
Despite this difficulty, this Rapid Communication aims at proving the existence of photoinduced η-pairing in the one-dimensional half-filled fermionic Hubbard model, directly in the thermodynamic limit and for finite temperatures. For this we exploit unbiased numerical techniques, specifically the infinite TEBD (iTEBD) technique [31] based on an infinite matrix-product-state (iMPS) representation [32] in combination with the purification method [33,34], which enables us to monitor the realtime evolution of thermal states at finite target temperature, as accessible by optical-lattice experiments.
Model.-Our starting point is the Hubbard Hamiltonian,Ĥ whereĉ † j,σ (ĉ j,σ ) creates (annihilates) a fermion with spin projection σ (=↑, ↓) at lattice site j, andn j,σ = c † j,σĉj,σ . The first term represents the kinetic energy (with nearest-neighbor particle hopping amplitude t h ) that acts against the Coulomb interaction (parametrized by the on-site Hubbard repulsion U ), which tends to lo- j,↑ denotes the on-site singlet pair creation operator, see [35].
As demonstrated in Ref. [15], photoinduced η-pairing states may appear when an external time-dependent field couples to the hopping term via a Peierls phase [36], where the vector potential describes a pump pulse with amplitude A 0 , frequency ω p and width σ p , centered at time t 0 (> 0). As a result the Hamiltonian becomes time-dependent,Ĥ →Ĥ(t), and the initial (equilibrium) ground state evolves (forward) in time: |ψ(0) → |ψ(t) . Numerically such a time evolution can be treated in an efficient manner by combining TEBD and second-order Suzuki-Trotter decomposition methods [24]. Hereafter we use t h (t −1 h ) as the unit of energy (time), and set the time step δt · t h = 0.01.
In fact, using the iTEBD technique, we directly examine the time evolution of the pair correlation function, in case that the number of lattice sites L → ∞. At r = 0, the pair correlation gives twice the number of double occupancy, i.e., P (0, t) = 2n d (t) = (2/L) j ψ(t)|n j,↑nj,↓ |ψ(t) . Most notably, the Fourier transform P (q, t) = r e iqr P (r, t) shows an enhancement after the pulse irradiation that was believed to be indicative of η-pairing in finite Hubbard clusters [15]. Since we are particularly interested in longer-range pair correlations, we will also analyze the modified structure factor P r>0 (q, t) = r>0 e iqr P (r, t), in which the contribution of the double occupancy n d (t) is excluded. Let us point out that P (q, t) obtained by iTEBD in iMPS representation fulfils the relation P (π, t) = 2 ψ(t)|η +η− |ψ(t) /L, which is not the case in any (finite-system) TEBD calculation with open boundary conditions (OBC), see Ref. [22] and the Supplemental Material [35].
ITEBD results at T = 0.-In a first step, we determine the optimal parameter set in view of an enhancement of P (π, t) at zero temperature. Figures 1(a)-(c) provide iTEBD contour plots for P (π, t), in dependence on A 0 and ω p , at different times t · t h . For t < t 0 , in the ramp-up regime of the pump field, the spectral intensity of P (π, t) is negligibly small (not shown). Noticeable pair correlations develop for t t 0 , albeit the signal is very broad [cf., Fig. 1(a)]. It becomes focused when the light   Figure 1(d) relates the time evolution of P (π, t), P r>0 (π, t) and 2n d (t) for the optimal parameter set marked by a cross in Fig. 1(c). All quantities show a clear response to pulse irradiation and will be strengthened as the system progresses in time until saturation is reached. Apparently, here, the nonlocal contributions P r>0 (π, t) have a stronger impact on P (π, t) than double occupancy.
A notable finding of previous ED calculations [15] was a peak structure of P (π, t) as a function of ω p which is essentially the same-for small A 0 -as those of the groundstate optical spectrum χ JJ (ω), folded with an appropriate Lorentzian of width η L (depending on 1/σ p ). The current-current spectral function χ JJ (ω) is given by where |ψ 0 is the ground state having energy E 0 , and the charge current operatorĴ, for the Hubbard model, takes the formĴ = it h σ (ĉ † ,σĉ +1,σ −ĉ † +1,σĉ ,σ ). The ED [15] and TEBD [22] calculations, which could be conducted for small lattices only, suffer from finite-size effects however. These give rise, inter alia, to stripe patterns in P (π, t), which makes it difficult to determine its maximum value. We demonstrate that a single peak structure evolves in the thermodynamic limit L → ∞ [see Fig. 1(c)], and therefore can address more seriously the question whether the χ JJ (ω) lineshape obtained by timedependent iMPS-based DMRG really agrees with that of P (π, t) for small A 0 and large t, where P (π) becomes time-independent. Figure 1(e) compares the iTEBD data, obtained for P (π, t) at various small A 0 and t · t h = 16, with the DMRG results for χ JJ (ω) (using η L /t h = 0.2), in dependence on ω p respectively ω. Here we show that P (π, t) divided by A 2 0 scales to the imaginary part of the optical spectrum Imχ(ω) [ P (π, t)/CA 2 0 with C ∼ 7.9] since the double occupancy n d is proportional to A 2 0 , for very small A 2 0 , in a wide ω p -range around the resonant frequency ω p U [37]. Close to the maximum in P respectively Imχ, at about ω 6.49, both quantities differ for larger amplitudes A 0 , because the nonlocal correlations contained in P r>0 (π, t) become increasingly important. Taking the relation Imχ JJ (ω) = ωσ 1 (ω) into account, where σ 1 (ω) is the real part of the optical conductivity, this behavior is in accordance with DMRG and fieldtheory results for the optical response in the half-filled Hubbard model [38].
ITEBD results for T > 0.-In a second step, we will investigate-under usage of the iMPS and purification approaches [33,34]-what happened to the η-pairing correlations at finite temperatures T = 1/β. Methodically, to obtain the equilibrium state |ψ T at some target temperature T , we first construct an iMPS representation of a state |ψ ∞ at infinite temperature, where each physical site is in a maximally entangled state with an auxil-   iary site, and then carry out the imaginary-time evolution e −βĤ/2 |ψ ∞ of the physical system. We note that combining the Suzuki-Trotter decomposition with swap gates [39], such a time evolution can be effectively implemented for any nearest-neighbor Hamiltonian. We start by checking the temperature dependence of the double occupancy n d in the pure Hubbard model (1) without optical pump Our iTEBD data in Fig. 2(a) reveal the well-known minimum in n d [40], which is shifted to higher temperatures as U increases and is related to the maximum in the local magnetic moment L 0 = 3 4 (n j,↑ − n j,↓ ) 2 [= 3 4 (1 − 2n d ) at half filling]. At T = 0, L 0 interpolates between the atomic limit (U = ∞) with L 0 = 3/4 since n d = 0 and the band limit (U = 0) where L 0 = 3/8, i.e., n d = 1/4, which is also the value for T → ∞ since empty, spin-up/down and double occupied sites are equally likely. Figure 2(b) shows the temperature dependence of P (π), together with those of P r>0 (π) and n d . At T = 0, on-site [n d ] and nonlocal [ P r>0 (π)] contributions cancel each other, so that P (π) = 0. Clearly the pairing correlations vanish in the opposite T → ∞, expressed by the fact that P r>0 (π) → 0 and the P (π)-curve tends to 2n d , see also Ref. [20]. As a result, strong η-pair correlations can be expected in the low-temperature region at best. Now, we take into consideration a time-dependent external field and carry out the real-time evolution ofĤ(t) to a thermal equilibrium state |ψ T . This allows us to discuss the development of η-pairing correlations as a function of time at T > 0. Figures 3(a)-(c) provide iTEBD contour plots of P (q = π, t), P (q = π, t) and 2n d (t) in the ω p -A 0 plane for T /t h = 1.0, at t · t h = 20, following the pulse exposure. We find a persistent enhancement of P (π, t). The crucial question is whether this enhancement can be related to the nonlocal part of the pairing correlation function, or simply stems from the on-site (double occupancy) contribution to P (r, t). The answer can be read off from the contour plot of P r>0 (π, t) [ Fig. 3(b)], which demonstrates its noticeable contribution. Figure 3(c) gives the corresponding values of double occupancy 2n d (t). Here we find two maxima at about ω p ∼ U and 2ω p ∼ U which can be assigned to resonant driving, i.e., to the existence of a Floquet virtual state [41]. How P r>0 (π, t) and 2n d (t) will influence P (q = π, t) over time can be seen in more detail in Fig. 3(d) for ω p /t h = 6.6 and A 0 = 0.5 [×position in Fig. 3(c)]. Apparently, all these quantities are growing when the light pulse acts on the correlated system [around t 0 · t h (= 10)]. Here the (photoinduced) nonequilibrium physics emerges. Note that the lineshape of P (π, t) (and especially its decay at larger times) is largely determined by P r>0 (π, t). At t · t h 20 saturation is reached. The comparison with the pure Hubbard model results shows the predicted dynamical generation of double occupancy [42,43] after pulse irradiation.
Finally, we look at the system response to the pulse at higher temperatures (T ∼ U ). Figures 4(a) and (b) display the contour plots of P (π, t), after pumping (t · t h = 20), for U/t h = 8 and U/t h = 10, respectively. Again we observe pronounced maxima when the pulse frequency is close to ω p U/m, which comes to light for m = 1 , 2 in (a) and m = 1 , 2 , 3 in (b). Figure 4(c) elucidates the origin of this multi-peak structure and the significant differences to the behavior at low-temperatures shown in Fig. 3(d). Before pulse irradiation and at long times (where P (π, t) reaches its saturation value), P (π, t) is completely determined by 2n d (t). The pure Hubbard model result is maintained up to t · t h 7.5 (cf. the dotted line), which can be considered as the linear response regime [41]. The nonequilibrium dynamics is evidenced at intermediate times 7 t·t h 13, when the irradiation is strong. In contrast to low temperatures, the contribution of P r>0 (π, t) is negligible after pulse irradiation for t·t h 18. This shows that the peak structure observed in Figs. 4(a) and (b) can be attributed to the enhanced double occupancy. The high-frequency expansion in the Floquet picture reveals the underlying mechanism: Performing a Schrieffer-Wolff transformation [44,45] for a periodically driven Hubbard model in the strong-coupling regime, an effective (Heisenberg) Hamiltonian is obtained 4 6 8 10 12 0  (see, e.g., Ref. [46]), containing an effective exchange interaction J eff , which diverges at the resonant frequencies U mω p [47]. Since time periodicityĤ(t + τ ) =Ĥ(t) with τ = 2π/ω p is absent in our model (1) and (2), the photoinduced double occupancy appears as a Floquet virtual state as in the nonequilibrium dynamics of pumped Mott insulators [41]. This effect can be observed at any temperature, see, e.g., Fig. 3(c).
Conclusions.-To sum up, we have demonstrated light-pulse photoinduced η-paring in the one-dimensional half-filled Hubbard model at both zero and finite temperatures by means of a defacto approximation-free numerical approach. For zero temperature, we carved out finite-size effects of previous exact diagonalization studies, but confirmed the basic relation between the pair correlation function and the ground-state optical spectrum for the infinite system. With a view to experiments, also the optimal pulse for an enforcement of η-pair correlations P (π, t) is determined. For finite but low temperatures, nonlocal pairing correlations P r>0 (π, t) were detected within the applied iTEBD-purification scheme. After pulse irradiation a dynamically generation of double occupancy is proved for finite temperatures. Overall, our results support a scenario where optical excita-tion of a Mott insulator may lead to a (nonequilibrium) state with very slowly decaying pairing correlations. If fermionic optical lattices will be cooled to temperatures T J ex = 4t 2 h /U < t h in the strong-coupling regime (U t h ) [29], our findings should be detected in the laboratory.

η-pairing symmetry of the Hubbard Hamiltonian in case of periodic boundary conditions
For a one-dimensional lattice with L sites and periodic boundary conditions, the fermionic Hubbard model (1) has two SU(2) symmetries [S1]. Besides the obvious spin symmetry a so-called η-paring symmetry emerges if we look at the operatorŝ which obey the SU(2) commutation rules and satisfy the relationships meaning that any eigenstate of the Hubbard Hamiltonian H is also an eigenstate |η, η z ofη 2 andη z with eigenvalue η(η + 1) and η z , respectively. The presence of η-pairing states in the Hubbard model was first recognized by Yang [S2]. He showed that the states |φ Nη ∝ (η + ) Nη |vac , with |vac being the vacuum state and N η denoting the number of η-pairs, are eigenstates ofĤ which possess off-diagonal long-range order. Because these states are excited states, the long-range order does not show up in the ground state or thermal states ofĤ as shown in this study and also Ref. [S3].

Pair correlation function in case of open boundary conditions
For open boundary conditions (OBC), the pair correlation function in real space can be expressed as where the summation extends to the number of pairs of sites, N b = L − r, with sites separated by r in an open chain with L sites [S4]. According to [S5], the Fourier transform P OBC (q, t) = r e iqr P OBC (r, t), also shows the characteristic enhancement after irradiation which, however, will continue to grow for t > t 0 , in contrast to what is found for periodic boundary conditions. Most likely this is caused by the definition of the (quasi-)momenta in Eq. (S5), in particular at the boundaries. Instead of addressing this issue directly, one can use better the relation which holds for periodic boundary conditions. Figure S.1 shows the time-evolving block decimation (TEBD) results for P OBC (π, t) and η +η− in the halffilled Hubbard model, where we have parametrized the pump by A 0 = 0.38 and ω p /t h = 6.8 (this corresponds to the peak position of the A 0 -ω p contour plot at t · t h = 30, see Fig. 2 of Ref. [S5].) After pulse irradiation, for t · t h 12, the magnitude of 2 η +η− /L saturates to a constant value, reflecting the conservation law of the η-pair numbers, even though P OBC (π, t) is still weakly growing. Note that the relation (S6) will be fulfilled within an infinite TEBD (iTEBD) calculation: Here, P (π, t) saturates to a constant value, provided the bond dimension ia large enough, see Fig. 1(d) of the main text.
Nonlocal pairing correlation Pr>0(π, t) at T /t h = 8 In this section we are presenting further details about the on-site and nonlocal pair correlations at high temperatures. Figures S.2(a)-(d) give the contour plots of P (π, t) and P r>0 (π, t) for U/t h = 8 and T /t h = 8. At t 0 · t h = 10, both pair correlation functions are finite in a large range of parameter space, but the enhancement is rather weak compared to those at zero temperature. After pulse irradiation (t · t h = 15), the spectral intensity of P (π, t) and P r>0 (π, t) becomes concentrated in two spots with driving frequencies ω p close to the Hartree energy U/2 and to the Hubbard interaction energy U . Interestingly, the peak with ω p U (U/2) has positive (negative)  show P (π, t), P r>0 (π, t) and 2n d (t) at ω p ≈ U/2, marked by the white asterisk in Fig. S.2(d), where ω p = 4.2 and A 0 = 1.56. As in Fig. 4(c), P r>0 (π, t) becomes zero for long times, but now it approaches its limiting value from below. Thus double occupancy n d (t) dominates P (π, t) after pumping at high temperatures.

Accuracy of the iTEBD simulations
In Ref. [S5] we compared the results of TEBD simulations with OBC with exact diagonalization (ED) data, and demonstrated good agreement up to some time t · t h , depending on the maximum bond dimension χ max used. On the other hand, keeping the cutoff less than 10 −7 during the TEBD simulation for system size L = 12 and OBC, a perfect agreement with ED data is also achieved for t · t h 30, if χ max is of the order 10 4 . Since an analytical solution for the time-evolution of the infinite Hubbard model (1) after pulse irradiation is lacking, we discuss the accuracy of our iTEBD approach in the latter way. Figure S.3 presents iTEBD results for P (π, t) at zero temperature, obtained enforcing various maximum truncation errors. The discrepancies between the iTEBD data with maximum truncation errors 1 × 10 −7 and 1 × 10 −8 are negligible [see panel (a)], albeit the simulations have to be performed up to different t · t h = 16.0 and 14.4 because of the rapid increase of χ max needed [see panel (b)]. Performing the iTEBD simulations with a cutoff less than 10 −7 may not always be realistic, however, in view of limited computational resources. Fortunately, Fig. S.3(c) demonstrates that a reasonable accuracy can be obtained quite often using smaller χ max (see the results for χ max = 400).