Mechanism for Synchronization of Charge Oscillations in Dimer Lattices

We discuss the mechanism and the conditions for the appearance of synchronized charge oscillations which have been observed experimentally and theoretically after strong photoexcitation of dimerized systems. In the Hubbard model with on-site repulsion, the Bloch equations for a wave-number-dependent pseudospin -- whose components describe the charge-density difference, current density, and bond density between the two sublattices -- involve an alternatingly tilted pseudomagnetic field, which assists the synchronization of pseudospins with different wave numbers, irrespective of the initial condition. This fact is numerically confirmed by the dynamics in finite lattices based on the exact diagonalization method. In the presence of nearest-neighbor repulsion, however, the synchronization can be hindered by excitons. Therefore, the excitation of a sufficiently large density of free electron-hole pairs, but low density of excitons, is needed to achieve synchronization.


Introduction
Among various nonequilibrium phenomena, the manyelectron dynamics in far-from-equilibrium lattice systems has attracted much attention. [1][2][3][4][5][6][7] Many interesting effects, including photoinduced phase transitions, are nonlinear and caused by interactions.  With progress in the experimental techniques, ultrafast phenomena that necessarily involve highenergy processes become increasingly important. Which electronic phase is realized in a correlated electron material is often determined by a subtle competition between the kinetic and interaction terms in the Hamiltonian. Thus, low-energy properties are linked to high-energy processes, which makes ultrafast control possible in some cases.
During photoexcitation, high harmonic generation is observed in solids and reflects properties of their electronic states. [32][33][34][35][36][37][38][39][40] Recently, second harmonic generation has been observed in a centrosymmetric organic superconductor κ-(bis[ethylenedithio]tetrathiafulvalene) 2 Cu[N(CN) 2 ]Br [κ-(BEDT-TTF) 2 Cu[N(CN) 2 ]Br] through a nonlinear petahertz current before substantial scattering processes occur. 41) In this compound, [42][43][44][45][46][47][48][49][50][51] which is known to have a dimer lattice, stimulated emission is observed after strong photoexcitation. 52) This stimulated emission is also caused by a nonlinear charge oscillation -an electronic breathing mode -as has been shown by numerical calculations based on the exact diagonalization method. 53) When the photoexcitation is weak, i.e., when the optical field amplitude is small, different charge oscillations appear whose frequencies correspond to the peak energies in the optical conductivity spectrum. When photoexcitation is strong, however, these charge oscillations synchronize to produce the electronic breathing mode. The synchronization can be demonstrated numerically by introducing randomness into the transfer integrals to suppress it and showing that sufficiently strong on-site repulsion overcomes the randomness and recovers the effect. 54) It has also been discussed that strong on-site attraction produces a pair analog of the electronic breathing mode. 55) A widely studied lattice with two atoms in a unit cell is the honeycomb lattice, on which a synchronization transition and * E-mail: kxy@phys.chuo-u.ac.jp resultant coherent oscillations have been reported to occur in a recent mean-field investigation. 39) Here, it was pointed out that the equations of motion are somewhat similar to those in the Kuramoto model, [56][57][58] where the presence of a synchronization transition is established. This analogy is highly suggestive, but it implies that whether or not the in-phase synchronization is realized depends on the sign of the interaction. Hence, the detailed form of the equations of motion merits further investigations. In the present work, we will adopt a pseudospin representation and show that the synchronization is indeed sensitive to the sign of the interaction, a conclusion that will be confirmed numerically using the exact diagonalization method. We hope that these insights into the mechanism and conditions for the synchronization of charge oscillations will be helpful for future experiments related to ultrafast control of charge motion.
Within a mean-field approximation, charge-oscillation dynamics is described by the Bloch equations for pseudospins in a similar manner to those for photoinduced magnetization 59) and pairing 60) dynamics, or photoinduced excitonic condensation dynamics, 25) which allows us to intuitively grasp the essence of the collective evolution. The time-dependent BCS pairing problem is known to be integrable, 60) and its similarity to and difference from the present problem needs to be clarified. On short timescales, the mean-field picture basically holds even when electron correlations are taken into account. 31,59) The numerical calculations in this paper consider onedimensional dimerized lattices. They are bipartite and consist of two sublattices. However, the synchronization mechanism described by the Bloch equations is independent of the dimensionality of the system, and the synchronization phenomena have been numerically observed in one-and two-dimensional systems. [53][54][55] As representative systems, we treat the Hubbard and a spinless fermion model. The former represents a model with repulsive interactions within a sublattice, while the latter possesses repulsive interactions between the sublattices, and their synchronization conditions are different. To show the sensitivity to the initial condition in the latter model, we prepare initial nonequilibrium states in different ways, by photoexcitation and quenching.

