Charge Oscillations Emerging after Application of an Intense Light Field to Superconductors on a Dimer Lattice

Motivated by recent experimental findings for an organic superconductor, charge oscillations that emerge after a strong pulse of an oscillating electric field is applied are studied in electron systems in a superconducting phase on a lattice with a dimerized structure. They are analyzed using Fourier spectra of charge densities whose time profiles are obtained by numerically solving the time-dependent Schr\"odinger equation within a mean-field approximation. Depending on the strengths of attractive interactions, different charge-oscillation modes appear. For weak attractions, the charge-oscillation mode is an electronic breathing mode, which was previously found for repulsive interactions. For strong attractions, it is a pair analog of the electronic breathing mode (a `pair' breathing mode). For intermediate attractions, it is another mode whose transient current distributions are considerably different from those of the breathing modes. We investigate how their frequencies and amplitudes depend on interactions and transfer integrals.


Introduction
In solid-state materials, various cooperative and dynamic phenomena are known. Photoinduced phase transitions are such phenomena, and they are achieved between different phases and on different timescales. [1][2][3][4][5] In many cases, electronic orders are melted or weakened.
However, in some cases, orders are transiently constructed or enhanced. [6][7][8][9][10][11] As a possible mechanism for the enhancement of an electronic order, dynamical localization may work, [12][13][14] which suppresses the itinerancy of electrons and enhances the relative importance of interactions during and even after photoexcitation. [15][16][17][18][19][20][21] In some cases, however, it has been shown that dynamical localization is insufficient and electron correlations are necessary to explain experimentally observed details in the enhancement. 7,[22][23][24] For excitonic insulators, 11) the excitonic order has been theoretically shown to be enhanced by dif-ferent mechanisms. [25][26][27] For photoinduced superconductivity, 6) electron-phonon interactions are considered to play an essential role, but their effects are still controversial. 28,29) If the coherence among pairs is reduced, conductivity is suppressed. 30) Thus, in many cases, electron motion needs to be driven coherently by an oscillating electric field. 31,32) Among the photoinduced phenomena, those that are realized after an intense light field is applied to solid-state materials have attracted much attention. Resultant nonequilibrium states are not simply far from equilibrium ones but may show some distinctive dynamic behaviors.
High-harmonic generation, which is now being extensively studied experimentally [33][34][35][36] and theoretically, [37][38][39][40] is such a phenomenon, although it appears during photoexcitation and is regarded as a time-periodic steady-state response. In some cases, the synchronized motion of electrons has been suggested to be caused by interactions during 39) and after 39,41,42) an intense light field is applied. This fact is reminiscent of the situation in a discrete time crystal, 43,44) which is realized in many-body-localized periodically driven systems: 45,46) interactions are essential for collective synchronization in strongly disordered systems. 44) Quite recently, a nonlinear charge oscillation and a resultant stimulated emission have been observed in the organic superconductor κ-  2 ]Br], 42) whose lattice has a dimerized structure. This and related materials are known for their phase diagram with dimer-Mott-insulator-metal and superconducting transitions and an unconventional critical behavior. 47) Photoinduced insulator-metal transitions are also known. [48][49][50] Because of the dimerized structure, its dielectric permittivity shows an anomaly 51) owing to polar charge distributions. [52][53][54][55][56] Charge fluctuations associated with such distributions may be relevant to the superconducting phase transition. 57,58) For the nonlinear charge oscillation to appear, a dimerized structure has been shown to be important by theoretical calculations on the basis of the exact diagonalization method. 41) However, the relationship between this oscillation and superconductivity has not yet been theoretically discussed. Experimentally, it has been observed that the resultant increase in reflectivity is enhanced by superconducting fluctuations and weakened as the temperature decreases below the superconducting transition temperature. 42) Thus, it is important to clarify their relationship theoretically, which would be useful in characterizing the superconducting phase in organic conductors.
Here, we theoretically study what kind of charge oscillations appear after an intense electric-field pulse is applied to superconducting states on a dimer lattice of a simpler form.
To treat the time profiles of charge densities and pairing order parameters, we employ the 2/23 Hartree-Fock-Gor'kov approximation. In general, without dimerization, time-dependent pairing order parameters have been extensively studied experimentally 59) and theoretically. [60][61][62][63][64] Dynamic relationships with charge density waves have also been studied. 60,64) On the other hand, without superconductivity, the synchronized motion of electrons has also been theoretically realized on a honeycomb lattice, 39) which also has two sites per unit cell. In this paper, we show that novel charge-oscillation modes appear through a dimerization-induced coupling between superconductivity and charge-density modulation after strong photoexcitation.

