Correlation-induced superconductivity dynamically stabilized and enhanced by laser irradiation.

Studies on out-of-equilibrium dynamics have paved a way to realize a new state of matter. Superconductor-like properties above room temperatures recently suggested to be in copper oxides achieved by selectively exciting vibrational phonon modes by laser have inspired studies on an alternative and general strategy to be pursued for high-temperature superconductivity. We show that the superconductivity can be enhanced by irradiating laser to correlated electron systems owing to two mechanisms: First, the effective attractive interaction of carriers is enhanced by the dynamical localization mechanism, which drives the system into strong coupling regions. Second, the irradiation allows reaching uniform and enhanced superconductivity dynamically stabilized without deteriorating into equilibrium inhomogeneities that suppress superconductivity. The dynamical superconductivity is subject to the Higgs oscillations during and after the irradiation. Our finding sheds light on a way to enhance superconductivity that is inaccessible in equilibrium in strongly correlated electron systems.

INTRODUCTION

After the highest superconducting critical temperature Tc was recorded in cuprate superconductors around 130 K decades ago, the highest Tc at ambient pressure has essentially remained unchanged (). One reason for this would be that the strong effective attraction between electronic carriers required for higher Tc unavoidably drives the tendency for charge inhomogeneities such as the phase separation in a mesoscopic scale and charge ordering (–). For this reason, materials may not be stably synthesized even if a sufficiently large effective attraction could be reached. Here, we pursue a way to circumvent this difficulty by using a new axis of control out of equilibrium.

Superconducting-like properties, such as gap formations, were reported above the room temperature in pump-probe laser measurements with the intension to excite coherent phonons in cuprates (, ). Such an approach was also applied to alkali-doped fullerides, and a similar phenomenon was observed above the critical temperature at equilibrium ().

Here, we focus on an alternative possibility ascribed to a mechanism more generic to strongly correlated electron systems. Our strategy is to control electron correlations by using the dynamical localization effect (–), where the electron hopping is effectively suppressed by the intense alternating current (ac) electric field at high frequency. The dynamical localization was recently observed experimentally (). We are also inspired by a recent finding that the effective attraction of carriers increases from intermediate to strong coupling regions in strongly correlated electron systems ().

By numerically calculating the time evolution of a standard correlated electron model for the cuprates—Hubbard model on square lattice containing the electronic transfer thop and the on-site Coulomb repulsion U—we show that the long-range d-wave superconducting correlation is substantially enhanced without causing charge inhomogeneity during a strong laser irradiation modeled by a time-dependent vector potential A(t) with the amplitude A0 and frequency ω (see Materials and Methods).

RESULTS

In tight-binding models, the ac laser irradiation causes a dynamical localization. Namely, the transfer integral thop is renormalized asat long time t ≫ 1/ω, where is the zeroth-order Bessel function (, ). If this reduction works in correlated electron systems as well, the effective interaction increases, because is satisfied until the first zero, , at A0 ~ 2.4 is reached (). The use of reduced hopping amplitude has been recently theoretically proposed for the purpose of enhancing superconductivity in different contexts (, ).

To analyze the model, we use the variational Monte Carlo (VMC) method (–). In particular, an innovation in the time-dependent VMC method for the fermion systems () has enabled the long-time nonequilibrium simulation of correlated electron dynamics in the present system (see Materials and Methods).

Ground state

Before showing our results for nonequilibrium dynamics, we briefly summarize the results of the ground states for the Hubbard model. Upon doping hole carriers to a half-filled antiferromagnetic Mott insulator, the superconducting order emerges with a dome-like order parameter amplitude as a function of the hole doping concentration δ around 0.1 in the intermediate to strong coupling regions (roughly for U/thop > 5) (, , , –), in essential agreement with the experiments in the cuprates, which is, however, strongly competing with charge inhomogeneous states such as phase separation and stripe order (–).

Figure 1 shows physical quantities of uniform and inhomogeneous superconducting states as functions of U/thop for the two-dimensional Hubbard model obtained by the VMC method. The detailed definitions of the physical quantities are given in Materials and Methods.

Fig. 1

Interaction dependence of physical quantities for superconducting states at zero-temperature equilibrium.