Pseudospin Representation
To describe collective dynamics, the pseudospin representation is often useful. Here, following Ref. 39, we employ a mean-field approximation to treat charge-oscillation dynamics in this representation and investigate the mechanism for synchronization. We consider a Hubbard model with onsite repulsion U and a spinless fermion model with nearestneighbor repulsion V on dimer lattices. Although the discussion in this section does not depend on the dimensionality of the system, we use one-dimensional lattices with alternating transfer integrals, t 1 and t 2 , and periodic boundary conditions when the kinetic term is explicitly shown, to make the relationship with the results in the next section clear. The distance between neighboring sites is set to be equal and unity when the system is photoexcited. 53) The Hubbard model can be written as where e † jσ (e † j ) creates an electron with spin σ (a fermion) at site j of the "even" sublattice and o † jσ (o † j ) creates an electron with spin σ (a fermion) at site j of the "odd" sublattice.
In what follows in this section, we mainly discuss the Hubbard model and its charge-oscillation dynamics. When nontrivial differences between the two models exist, they will be pointed out. We use a mean-field approximation corresponding to the Hartree approximation, which allows us to capture the essence of interaction effects. Via the Fourier transforms e kσ = (1/ √ N) j e −ik j e jσ and o kσ = (1/ √ N) j e −ik j o jσ with N denoting the number of dimers, the mean-field Hamiltonian can be written as with for the Hubbard model and as with for the spinless fermion model, where constant terms have been omitted. On the one-dimensional lattice, the off-diagonal element is given by Following Ref. 39, we define and rewrite the meanfield Hamiltonian in Eq. (3) as with h ′ (k) = sgn Re (h(0)) |h(k)| and δn in Eq. (4) as As in Ref. 60, we define the pseudospin components as Here the normalization condition r 2 1kσ + r 2 2kσ + r 2 3kσ = 1 is always satisfied for kσ for which one state is initially occupied and the other is initially unoccupied, and r 2 1kσ + r 2 2kσ + r 2 3kσ = 0 is satisfied otherwise. In equilibrium with δn = 0 and h ′ (k) < 0 (i.e., t 1 + t 2 < 0), r 3kσ = −1 for kσ in the former case. Further, we define Ω 60) as Ω = (U/2)δn. The self-consistency condition leads to The pseudospin dynamics is described by the Bloch equations,ṙ which are consistent with the normalization condition. For For the spinless fermion model with Ω ≡ 2Vδn = (V/N) k r 1k , the Bloch equations can be obtained by replacing Ω by −Ω in Eqs. (15)- (17). For V larger than a critical value V c , there is a static solution with Ω 0, and the dynamical problem becomes equivalent to the time-dependent BCS pairing problem. 60) Equations (15)-(17) describe the Larmor precession of the pseudospinṙ in a pseudo magnetic field where the time dependence is made explicit here, and B 3k = −2h ′ (k) > 0 for h ′ (k) < 0. The x, y, and z components of the pseudospin correspond to the charge-density difference, current density, and bond density between the "even" and "odd" sublattices. The rate of rotation around the z axis depends on the wave number. In what follows, we assume that U > 0. The rotation axis is tilted by Ω(t) whose sign is that of δn(t), i.e., that of the majority of r 1kσ (t). When the majority of r 1kσ (t) take positive (negative) values, Ω(t) is positive (negative), and the rotation axis is tilted to have a positive (negative) x component, as shown in the left (right) panel of Fig. 1. By tilting the rotation axis in this way, most of the pseudospins are located on the same side of the rotation axis, those near the bottom t (0, 0, −1) are accelerated, and consequently their dynamics become similar, which assists the synchronization of charge oscillations. These considerations suggest that charge oscillations are synchronized regardless of the initial distribution of r 1kσ . On the right hand side of Eq. (16), the sign of h ′ (k) and that of r 3kσ are the same, as is obvious from Eqs. (9) and (13), at least near equilibrium and the sign of the majority of r 1kσ (t) and the sign of Ω(t) are the same, as already pointed out. Thus, a U-driven (see Eq. (14)) but k-independent force component is applied in the direction of the majority of the k-dependent current flow, which is equivalent to the above picture based on the tilted rotation axis.
A solution to the Bloch equations can be obtained, according to Ref. 60, by the ansatz The Bloch equations and the normalization condition lead to If the above equation is independent of k, and of the forṁ it has a solution described by the elliptic cosine function, When we set the amplitude ∆ + = δ and the frequency 53) the elliptic modulus becomes ∆ + / ∆ 2 + + ∆ 2 − = δ/(−2h ′ (0)). In order for Eq. (21) to become Eq. (22), the k-dependent factors must satisfy The negative sign in Eq. (24) is due to the self-consistency condition, Eq. (14), which can be rewritten as 1/U = 1/(4N) ′ kσ 2h ′ (k)C k , where the summation is over the kσ for which the pseudospin has a nonzero magnitude. In the limit of vanishing oscillation amplitude δ → 0, D k → 1, and r 3kσ → −1, i.e., the solution indeed approaches the equilibrium state. Although the initial condition must satisfy Eqs. (24) and (25) for the pseudospin dynamics to be exactly described by the elliptic cosine function and this is not generally achieved by a photoexcitation, this solution strongly suggests that charge oscillations are synchronized by a sufficiently large U, which has indeed been observed previously 53,54) and will also be demonstrated in the next section.
For the spinless fermion model, the situation is quite different from the Hubbard model. For V > V c , a static charge ordered state is realized in equilibrium. As mentioned before, its nonequilibrium dynamics is known 60) to be described by another elliptic function, the delta amplitude, whose smallamplitude limit corresponds to the Higgs amplitude mode. For 0 < V < V c , the excitonic effect stems from the coupling among electron-hole pairs a † k b k with different k's through the element 2Vδn of the mean-field Hamiltonian. An electron and a hole are bound by V to form an exciton, which is described by a linear combination of a † k b k . If the initial state contains such excitons, the charge-oscillation dynamics contains a slow component and the synchronization is difficult to achieve. In the pseudospin representation, the rotation axis is tilted in the direction opposite to the Hubbard case. Thus, a V-driven k-independent force component is applied in the direction opposite to the current flow, which allows the existence of a charge-ordered state for V > V c . For 0 < V < V c , however, this k-independent force component may lead to a synchronization if the initial state does not have a significant excitonic component. This fact is numerically demonstrated in the next section.