Model with On-Site and Intradimer Interactions
In a previous study, 41) the nonlinear charge oscillation was shown to appear in a onedimensional spinless-fermion "t 1 -t 2 -V" model at half filling and a two-dimensional extended Hubbard model for κ-(BEDT-TTF) 2 X at three-quarter filling. It is briefly mentioned that the appearance is not limited to these fillings. To obtain a hint on how general these models are, which produce nonlinear charge oscillations, and to study the effect of superconductivity on them, we consider a toy model whose lattice structure is much simpler than that of κ-(BEDT-TTF) 2 X, where the kinetic term H kin is a tight-binding model on the dimer lattice shown in Fig. 1(a) with Peierls phase factors [ ea c A x(y) with A x(y) being the x(y) component of the vector potential and a being the lattice constant is simply written as A x(y) below], Here, c † i,α,σ creates an electron with orbital α and spin σ at dimer site i. The intradimer transfer integral t 1 and the interdimer ones t 2 and t y are configured as shown in Fig. 1 where e 12 (k, A) is defined as e 12 (k, A) = t 1 e iA x /2 + t 2 e −ik x −iA x /2 + 2t y e −ik x /2 cos(k y + A y ) .

4/23
For the interaction term H int , we consider on-site and intradimer interactions, where We employ the Hartree-Fock-Gor'kov approximation and replace H int by H MF int , assum- and c i,α,↓ c i,β,↑ = ∆ αβ .
Using V ± ≡ 1 2 V sp ± V ch ,ᾱ = 2 for α = 1, andᾱ = 1 for α = 2, H MF int is written as where and The order parameters are calculated using the mean-field Hamiltonian H MF = H kin + H MF int with the formulae and which are numerically iterated until convergence in obtaining the ground state. For the model The chemical potential µ is set so that the system is at three-quarter filling (i.e., three electrons per dimer) throughout this paper. Even at half filling, which is insulating without interactions, nonlinear charge oscillations appear. To study the effect of superconductivity on them, however, the filling is set to be away from half filling. For the system size, we use N=100×100 and periodic boundary conditions. This is sufficiently large for the spectra and 6/23 time profiles shown in this paper: the finite-size effects are smaller than the symbols.
The initial state is the mean-field ground state. For U < 0, it is always a superconducting state. Photoexcitation is introduced through the Peierls phase, where different time-dependent vector potentials are used depending on the purpose. For absorption spectra, we use slowly decaying oscillating electric fields with a small amplitude, 66) where ω is the frequency and γ is the decay constant. For photoinduced charge oscillations and associated Fourier spectra, we use symmetric one-cycle electric-field pulses, 18,19,23,41) where the central frequency ω c is chosen to be ω c = 0.7 throughout the paper because the qualitative results are independent of its value as in the previous study. 41) In both cases, the maximum electric field F is polarized along (1, 1), F = (F, F), unless stated otherwise. The time-dependent Schrödinger equation is numerically solved by directly extending the method that uses the Hartree-Fock approximation. 26,[67][68][69] In contrast to the previous study, 41) every dimer is always equivalent for any polarization of F as assumed in Eqs. (6)- (8). In the ground state, the system satisfies σ n 1,σ = σ n 2,σ .
Photoexcitation with F x 0 leads to σ n 1,σ σ n 2,σ . Thus, the time profile of σ n 1,σ (or σ n 2,σ ) has information that is shared with the optical conductivity or absorption spectrum. Then, as in the previous study, 41) we calculate the absolute values of the Fourier transforms of these time profiles and refer to them as Fourier spectra. The time span used for the Fourier spectra is mainly T < t < 50T with T = 2π/ω c , but we use T < t < 500T when oscillation frequencies with a higher resolution are required.
To facilitate the characterization of charge oscillations, we calculate current order parameters, which are defined as which are independent of i. They are related to the current density