(A) Long-range average of the superconducting correlation in two-dimensional Hubbard model. The open and solid symbols represent the results for size L × L = 8 × 8, δ ≈ 0.06 and 16 × 16, δ ≈ 0.1, respectively. Green triangles (upside-down triangles) and blue circles represent the results of charge-uniform and inhomogeneous ground states (GS), respectively. The red square (diamond) symbols represent the results for charge-uniform metastable excited states (ES). The shaded regions are guides for the eyes. The blue and red arrows represent equilibrium path and hypothetical nonequilibrium process toward uniform excited states with enhanced superconductivity by the laser irradiation A(t), respectively. Error bars indicate the statistical errors arising from the Monte Carlo sampling, but most of them are smaller than the symbol sizes here and in the following figures. (B) Energy difference per site between the uniform state with strong superconductivity and the inhomogeneous state with the weak superconductivity for L = 8 and δ ≈ 0.06. Insets represent the charge distribution n of the uniform state for U/thop = 5 (left top), that of the inhomogeneous state for U/thop = 10 (right top), and the amplitude of the density inhomogeneity in the ground states (right bottom) respectively.

In Fig. 1A, we see that the long-ranged part of the d-wave superconducting correlation function for the charge uniform states shows a sharp crossover around U/thop ~ 5 (yellow region). Note that is basically the square of the superconducting order parameter. However, because the energy of the inhomogeneous state is slightly lower than that of the uniform state for U/thop ≳ 5 (see Fig. 1B), the uniform superconducting states are not the true ground states in the strong coupling region (blue arrow in Fig. 1A), indicating that the superconductivity is suppressed by charge inhomogeneities in equilibrium ground states. This is naturally understood, because charge inhomogeneous states can be regarded as superconducting regions bridged and separated by non- or suppressed-superconducting areas similarly to the Josephson network. This is consistent with previous numerical results that the phase separation (or charge inhomogeneities) occurs above U/thop ~ 6 and preempts most of the dome structure of superconductivity for δ ≲ 0.2 attained under the charge uniform condition (, ).

We add a note of caution on the quantitative aspects. The exact answer is not available even for the simple Hubbard model, and the quantitative range of the superconducting dome as well as the charge inhomogeneous region depend on details of the accuracy in numerical methods (–, –). However, the tight connection between the superconductivity and the charge inhomogeneities is a widely observed tendency (–, ) in accordance with the tendency in the experimental observations (). The strong-coupling superconductivity can be achieved only when the effective attractive interaction becomes strong, whereas too strong static attraction results in the abovementioned aggregation of electrons with charge inhomogeneity, which constitutes a “double-edged sword.”

Dynamically enhanced superconductivity

Figure 2 shows the dynamics of induced by a strong laser with the scaled frequency ω/thop = 15 simulated after the rising time tc = 60. Here, we confirmed that the long-ranged part of the superconducting correlation becomes nearly independent of the distance and reaches an enhanced nonzero constant after the rising time, supporting the emergence of the d-wave long-ranged superconducting order (see Materials and Methods). The initial state is a weakly superconducting ground state for U/thop = 5. This ratio U/thop is close to the values for the cuprates estimated by the recent first-principles calculations (). By the laser irradiation, we see that the time-averaged for t ≥ tc is substantially enhanced from the initial state. The averages are consistent with those of the charge-uniform excited states for the expected effective interaction, U/thop ~ 5.78 for A0 = 0.75 and ~6.53 for A0 = 1.0 inferred from Eq. 1. See also figs. S1 and S2 for results of different system sizes.

Fig. 2

Time dependence of induced by strong laser in the two-dimensional Hubbard model.

Parameter set is (δ, L, U/thop, ω/thop, t) = (0.0625, 8, 5.0, 15, 60) on an 8 × 8 lattice. Red circles and green squares represent the results for A0 = 1.0 and A0 = 0.75, respectively. Red solid and broken green lines indicate corresponding values of obtained from VMC for the charge-uniform superconducting states by assuming the effective enhancement of U/thop. Black dash dotted line is the equilibrium value in the absence of the laser irradiation. Insets represent the charge distributions n(r) for A0 = 1.0 at t = 130 and 210, corresponding to a peak and a bottom of the oscillation, respectively.

These results strongly support that this enhancement is caused by an intrinsic effect of the dynamical localization beyond the finite-size artifact. Other physical quantities also consistently indicate that the laser-driven system is well represented by the corresponding uniform excited state at the effective interaction (see section SB and figs. S3 to S5).