Dynamics from Charge Disproportionation
Employing one-dimensional lattices, Eqs. (1) and (2), with periodic boundary conditions, we numerically solve the timedependent Schrödinger equation to investigate the chargeoscillation dynamics from far-from-equilibrium states using the exact diagonalization method. 61) For the transfer integrals, we use t 1 =−0.3 and t 2 =−0.1, from which the frequency of the electronic breathing mode is given by 53) ω osc = 2(|t 1 | + |t 2 |) = 0.8.
The time slice employed in the numerical solutions of the time-dependent Schrödinger equation is dt=0.02 for small U or V, while smaller dt values are used for large U or V to ensure the conservation of the total energy and of the norm. The number of sites is denoted by L (L = 2N). We use system sizes up to L=16 for the Hubbard model at quarter filling and up to L=24 for the spinless fermion model at half filling. When photoexcitations are used to produce nonequilibrium states, the procedure is the same as in previous studies, 53,54) i.e., we use symmetric one-cycle electric-field pulses, A(t) = (cF/ω c ) [cos(ω c t) − 1] θ(t + (2π/ω c ))θ(−t), with central frequency ω c and field amplitude F. In the quench calculations, we add a staggered potential −∆ jσ e † jσ e jσ − o † jσ o jσ to Eq. (1) or −∆ j e † j e j − o † j o j to Eq. (2) to obtain the ground state and use this state as the initial state at t=0. The charge-oscillation dynamics is then computed for t > 0 with ∆=0. Fourier spectra of the charge density denote the absolute values of the Fourier transforms of the time profiles (0 < t < 10 3 ) of the charge density immediately after setting ∆=0 or immediately after the photoexcitation.
Before showing numerical results, we consider the Hubbard model in the limit of infinitely large repulsion, U = ∞. In this limit, it becomes equivalent to a spinless fermion model at doubled filling without interaction, with h ′ (k) = − t 2 1 + t 2 2 + 2t 1 t 2 cos k, i.e., to an integrable system; thus, independent of the frequency, the charge oscillations do not decay. At quarter filling, (N/2) spin-up and (N/2) spin-down electrons are present. When an electron is virtually moved through a distance of L sites, it passes over (N/2) fermions with the opposite spin. If (N/2) is odd, the outcome is different from the similar process in the noninteracting spinless fermion model or in the fully polarized Hubbard model by a factor (−1). Thus, the allowed wave numbers are k = (2π/N) j with j being an integer for even (N/2) and k = (π/N)(2 j + 1) for odd (N/2). On the other hand, the Hubbard model with U=0 is equivalent to the noninteracting spinless fermion model at the same filling; thus, the allowed wave numbers are always k = (2π/N) j. The charge-density difference 2δn can be rewritten as Thus, it is maximized by the state k , which is set to be the initial state of |Ψ(t) , so that In this state, the time evolution of 2δn is given by Its Fourier spectrum has peaks at ω k = 2 t 2 1 + t 2 2 + 2t 1 t 2 cos k with the k values described above, i.e., k=0, π for U=0 and ±π/2 for U=∞ if L=4 (N=2), k=0, ±π/2, π for U=0 and U=∞ if L=8 (N=4), etc.