Absorption spectra
The mean-field Hamiltonian H MF is described by a 4 × 4 matrix at each k, so that two Because the Higgs mode is not excited in the first order with respect to A, 63) it does not appear in a linear absorption spectrum. Then, we make the field amplitude ten times larger (F = 0.02) and show the resultant spectrum in Fig. 2(a). Because the energy supplied by the field basically becomes 100 times larger, the spectrum for F = 0.002 is multiplied by 100 for comparison. Since the difference between the lowest unoccupied and highest occupied energy levels is 0.153, the so-called Higgs mode is responsible for the peak at ω = 0.153 for the larger F. In fact, for the initial state without imaginary parts of ∆ αα , the photoexcitation with this frequency oscillates not only the real parts of ∆ 11 and ∆ 22 (in phase) but also their imaginary parts (antiphase) and the charge-density difference because the latter is coupled to the pairing order parameters for nonzero dimerization. Because the imaginary parts of ∆ αα are smaller than the real parts, their magnitudes mainly oscillate with this frequency. In this sense, this mode is regarded as the Higgs mode. However, we will no longer discuss the Higgs mode since the absorption due to the Higgs mode is generally very weak and it barely affects charge oscillations.
Those for U = −0.4 are shown with different F in Fig. 2(b). For small F (F = 0.1), it has spectral weights at energies where the linear absorption spectrum [corresponding to the smaller-F case in Fig. 2(a) except for the tail continuing from zero energy due to γ > 0] has spectral weights. The lower-energy peak at ω ≃ 0.4 is conspicuous, and small but nonzero 9/23 weights exist in the range of 0.6 < ω < 1.2. Because | ∆ αα | is coupled to σ n 1,σ − σ n 2,σ for nonzero dimerization, | ∆ αα | also oscillates with ω ≃ 0.4. The conspicuous peak becomes slightly blueshifted and lowered with increasing F (from F = 0.1 to 0.5). It is suppressed for large F. That of ω = 0.54 also appears in the Fourier spectrum of | ∆ αα | (not shown); thus, | ∆ αα | is coupled to σ n 1,σ − σ n 2,σ in this mode. It is referred to as a middle-frequency chargeoscillation mode. These modes are analyzed later in detail. The charge-oscillation mode of ω = 0.33 also appears in the Fourier spectrum of | ∆ αᾱ | (not shown); thus, | ∆ αᾱ | is coupled to σ n 1,σ − σ n 2,σ in this mode. Its weight in the Fourier spectrum of σ n α,σ is relatively small as long as the on-site attraction | U | is a dominant interaction. The main results for large F here and below are independent of the polarization of F as long as F x 0. This is in contrast to the two-dimensional case in the previous study, 41) where every dimer becomes equivalent only for specifically polarized fields.

Interaction-strength dependence
Hereafter, we focus on charge-oscillation modes that appear for large F and observe their behaviors with different model parameters. Fourier spectra for weak, intermediate, and strong on-site attractions are shown in Figs. 3(a)-3(c), respectively. For weak interactions, the highfrequency charge-oscillation mode always appears at ω = 1.2. Its peak height decreases with increasing | U |. For the middle-frequency charge-oscillation mode, its frequency and peak height increase with | U | as shown in Fig. 3(a). The charge-oscillation mode coupled to | ∆ αᾱ | appears at ω = 0.33 for weak to intermediate on-site attractions. However, it does not dominate the Fourier spectra of σ n α,σ . This mode will not be discussed hereafter.