More importantly, this enhancement is numerically achieved without causing substantial charge inhomogeneities in contrast to the equilibrium ground state (see insets in Fig. 2). The enhancement without the inhomogeneity is not very sensitive to the irradiation process and may partly be associated with larger transition elements of the irradiation term between charge homogeneous states than those between homogeneous and inhomogeneous states, in addition to the fast evolution of the pairing. The laser irradiation dynamically stabilizes a charge-uniform state with enhanced superconductivity that is not possible in the equilibrium ground state. We discuss the mechanism of the dynamical stabilization later in the paper.

Higgs oscillation

, which is essentially the squared superconducting order parameter, shows long-period oscillations in time in Fig. 2. The oscillations continue after the laser irradiation (see section SC and figs. S7 to S9), signaling the Higgs amplitude oscillation mainly studied in Bardeen-Cooper-Schrieffer or phonon-mediated superconductors with the s-wave symmetry (–). Depending on the irradiation process, the system may fall into a linear combination of several charge-uniform excited states at the effective , which have different superconducting correlations.

Figure 3 shows the optical conductivity of the uniform superconducting state for U/thop = 5.78 obtained by the VMC method combined with a procedure used in the study by Shao et al. () as detailed in Materials and Methods. A shoulder around ω/thop ~ 0.09 coexisting with the delta function at ω/thop = 0 smeared by an artificial damping η = 0.04 is in accordance with the inverse oscillation period of the squared superconducting order for A0 = 0.75 seen in Fig. 2, which is consistent with the interpretation by the gapped Higgs mode.

Fig. 3

Real part of optical conductivity of the superconducting state.

The parameters are L = 8, δ = 0.06, and U/thop = 5.78. The results (red curves with symbols) are obtained by using the VMC method combined with a numerical procedure used in the study by Shao et al. (32). Error bars indicate the statistical errors arising from the Monte Carlo sampling. The inset shows a wider range of ω dependence. The blue solid line is a fitting obtained by a sum of two Lorentzian curves (dashed and dotted lines).

DISCUSSION AND SUMMARY

For an intuitive and deep understanding of the dynamical stabilization of the superconductivity, further studies are needed. However, even at this moment, a possible mechanism of the dynamically stabilized superconductivity consistent with our results can be argued in the following way. In the driven electron systems by laser, carriers coupled to electric fields are shaken in the direction of laser irradiation. Because the frequency of the vector potential modulation is larger than the electron hopping thop, electrons partially follow the modulation, and the retardation causes weakening of the stripe charge density modulation by the motion of electrons from electron-rich to electron-poor regions, in addition to the dynamical localization effect. The shaking also increases the double occupation of electrons especially around the point of peak absolute value of the electric field oscillation as one can see in Fig. 4. In the equilibrium condition, the ground state has a larger amplitude of stripe modulation and a smaller double occupation (by an increment of ΔD ~ − 0.0007) than the strongly superconducting metastable state at U/thop ~ 6.53. Because the energy difference between these two states is tiny (at most less than 0.01thop), the laser irradiation therefore plays the role of a bias field to favor the metastable state over the ground state. In fact, the increment in the double occupation at the peak (bottom or top) of the electric field oscillation is comparable to the above increments of the superconducting metastable state relative to the stripe ground state in equilibrium. Then, the irradiation is expected to dynamically stabilize the uniform superconducting state, which is only metastable when the laser irradiation is absent.

Fig. 4

Time dependence of double occupation in short time scale around t = 100 for the same case as Fig. 2 for A0 = 1.0.

The black curve shows time evolution of electric field . The dotted line represents the zero line of the electric field.

Next, we discuss realistic possibilities of enhancing superconductivity in the cuprates by using the present mechanism. The most striking enhancement of the superconductivity is expected if U/thop of real materials is in the crossover region in equilibrium. However, despite first-principles calculations (), it is not yet clear which region of U/thop is adequate for the cuprates when the single-band Hubbard model representation is applied, because the derivations of low-energy effective Hamiltonians call for further improvements such as in the disentanglement procedure of bands (, ).

In our example of A0 = 0.75 and ω/thop = 15, by using realistic values of the lattice constant a ~ 0.4 nm and thop ~ 0.45 eV, the maximum intensity of the electric field Fmax = A0ω is required to be ~120 MV/cm, which is relatively strong in terms of so far achieved intensity. One way to alleviate the intensity condition would be to use short-width pulses with small frequencies, which contains high-frequency components because of the energy-time uncertainty relation.