Hubbard model at quarter filling
First we show in Fig. 2 that the Fourier spectrum of the charge density, σ e † jσ e jσ or σ o † jσ o jσ , evolving from a near-equilibrium initial state, has peaks at energies where the optical conductivity spectrum has peaks. Here, the optical conductivity spectrum is calculated for the ground state as in previous studies. 49) It is also shown in this figure that the Fourier spectrum of the charge density evolving from the near-equilibrium state with small charge disproportionation (∆=10 −2 ) is quite different from that evolving from a state with large charge disproportionation (∆=1). The latter is dominated by the electronic breathing mode at ω=0.8, although this mode is significantly broadened, reflecting its short lifetime.
The Fourier spectra of the charge density evolving from the state with large charge disproportionation (∆=1) are shown in Figs. 3(a) to 3(c) for different system sizes from L=4 to L=16. For very weak repulsion (U=10 −2 ), the spectra have discrete peaks whose energies are given by ω k with k = (2π/N) j. For very strong repulsion (U=8.0), the systems become close to the noninteracting spinless fermion model with periodic (for even N/2) or antiperiodic (for odd N/2) boundary conditions; thus, the spectra have discrete peaks whose energies are given by ω k with k = (2π/N) j or k = (π/N)(2 j + 1), respectively. In both cases, the systems are close to integrable ones, and all the charge oscillations are long-lived, i.e., the peaks are narrow and high. In the figures for larger systems, the corresponding peak heights are multiplied by factors smaller than unity, which allows them to be compared with intermediate repulsion (U=0.8) cases. Correlation effects become significant when the repulsion strength is intermediate. In this case, charge oscillations are synchronized 53,54) until they finally decay due to dephasing. The resultant electronic breathing mode has a broad peak at ω=ω osc . As the system size L increases, the difference between the very weak/strong and the intermediate repulsion cases, i.e., the difference between the nearly integrable and the strongly correlated cases, becomes apparent in Fig. 3. For the largest system (L=16) with intermediate repulsion (U=0.8) we have calculated the Fourier spectra with different initial states, as shown in Fig. 4. The initial states with large charge disproportionation are obtained either as the ground state with large alternating site energies (∆=1) or by the application of one-cycle large-amplitude (F=0.5) electric-field pulses with different central frequencies ω c . Although the detailed structures are different, all the spectra are similar and dominated by the electronic breathing mode at ω=0.8. Thus, the synchronization of charge oscillations is independent of the preparation of the initial state. It is universally observed  for states with large charge disproportionation in the Hubbard model with intermediate repulsion.