For intermediate on-site attractions, the high-frequency charge-oscillation mode is almost
invisible. In addition to the middle-frequency charge-oscillation mode, whose frequency increases with | U | and peak height shows a maximum as a function of | U |, a new mode appears on its low-energy side in Fig. 3(b). It is referred to as a low-frequency chargeoscillation mode. For this mode, its frequency decreases and its peak height steeply increases from zero as | U | increases; thus, this mode is not split from the middle-frequency charge-   For strong on-site attractions, only the middle-frequency and low-frequency chargeoscillation modes are visible in the energy range shown in Fig. 3(c). Here, their frequencies are separated further as | U | increases. The low-frequency charge-oscillation mode is now dominant, and its peak height shows a maximum as a function of | U |. Note that the present mean-field approximation always produces narrow peaks in Fourier spectra, whose widths are comparable to the frequency slice; thus, the peak heights shown in Fig. 4 understood as a guide. In the following subsections, we show the characteristics of the chargeoscillation modes in Fig. 4 one by one.

High-frequency charge-oscillation mode
For U = −0.1, where the high-frequency charge-oscillation mode is dominant in the Fourier spectrum, the time profiles of the current order parameters defined in Sect. 2 are shown in Fig. 5(a). Their colors (darkness) hereafter correspond to those of the bonds and the arrows in Fig. 1. Note that their plots are not shifted vertically: their centers of oscillations are nonzero. The system size used here is sufficiently large to ensure that this is not a finite- 13/23 size effect. The current density induced by the first half of the one-cycle electric-field pulse ( j x1 , j x2 , j y1 , j y2 < 0) is not completely canceled by its second half even after time averaging.
This means that the electrons drift after photoexcitation. The drift of electrons occurs only for weak interactions. As shown later, the centers of oscillations of the current order parameters are zero for intermediate to strong interactions. Here, the oscillating components of j x1 and j y1 are opposite in sign to those of j x2 and j y2 . Thus, the current order parameters oscillate back and forth between the left and right panels of Fig. 1(b). This mode was called an electronic breathing mode in the previous study. 41) These patterns of current distributions maximize the charge-density difference | σ n 1,σ − σ n 2,σ |. As derived in the previous study, 41) the frequency of the electronic breathing mode is given by the relation which is independent of U, and takes a value of 1.2 in the present case.
In Fig. 4(b), the peak height of the electronic breathing mode is shown to decrease with increasing | U |. This fact implies that the charge-density difference | σ n 1,σ − σ n 2,σ | competes with the pairing order parameter | ∆ αα | during the charge oscillation. We vary U and plot the peak height as a function of | ∆ αα | in Fig. 5(b). In addition, we fix U = −0.4, vary either V sp , V ch , or V ph from zero, and plot the peak height in the same figure. Note that | ∆ αα | is increased in the cases of V sp > 0 (antiferromagnetic spin-spin interaction), V ch > 0 (attractive charge-charge interaction), and V ph < 0 (where ∆ 11 cooperates with ∆ 22 ). In all cases, the peak height of the electronic breathing mode decreases as | ∆ αα | increases, which makes the competition clear.