We comment on the heating effect possibly induced by the laser irradiation. Contrary to the quench systems where the hopping parameter is suddenly changed from the initial value, the laser irradiation has a possibility of continuously heating up systems in general. However, in the present calculations, such heating is negligible within our simulated time scale because the time-averaged energy does not change during the laser irradiation, as seen in fig. S3. This is related with the fact that the single-band Hubbard model with U/thop < 10 does not have electronic density of states at ω/thop ≥ 10. We should be careful to treat the heating effect if we use lasers with small frequencies below ω/thop ~ 10. Note that although the time-averaged energy is increased by laser irradiation from about −0.8 to −0.5 in fig. S3, most of this increase is accounted by the increase in the ground-state energy for the system with the effectively reduced hopping because it causes a substantial increase in the kinetic energy. Therefore, during the laser irradiation, this increased energy has no way to relax.

Aside from technical problems, we need to avoid the heating of the system caused by the dissipation to lattice degrees of freedom. However, this will be delayed and we may be able to see the enhancement before the heating because the time scale of phonons in most materials is above picoseconds, which is much slower than that of electrons. In fact, for thop = 0.45 eV, our results show that the enhanced superconductivity is stable during the irradiation from roughly 100 fs after the rising time as seen in Fig. 2. In addition, the excitation energy ΔE from the ground state at the renormalized hopping is only a few kelvins for thop = 0.45 eV, as seen in Fig. 1B, and the averaged energy is conserved during the laser irradiation as mentioned above (see section SB and fig. S3). Therefore, the energy will hardly be relaxed to the lattice degrees of freedom during the irradiation. It is an intriguing future issue to examine longer-time sustainability of the superconductivity by tuning the repeated switching on and off process of the irradiation to suppress the heating.

Another problem is the multiband structure in real materials missing in the Hubbard model we considered. In our calculations, no optical transition is allowed because the frequency of light is in the order of 10, as mentioned above. However, in the realistic condition, interband transitions may happen because of multiband structures in materials, which we ignored. It may be important to use the frequency at which the optical transition is prohibited or at least the frequency at which the transition is allowed only between the energies with small density of states. More detailed studies on the heating by laser and the effects of multiorbital structure are left for future work.

The present approach is also useful to clarify the mechanism of high-Tc superconductivity in equilibrium, which is still an open issue. For instance, we can distinguish whether the U/thop of real materials, such as the cuprates and the iron-based superconductors, is in the weak or strong coupling regions from irradiation effects.

Finally, we comment on the possibility of enhancing s-wave superconductivity in alkali-doped fullerenes where light-induced superconductivity was proposed similarly to the cuprates (). In this system, negative Hund’s rule coupling arising from Jahn-Teller phonons was proposed to play an important role in the emergence of superconductivity (, ). It is an intriguing future issue to study whether the reduction of hopping integrals by the dynamical localization and the resultant relative increase of interactions enhances superconductivity.

In summary, in a model of laser irradiation to strongly correlated electron systems, we found dynamically stabilized and enhanced d-wave superconducting states that do not suffer from inhomogeneity. It will inspire further studies on high-Tc superconductivity that is not possible in equilibrium.

MATERIALS AND METHODS

Model

The Hubbard model is defined bywhere thop and U represent the hopping amplitude of electrons between nearest-neighbor sites and the on-site interaction, respectively. The creation and annihilation operators of electrons at site i with spin σ are expressed as and c, respectively. In the presence of a vector potential, the hopping amplitude is modified as thop → thop exp[iA(t) ⋅ (r − r)], where A(t) is a vector potential at time t and r is the position vector at site i. We considered a square lattice and assumed that the vector potential is polarized along the direction d = (1, 1), that is, A(t) = A(t)d. The time dependence of the vector potential A(t) was taken as

We set t = 3t and used thop and as units of energy and time, respectively.

Definitions of physical quantities