Spinless fermion model at half filling
In the previous section, we mentioned that the excitonic effect should be taken into account in the model with nearestneighbor repulsion. This implies that the charge-oscillation dynamics is sensitive to how the initial state is prepared. Hence, we mainly use one-cycle electric-field pulses with a central frequency near the absorption edge at ω=0.8, ω c =0.7, which allow excitations of many free electron-hole pairs. We  show in Fig. 5 that the Fourier spectrum of the charge density, e † j e j or o † j o j , evolving from a near-equilibrium state, has peaks at energies where the optical conductivity spectrum has peaks. Both spectra are dominated by an excitonic peak at ω ≃0.4 even if the near-equilibrium initial state is prepared by a weak photoexcitation (F=10 −2 ) with ω c =0.7. This figure also shows that, if the initial state is prepared by a strong photoexcitation (F=0.5) to be far from equilibrium, the Fourier spectrum is dominated by the electronic breathing mode at ω=0.8, which is in contrast to the Fourier spectrum after a weak photoexcitation.
The Fourier spectra of the charge density after strong photoexcitations (F=0.5) are shown in Fig. 6 for different repulsion strengths V below the critical value V c . For V=0, the system is integrable, and all the charge oscillations have infinitely long lifetimes. The corresponding discrete peaks appear at ω k with k=(π/6) j, and their peak heights are multiplied by 0.1 in order for them to be comparable with the nonintegrable cases. For small V (V=0.1), some structures in the Fourier spectrum reflect the peaks at V=0, but they are significantly  broadened. For V=0.1, 0.2, and 0.3, the Fourier spectra are dominated by the electronic breathing mode at ω=0.8, but its peak height decreases with increasing V. For V=0.4, the corresponding peak appears slightly blueshifted, although V is still smaller than V c . In general, the peak associated with the electronic breathing mode is considerably broadened when V approaches V c ≃ 2|t 1 |=0.6, and it disappears in the chargeordered phase at V > V c . 53) For V=0.2, we have also calculated the Fourier spectra with different initial states with large charge disproportionation, as shown in Fig. 7. They are sensitive to the initial condition. When the initial state is prepared by a photoexcitation with ω c =0.7, the spectrum is dominated by the electronic breathing mode. However, when ω c is lowered to allow the exciton at ω ≃0.4 (Fig. 5) to be excited, the Fourier spectra have large contributions from the exciton and smaller contributions from the electronic breathing oscillation. When the initial sate is prepared by ∆=1, the Fourier spectrum becomes more complex, presumably because various electron-hole pairs are involved. Thus, to synchronize charge oscillations, the initial state must be strongly photoexcited with a frequency that produces many free and few bound electron-hole pairs.

Conclusions
To clarify the mechanism for the previously observed synchronization of charge oscillations [52][53][54] and investigate the conditions for the appearance of this phenomenon, we analytically studied the Bloch equations in the mean-field approximation and numerically studied the charge-oscillation dynamics using the exact diagonalization method for the Hubbard and spinless fermion models on one-dimensional dimerized lattices.
In one of the Bloch equations, the time derivatives of current densities between the sublattices are determined by two terms, one of which is of kinetic origin and the other is of interaction origin. The kinetic term depends on the wave number, while the term derived from the interaction is proportional to a wave number independent and thus universal factor. The latter facilitates the synchronization of charge oscillations. In the Hubbard model with U > 0, more precisely speaking, when the interaction within a sublattice is repulsive, the kinetic term and the interaction term have the same sign and constructively work together. From the viewpoint of pseudospins, the rotation axis describing the Larmor precession is alternatingly tilted in such a way that it leads pseudospins of different wave numbers to be synchronized. In the spinless fermion model with V > 0, more precisely speaking, when the interaction between the two sublattices is repulsive, the kinetic term and the interaction term have opposite signs. This allows the existence of static charge order for V > V c . For V < V c , the synchronization of charge oscillations is possible but requires suitable excitation protocols.
The above insights based on the analytic form of the Bloch equations are all consistent with the Fourier analyses of the numerically obtained charge-oscillation dynamics. To prepare initial states with charge disproportionation, we either quenched a staggered potential to zero or applied one-cycle electric-field pulses. In the Hubbard model with intermediate repulsion, i.e., sufficiently away from the integrable limits of U=0 and U=∞, the synchronization is achieved irrespective of how the initial state is prepared. In the spinless fermion model with V < V c , a sufficiently large number of free electron-hole pairs must be excited to achieve the synchronization, which is hindered by excitons and is not achieved by the quenching. These results clarify the conditions for the synchronization of charge oscillations on dimer lattices, which may guide the study of synchronization phenomena in different classes of lattices.