Low-frequency charge-oscillation mode
For U = −2.0, where the low-frequency charge-oscillation mode is dominant in the Fourier spectrum, the time profiles of the current order parameters are shown in Fig. 6(a). For (b), the time span T < t < 500T is used for the Fourier spectra with higher resolution. Also shown is the second-order estimation with respect to transfer integrals [Eq. (27)], which is explained in the text.
| σ n 1,σ − σ n 2,σ |. The current order parameters oscillate back and forth between the left and right panels of Fig. 1(b). Because of the strong on-site attraction, an electron with σ =↑ is always accompanied by one with σ =↓.Thus, this mode is regarded as a pair breathing mode or a bipolaronic breathing mode.
For even larger | U | or smaller transfer integrals, the slowly varying component becomes even slower and the breathing motion becomes more evident in the time profiles of the current order parameters (not shown). In the limit of small transfer integrals, we can estimate this frequency. Two electrons with opposite spins are tightly bound and transferred in the second 15/23 order with respect to transfer integrals; thus, their effective transfer integral is given by the for b = 1, 2, and y. The factor 2 comes from the two second-order processes (an electron with spin up first or an electron with spin down first). However, in the present approximation, one-body wave functions with spin up and spin down are always degenerate, and the factor 2 is lost. In the strong-coupling limit, the difference between the lowest unoccupied and highest occupied energy levels in the mean-field ground state 2∆ becomes | U |. For large | U | but not in this limit, the estimation becomes better if the energy denominator | U | is replaced by 2∆. Thus, in the present approximation, we use for b = 1, 2, and y. In Fig. 6(b), we plot the frequency of the pair breathing mode and its estimation on the basis of t eff,MF b (b = 1, 2, and y) They coincide with each other in this limit. It is confirmed that when the exact diagonalization method is used for a 12-site chain in the present model with t y = V sp = V ch = V ph = 0 and periodic boundary conditions, the pair breathing mode appears after strong excitation for large | U | and its frequency is given by in the strong-coupling limit (not shown). The relationship between the frequencies determined from the exact-diagonalization spectra and the numerical values given by Eq. (28) is similar to that shown in Fig. 6(b). The appearance of the pair breathing mode in the exact diagonalization study implies that the symmetry breaking accompanied by superconductivity is not required for its appearance. for j x1 are reversed in Fig. 1(b) that describes the breathing modes. It does not maximize | σ n 1,σ − σ n 2,σ | but it wastes charge motion. The increase (decrease) in the local charge 16/23 density caused by j x1 is partially canceled by j x2 , j y1 , and j y2 . Sometimes the relative signs are not like this, but j x1 almost always has the same sign as j x2 , and j y1 and j y2 almost always have opposite signs for any time domain after photoexcitation. Now, we show the dependence of the frequency of the middle-frequency chargeoscillation mode on transfer integrals. The frequency increases as the magnitude of the intradimer transfer integral | t 1 | increases [ Fig. 7(c)], while it decreases as the magnitude of the interdimer transfer integral | t 2 | or | t y | increases [although | t 2 | and | t y | are simultaneously varied in Fig. 7(d)]. Similar dependences are also obtained for U = −0.4 and U = −0.7 (not shown). These facts are consistent with the fact that the increase (decrease) in the local charge density caused by j x1 (t 1 process) is partially canceled by j x2 (t 2 process), j y1 , and j y2 (t y processes). Note that as | t 2 | and | t y | decrease, the frequency of the electronic breathing mode (the so-called high-frequency charge-oscillation mode) decreases and that of the present mode (the so-called middle-frequency charge-oscillation mode) increases. Consequently, the frequency of the former can be lower than that of the latter, contrary to their naming.

Conclusions and Discussion
Following the previous study, 41) where a high-frequency charge-oscillation mode is shown to appear after an intense electric-field pulse is applied to electron systems on dimer lattices including those of organic conductors κ-(BEDT-TTF) 2 X, we have searched for charge-oscillation modes that appear in superconductors under similar conditions. Using For weak attractions, the charge-oscillation mode has a high frequency and is regarded as an electronic breathing mode, whose frequency is already known as a function of transfer integrals for on-site interactions. 41) Its amplitude becomes smaller when the pairing order parameter is increased by modifying interactions. This is due to the competition between the charge-density difference and the pairing order parameter. Although the present pairing is swave and interaction strengths are varied here, this fact seems consistent with the fact that the corresponding photoinduced increment in reflectivity decreases as the temperature decreases below the superconducting transition temperature. 42) For strong attractions, two electrons with opposite spins move together. Consequently, the 17 Although the two new charge-oscillation modes are concerned with s-wave pairing, we expect similar charge-oscillation modes for other pairings in addition to the electronic breathing mode if the lattice has a dimerized structure. If such a mode is experimentally observed, it would optically contribute to the characterization of superconductivity, e.g., where it is located in the BCS-BEC crossover. As discussed already, 39,41) the emergence of a strongfield-induced charge oscillation is regarded as a synchronization phenomenon. This aspect will also be studied in the future. (a)