To study physical properties in equilibrium and nonequilibrium, we measured the energy per site E/N, its average over an ac cycle , double occupancy D, amplitude of charge density Δn, charge structure factor N(q), and the -wave superconducting correlations Pd(r) defined asrespectively. Here, nmax(nmin) are those among {n, i = 0,⋯, N} and n = N/N, where N is the number of electrons. 〈Q〉 represents the expectation value of an operator Q, that is, where |ψ〉 represents a wave function. Δd(r) represents the -wave superconducting order parameter defined aswhereis the form factor, which describes the -wave symmetry, r = (r, r), and δ = 1 for k = l and 0 otherwise is the Kronecker’s delta. We also defined the averaged value of long-range parts of the superconducting correlation aswhere M is the number of vectors satisfying r = |r| > 5. In our calculated data, Pd(r) saturates to a nonzero constant if the superconducting order exists (fig. S6).

VMC method for correlated electron systems

We used a VMC method to obtain nonequilibrium states and ground states in strongly correlated electron systems (–, ). As with the VMC method for ground states, the recently developed time-dependent VMC method for fermion systems provided us with an accurate and efficient way to study nonequilibrium dynamics of correlated electrons in finite-dimensional systems without the negative sign problem (). Because of this method, we were able to simulate long-time electron dynamics for two-dimensional carrier-doped correlated systems, whose sizes can be taken sufficiently large to estimate the stabilities of various competing long-range orders on equal footing.

As a trial wave function, we adopted the generalized pair product wave function with correlation factors , where and is the correlation factor. Here, f are site-dependent variational parameters as noted below. The correlation factor consists of the Gutziwiller, Jastrow, and doublon-holon binding factors, which are defined asandrespectively. is 1 when a doublon (holon) exists at the ith site and m holons (doublons) surround at the lth nearest neighbor. Otherwise, is 0. Here, we treated g, v, , and f as variational parameters. Because many variational parameters were introduced not only to correlation factors but also to the pair product part, this type of trial wave functions can flexibly and accurately describe quantum states in correlated electron systems (fig. S2) (, , ). We note that our trial wave function reproduces the exact results for the time evolution of physical quantities well, such as energy, momentum distribution, spin structure factor, and superconducting correlations in the Hubbard model (). All the variational parameters were optimized simultaneously by using the stochastic reconfiguration (SR) method for imaginary-time or real-time evolution based on the time-dependent variational principle (–, –) to be described later in detail.

In actual calculations, we used a trial wave function that has an 8 × 2 sublattice structure by introducing 8 × 2 × N variational parameters for f, because it enables flexibly describing realistic and possible inhomogeneous states, that is, the inhomogeneous state is allowed to emerge dynamically. Note that because uniform-charge states do not exist for the equilibrium ground state above U/thop ~ 6, we imposed a 2 × 2 sublattice structure on a trial wave function for U/thop ≳ 6 to examine the charge-uniform lowest excited states.

Real-time and imaginary-time evolution in VMC method

In our VMC calculations, ground states and nonequilibrium states were obtained by using the SR method (–). This method can be regarded as the imaginary-time or real-time evolution of a trial wave function in the truncated Hilbert space based on the time-dependent variational principle (–). In the case of real-time evolution, by minimizing a distance between and ℋ|ψ〉 in the time-dependent Schrödinger equation in accordance with the time-dependent variational principle (–), the equation of motion for variational parameters α = {α|k = 1,⋯, N} was obtained aswhere and (, , , ). Here, we used the fourth-order Runge-Kutta method to solve Eq. 13. The imaginary-time evolution of variational parameters was obtained by replacing the real-time variable t in Eq. 13 with the imaginary-time τ = it. If we solve Eq. 13 for the imaginary-time evolution by using the Euler method, it is nothing but the conventional SR scheme for ground states (, ). In the VMC method, the quantities S and g were estimated by the Markov chain Monte Carlo method.

Optical conductivity

We also calculated the optical conductivity measurable by a pump-probe experiment by using a laser with weak intensity (). The details of the method are as follows: First, we obtained the time-dependent trial wave function in the system under a laser irradiation with weak intensity Aprobe. Here, Aprobe represents the probe laser polarized along the x direction defined as

Next, we calculated the expectation value of the current operator defined bywhere a is the lattice constant. To obtain the optical conductivity in Fig. 3, we set t = 0.05, t = 0.01 and A0 = 0.05. Finally, we calculated the optical conductivity σ defined aswhere the electric field is defined as . Note that the Fourier transformations in Eq. 16 is defined bywhere we set sufficiently large T and small η as T = 200 and η = 0.04.