Evolution of Charge-Lattice Dynamics across the Kuramoto Synchronization Phase Diagram of Quantum Tunneling Polarons in Cuprate Superconductors

Abstract

Because of its sensitivity to the instantaneous structure factor, S(Q,t = 0), Extended X-ray Absorption Fine Structure (EXAFS) is a powerful tool for probing the dynamic structure of condensed matter systems in which the charge and lattice dynamics are coupled. When applied to hole-doped cuprate superconductors, EXAFS has revealed the presence of internal quantum tunneling polarons (IQTPs). An IQTP arises in EXAFS as a two-site distribution for certain Cu–O pairs, which is also duplicated in inelastic scattering but not observed in standard diffraction measurements. The Cu–Sr pair distribution has been found to be highly anharmonic and strongly correlated to both the IQTPs and to superconductivity, as, for example, in ${\mathrm{YSr}}_2{\mathrm{Cu}}_{2.75}{\mathrm{Mo}}_{0.25}{\mathrm{O}}_{7.54}$$(T_c=84$ K). In order to describe such nontrivial, anharmonic charge-lattice dynamics, we have proposed a model Hamiltonian for a prototype six-atom cluster, in which two Cu-apical-O IQTPs are charge-transfer bridged through Cu atoms by an O atom in the ${\mathrm{CuO}}_2$ plane and are anharmonically coupled via a Sr atom. By applying an exact diagonalization procedure to this cluster, we have verified that our model indeed produces an intricate interplay between charge and lattice dynamics. Then, by using the Kuramoto model for the synchronization of coupled quantum oscillators, we have found a first-order phase transition for the IQTPs into a synchronized, phase-locked phase. Most importantly, we have shown that this transition results specifically from the anharmonicity. Finally, we have provided a phase diagram showing the onset of the phase-locking of IQTPs as a function of the charge-lattice and anharmonic couplings in our model. We have found that the charge, initially confined to the apical oxygens, is partially pumped into the ${\mathrm{CuO}}_2$ plane in the synchronized phase, which suggests a possible connection between the synchronized dynamic structure and high-temperature superconductivity (HTSC) in doped cuprates.

1. Introduction

The absence of an accepted theory for superconductivity in cuprates has recently resulted in calls for new ideas [1,2]. One possibility that was observed soon after their initial discovery but was subsequently neglected is their dynamic local lattice structure, which, for these mixed valence oxides, will be coupled to their charge distribution [3]. The results from X-ray absorption fine structure (XAFS), inelastic neutron scattering (INS) and inelastic X-ray scattering (IXS), Raman spectroscopy (RS), infrared absorption spectroscopy (IR), time-resolved optical spectroscopy and nonlinear magnetic response measurements [4,5,6,7,8,9,10,11,12,13] have pointed to an inhomogeneous ground state in different hole-doped cuprates as a precursor to the superconducting phase. Anomalous isotopic effects [14] and short-range structural fluctuations [15] are among the manifestations of this inhomogeneity. In particular, inelastic neutron scattering [16,17] and X-ray absorption spectroscopy [18], which are sensitive to the dynamic structure [19,20], revealed two-site distributions of certain Cu–apical O pairs [21,22,23,24]. Since the static structure determined by crystallography only identified single positions for these atoms [25], a proposed interpretation of these anomalous pairs was that they were the result of the oxygen atoms and their associated charges undergoing quantum tunneling between the two sites [3,19,20,26]. In addition, the difference between the two Cu–O distances indicated that these sites were associated with different charges on the oxygen. We have recently described this phenomenon as internal quantum tunneling polarons (IQTPs) [27], with the “internal” referring to the internal tunneling frequency being much larger than the hopping rate that is a characteristic of small polarons. On this time scale, the lattice distortion and charge inhomogeneity are therefore substantially confined within the target pair of atoms, the small polaron core. The internal dynamics of the small, substantially immobile polaron can also be viewed as lattice-assisted dynamic charge-transfer.

The potential relevance of IQTPs to superconductivity follows from the original finding that the separation between the two sites and associated tunneling frequency are observed to be modified through the superconducting transition [11,18,24,26]. Subsequent measurements on compounds in which a fraction of the Cu was substituted by other metals were interpreted as showing that the IQTPs were, at a minimum, coupled to the superconductivity [28]. Although early neutron scattering experiments were essential in corroborating the IQTPs in the dynamic structure [16], supplemented by anomalies in vibrational spectra, the preponderance of the experimental results have been performed by XAFS. Its sensitivity to the local dynamic structure via its measurement of the instantaneous structure factor in samples of only a few milligrams, its precision with respect to both the metrical parameters of the local structure and the actual pair distributions and its ability to separate specific atom pairs via its element dependence supplemented by, e.g., magnetic orientation of the samples, give it an unrivaled capability to probe IQTPs and their coupling to the superconductivity.

Enhanced by improvements in the technique that allow the measurement of accurate spectra five times faster than originally possible, recently renewed XAFS experiments on cuprates have been performed, specifically highly overdoped superconducting compounds prepared by high-pressure oxygen methods [29,30,31,32]. These materials have provided surprising revelations on the importance of IQTPs. They have also provided counterexamples to the assumptions that higher transition temperatures are correlated with longer Cu–apical O distances, of the requirement for fully ordered ${\mathrm{CuO}}_2$ planes [33,34], of the prolate Cu geometry and resulting ordering of the Cu 3D states [35] and have even identified a very large structural transformation concomitant with the superconducting transition [27]. The remaining common property of superconductivity in hole-doped cuprates is therefore the IQTPs, as revealed by XAFS. Notably, these have also been observed in other exotic superconductors [36,37,38]. While these data do not establish a mechanism for superconductivity, they certainly demonstrate a coupling of the superconductivity and IQTPs. We have therefore renewed efforts to understand the properties and consequences of IQTPs.

Stimulated by these new results, we have resumed calculations exploring the behavior of the dynamic structure and its coupling to other properties. The original calculations were performed on a three-atom, O–Cu–O cluster, limited to examining the effects of the dynamic structure on the structure factors and the electronic states [19]. We have now enlarged the cluster to include an oxygen atom in the ${\mathrm{CuO}}_2$ superconducting plane between two Cu–Oap pairs (IQTP environment) and the anharmonically coupled Sr atom that bridges the apical oxygen atoms (they are bridged by Ba in the majority of compounds) [39]. This final addition is the result of an early EXAFS finding of strong anharmonicity at the divalent alkaline earth site inducing major structural distortions in neighboring oxygen atoms [40], corroborated by our results and expanded to show that this anharmonicity is also coupled to the superconductivity [27,34,41]. The minimal cluster incorporating these new effects is the six-atom one depicted in Figure 1. This generic moiety occurs in various hole-doped cuprates. Although we use the specific structure from ${\mathrm{YSr}}_2{\mathrm{Cu}}_3{\mathrm{O}}_{7.54}$ [32], our objective is not to duplicate or even examine the XAFS results and their connection to the superconductivity. At this stage of our research, our goal is to elucidate whether this larger cluster that incorporates additional elements of the actual structure exhibits any properties or behaviors beyond those of the static crystallographic structure that stimulate future exploration of the coupling of the IQTPs to anharmonicity and to the superconductivity.

After rigorous numerical analysis, a phase diagram is obtained in terms of the electron–phonon coupling $\lambda$ at the apical oxygen sites and the anharmonic phonon coupling K brought by the presence of a triatomic Oap–S-r-Oap molecular structure (see Figure 1). In view of the Kuramoto model for sycnhronization of quantum oscillators [42], a synchronized phase for the IQTPs is found, and the evolution of charge-lattice dynamics is studied across the phase diagram. Finally, we introduce an unitary transformation that gives rise to a novel electron-lattice coupling in the central, planar site of the cluster that is only finite in the synchronized IQTP regime, proportional to the Kuramoto order parameter, which then connects the numerical results with the theoretical analysis based on the synchronization of the IQTPs.

2. Model and Methods

Soon after the discovery of high-temperature superconductivity in doped cuprates, the initial observations of IQTPs found their signature two-site distributions in the apical oxygens (Oap) [5,18,21,23,43,44]. Remarkably, however, the same type of two-site signature of IQTPs was also found in some other compounds, but for the O atoms in the ${\mathrm{CuO}}_2$ planes (Opl) [6,12,22,45,46], and even in non-copper-based HTSC materials [36,37,38]. The application of exact diagonalization procedure for a minimal three atom cluster, based on an O–Cu–O chain, captured the minimal physics of an IQTP and predicted intrinsic local physical properties and at the same time validated this numerical technique for incorporating quantum dynamics [3,19,26,43,44]. More recent calculations examining fluctuating stripes within the ${\mathrm{CuO}}_2$ planes [47] and the role of the cation on the opposite side of the Oap from the Cu [48] also included dynamical contributions.

Although these reports represent an advance in the treatment of the lattice, they do not address the anharmonicity and its manifestation in the double well potential of the Oap atoms of the IQTPs. By expanding to our six-atom cluster (see Figure 1), we are able to study the effects of a triatomic molecule, namely, the chain Oap–A–Oap, where A is an atom in the dielectric layer on the IQTPs and in the ${\mathrm{CuO}}_2$ plane. Here, we describe the atom in the inert layer as the Sr atom, as in Figure 1, based on the ${\mathrm{YSr}}_2{\mathrm{Cu}}_{2.75}{\mathrm{Mo}}_{0.25}{\mathrm{O}}_{7.54}$ structure.

2.1. The Six-Atom Cluster

In order to capture the intriguing complexity of charge and lattice dynamics in the reduced cluster, we introduce a Hamiltonian that captures the interplay between the charge, harmonic phonons and anharmonically coupled molecular phonons resulting from the presence of a triatomic molecular structure, Oap–Sr–Oap, within the cluster. The construction of this Hamiltonian is based on an extension of the previous three-atom cluster but now with the novelty of the anharmonicity term. Namely, we write $H=H_{el}+H_{ph}+H_{el-ph}+H_{anh}$, where:

$$\begin{array}{ccc}\hfill H_{el}& =& \sum _i{\epsilon }_i{n}_i+t\sum _{\langle ij\rangle ,\sigma }(c_{i\sigma }^{†}{c}_{j\sigma }+h.c.)+U\sum _i{n}_{i\uparrow }{n}_{i\downarrow }\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill H_{ph}& =& \hslash \omega b_L^{†}{b}_L+\hslash \omega b_R^{†}{b}_R+\hslash \Omega {\beta }^{†}\beta \hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill H_{el-ph}& =& \lambda n_L\left(b_L^{†}+b_L\right)+\lambda n_R\left(b_R^{†}+b_R\right)\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill H_{anh}& =& K\left({\beta }^{†}{b}_L{b}_R+\beta b_L^{†}{b}_R^{†}\right)\hfill \end{array}$$

(4)

In the first term (1), $c_{i\sigma }^{†}\left(c_{i\sigma }\right)$ are the usual creation (annihilation) fermionic operators, for charges with spin projection $\sigma =\uparrow ,\downarrow$ at sites $i=1,\dots ,5$ and $n_i={\sum }_{\sigma }{c}_{i\sigma }^{†}{c}_{i\sigma }$. It describes the electronic contributions to the cluster, via on-site energies ${\epsilon }_i$ at each of the five charge transfer sites i, namely, the Oap–Cu–Opl–Cu–Oap chain, where we emphasize that the Sr site is purely structural, giving the cluster its six-atom characteristic. t is the spin-preserving, nearest-neighbor hopping amplitude, which will be our standard energy scale throughout this work, and U is the amplitude of the on-site Hubbard repulsion term, preventing double occupancy ($U/t\gg 1$). The bosonic term (2) describes the harmonic contribution of the three phonon modes in the cluster: $b_L^{†}\left(b_L\right)/b_R^{†}\left(b_R\right)$ creates (annihilates) the infrared left/right ($i=1$ and $i=5$) apical-oxygen-related phonon modes and ${\beta }^{†}\left(\beta \right)$ creates (annihilates) the triatomic molecular vibration associated with the Oap–Sr–Oap chain. In other words, it locks the vibration of the two isolated apical oxygens into a composite molecular mode of vibration. Equation (3) describes the coupling between the electronic degrees of freedom to the lattice deformations at the apical oxygen sites in the form of an electron–phonon coupling $\lambda$, and $n_{L,R}$ are the usual electronic occupations at these sites ($i=1$ and $i=5$). Finally, the last and novel term in Equation (4) provides the cluster with a nontrivial Oap–Sr–Oap triatomic molecule structure, motivated by the coupling of the anharmonicity brought by the Sr ion to the behavior of the Oap deformations. This last term is anharmonic in the sense that it describes a three-phonon interaction, controlled by the anharmonic coupling K, and it can be interpreted as a chemical potential that control processes where left/right independent apical oxygen vibrations are destroyed to created a composite, triatomic molecular vibration, ${\beta }^{†}{b}_L{b}_R$, as well as its conjugate, when the molecular vibrations decays into two independent vibrations of left/right phonons in the cluster. As we shall demonstrate below, the anharmonic amplitude K promotes the locking of phases at each apical oxygen, or equivalently, the synchronization of Kuramoto oscillators [49]. For this reason, from now on, K will be referred to as the Kuramoto coupling of the cluster.

The exact diagonalization applied to the full Hamiltonian, Equations (1)–(4), was perfomed using the basis wave functions written as:

$$\begin{array}{c}\hfill |\Psi \rangle =\sum _{k\ell rs}{\alpha }_{k\ell rs}|n_k\rangle |n_{\ell }\rangle |n_r\rangle |n_s\rangle ,\end{array}$$

(5)

in an $n_{el}\times n_L\times n_R\times n_{\beta }$ dimensional Hilbert space. For the number of electronic states $n_{el}$, we consider a total of three holes added to the five charge transfer sites and work out all possible spin-up and spin-down configurations, for each $|n_k\rangle$. However, only one extra hole is allowed to hop around the cluster while the other two are fixed at the Cu sites, giving one additional hole per planar copper–oxygen unit [50]. This results in $n_{el}=28$ electronic states. Furthermore, $n_L$, $n_R$ and $n_{\beta }$ are the number of phonon modes related to the left and right apical oxygen sites, as well as to the triatomic molecular phonon, respectively, and the bosonic states $|n_{\ell }\rangle |n_r\rangle |n_s\rangle$ are written in an occupation number representation. By continuously enlarging the phonon Hilbert space, through a systematic increase in the number of phononic modes, we have found convergence of the exact diagonalization procedure with 5 modes for each type of boson in the cluster, totalizing a 6048-dimensional electron–phonon Hilbert space. For the numerical calculation, derived from previous works [26], we used the following values: ${\epsilon }_i=0.5$ eV $(i=1,3,5)$ at the oxygen sites and ${\epsilon }_j=-0.5$ eV $(j=2,4)$ in order to favor charge localization at the copper positions, $t=0.5$ eV as our standard energy scale, $U=7.0$ eV as the onsite Coulomb repulsion, $\hslash \omega =0.06$ eV for the infrared phonon mode of the apical oxygens and, finally, $\hslash \Omega =0.12$ eV $=2\hslash \omega$, constrained by energy conservation, as a collective molecular mode for the triatomic chain and built out of two independent oscillators.

2.2. The Kuramoto Model

The spontaneous emergence of collective behavior in large networks of coupled oscillators is a phenomenon that appears in many different branches of science, and the model that has become the simplest paradigm for describing the synchronization phenomena is the Kuramoto model [42]. For example, in the classical picture of synchronization, it is used to describe the alpha rhythms in human brains [51] and the locking of vibrations for metronomes sharing a commom base [52], just to name a few applications. As for its quantum analogue, the dynamics of fast spins coupled to slow exchange interactions in XY spin glasses [53], the frequency locking in superconducting Josephson junctions arrays [54], the emergence of bulk $d-$wave high-temperature superconductivity in inhomogeneous cuprates [55] and its manifestation in momentum space localization for ultracold atoms in a lattice [56] are examples of how synchronization, described by the Kuramoto model, is an important feature to be included in diverse scientific models.

The Kuramoto model relies on two basic properties: the couplings between the node oscillators, a superfluid density $K_{ij}$, and the existence of white, quenched, $\delta$ or thermal, $k_{B}T$, noises provided by the network embedding environment. Partial or full synchronization can be achieved whenever the couplings dominate the noise. This process is described by a complex order parameter, $r{e}^{i\psi }$, where the real part $r=0$ for the unsynchronized phase, while $0<r<1$ for partial and $r=1$ for full synchronization. In order to show that our cluster can be described by this model, with the anharmonicity parameter K playing the role of the coupling between oscillators, we first rescale $K\to K/N$, where N is the total number of mode oscillators, namely phonons, at each independent oscillator community at the left or right in the cluster. Second, since all oscillators are phonons, we introduce a mean field approximation to treat the anharmonicity term (4), in which the triatomic molecular phonon population can be written as $\langle {\beta }^{†}\beta \rangle ={\left|R\right|}^2=\langle {\beta }^{†}\rangle \langle \beta \rangle =R^{\ast }R\ne 0$, where R is a complex order parameter. Thus, the Heisenberg equations of motion for $i=L,R$ phonons derived from the full Hamiltonian read as follows:

$$\begin{array}{c}\hfill i\frac{d{b}_i}{dt}=\omega \left(b_i-\frac{\lambda n_i}{\hslash \omega }\right)+\frac{KR}{N\hslash }{b}_j^{†},\end{array}$$

(6)

together with its Hermitian conjugate, for $b_i^{†}$. Taking the quantum mechanical expectation value on each side of the last equation and introducing phonon coherent states, $b_{i,j}|z_{i,j}\rangle =z_{i,j}|z_{i,j}\rangle$, we can drop constant terms that are simply overall shifts for the equilibrium position of each oscillator. We write the complex number, $z_{i,j}=\left|z_0\right|e^{-i{\theta }_{i,j}\left(t\right)}$, where we make the reasonable assumption that the amplitude for vibrations at the left L and right R apical oxygen atoms, $z_0$, are equal and only their phases change. Since the claim is that the synchronization transition appears due to the anharmonic triatomic molecular vibrations, we introduce the Kuramoto order parameter as $R=r{e}^{i\psi }$. Then, motivated by our numerical results, we switch from the laboratory reference frame to the center of mass of the oscillator network, ${\theta }_i\left(t\right)+{\theta }_j\left(t\right)=0$, where all phases evolve opposite to one another in Kuramoto’s circle. In our calculations, this is equivalent to setting $\psi =-\pi /2$, or equivalently, $R=-ir$, which allows us to write:

$$\begin{array}{c}\hfill \frac{d{\theta }_i\left(t\right)}{dt}=\omega +\sum _{j=L,R}\frac{Kr}{N\hslash }sin\left[{\theta }_i\left(t\right)+{\theta }_j\left(t\right)\right].\end{array}$$

(7)

The above equation can be interpreted as a Kuramoto’s differential equation for anti-phase synchronization ${\theta }_i\to -{\theta }_j$. The mean field limit of a large number of oscillators at each left or right communities is then obtained by introducing the definition of the Kuramoto order parameter [42], $r{e}^{i\psi }=\frac{1}{N}{\sum }_{j=1}^N{e}^{i{\theta }_j}$ and by adding thermal noise into the problem, ${\xi }_i\left(t\right)$. The result is the following Kuramoto equation of motion for the phases of each community of oscillators:

$$\begin{array}{c}\hfill \dot{{\theta }_i}\left(t\right)=\omega +K{r}^{2}sin\left[{\theta }_i\left(t\right)+\psi \right]+{\xi }_i\left(t\right),\end{array}$$

(8)

with $\langle {\xi }_i\left(t\right)\rangle =0$ and $\langle {\xi }_i\left(t\right){\xi }_j\left(t^{\prime }\right)\rangle =2\gamma k_{B}T{\delta }_{ij}\delta (t-t^{\prime })$, where $\gamma$ is a damping constant associated with the normal modes of the environment. The anti-phase synchronization transition is obtained because, according to (8), phase locking occurs for ${\theta }_i\to -\psi$, and as such, the anti-phase synchronization favored by Equation (7) implies ${\theta }_j\to \psi$. Since we have chosen to work on the center of the mass reference frame, where $\psi =-\pi /2$ (such that $R=r{e}^{i\psi }=-ir$), we end up with ${\theta }_i+{\theta }_j=0$ and ${\theta }_i-{\theta }_j=-\pi$, which reflects the anti-phase synchronization transition. Once established that the phase locking is favored in anti-phase character caused by the anharmonicity, K, we can show that the nature of this synchronization transition is of the first order. Using consistency equations such as in Ref. [57], the solution to the first order differential Equation (8) yields for the Kuramoto order parameter:

$$\begin{array}{c}\hfill r=\sqrt{1-\frac{K_c(\delta ,T)}{Kr}},\end{array}$$

(9)

where $\delta$ is a quenched spread of frequencies due to disorder that enters in the distribution of frequencies of the system, as for example, a Lorentzian, $g\left(\omega \right)=\frac{\delta }{\pi ({\omega }^2+{\delta }^2)}$. Therefore, the first-order, anti-phase synchronization transition is obtained whenever the Kuramoto coupling of the system, K, exceeds a critical value, $K>K_c(\delta ,T)$, which is determined by disorder and temperature [57].

3. Results

The exact diagonalization procedure applied to the six-atom cluster was reported in Ref. [39]. In this work, we provide a deeper understanding of the regimes where small polarons are formed and when they are allowed to tunnel in the form of an IQTP as a function of the electron-phonon coupling, $\lambda$, and the Kuramoto–anharmonic coupling, K. We take different crossover couplings, ${\lambda }_c$, for the formation of polarons centered at the left/right apical oxygens for different values of the Kuramoto coupling, K, searching for a transition between adiabatic and non-adibatic regimes for polarons. We also search for different values of $K_c$ for the onset of IQTP synchronization within the cluster for different values of $\lambda$. These quantities are extracted by examining the evolution of the ground state energy $E_0$, as shown in Figure 2, bottom panel, where the critical couplings are determined by the minimum in the second derivative of $E_0$ with respect to each coupling, as shown in the insets. It is clear that the transition at ${\lambda }_c$ is defined by a crossover region, where polarons are being formed, thus corresponding to a second order phase transition. On the other hand, the transition at $K_c$ is a first-order phase transition, evidenced by a clear kink at $K_c$ in the ground state energy $E_0$ as it evolves with K. We have here numerical confirmation of the first order character of the synchronization transition predicted by the mean-field Kuramoto analysis.

After continuously evaluating all transitions, a set of different $(K,{\lambda }_c)$ and $(\lambda ,K_c)$ points were obtained and plotted in the $\lambda \times K$ plane to give the phase diagram as shown in Figure 2. This diagram displays four distinct phases: (I) a weak coupling, adiabatic regime, where no polarons or collective molecular vibrations are formed, $\lambda <{\lambda }_c$ and $K<K_c$. Here, charge is decoupled from all sorts of vibrations; (II) a polaronic, non-adiabatic regime ($\lambda >{\lambda }_c$), where holes are localized either at the left/right parts of the cluster (see Figure 3). Since K is not strong enough to ensure IQTP formation and since synchronization is not achieved ($K<K_c$), no quantum tunneling within the cluster is observed; (III) an IQTP and non-adiabatic regime, where synchronization is set up $(K>K_c)$ and holes are now coupled to the collective molecular vibration corresponding to the triatomic molecule, $O_{ap}-Sr-O_{ap}$, in such a way that quantum mechanical tunneling is now allowed, and the hole becomes delocalized. In this phase, a split-polaron is formed, and the cluster dynamic is of the polaronic charge plus collective molecular vibration type. Finally, (IV) back to an adiabatic molecular regime, where the electron–phonon coupling is not strong enough for polaron formation, although a collective molecular vibration is set up, since $K>K_c$. Here, the excess charge becomes decoupled from the lattice dynamics and is only weakly coupled to the collective molecular vibrations via a coupling strength of ${\lambda }_{eff}\equiv {\lambda }^{2}K$, where the apical oxygen vibrations are locked.

In order to understand how the electronic degrees of freedom behave at the different phases of the system, in Figure 3, we show the electronic occupations along the numbered path shown in the phase diagram, Figure 2, where we can see the evolution of the extra hole occupation at the left/center/right oxygen positions. It is clear that the synchronization transition (path $2\to 3$) drives the system to a delocalized state for the extra hole, pumping it to the copper oxide plane, together with the formation of IQTPs, since now tunneling is allowed (See Figure 4), in a clear evidence that the dynamic structure brought by the anharmonic Sr ion can influence the parameters related to the superconductivity in the system. Although a finite central occupation appears in the weak coupling regime, the hole is decoupled from the bosonic degrees of freedom, therefore, no dynamics associated to the displacement of the apical oxygens, nor the molecular vibration, can affect the hole delocalization. This is clear when we go through the path $3\to 4$, where the hole delocalization is insensitive to the the polaronic transition, since it is controlled solely by the molecular vibration ($K>K_c$), which pumps it to the planar site. When going back from phase $IV)$ to phase $I)$, we see a clear discontinuity at the transition, a consequence of the hole decoupling from the molecular vibration.

4. Discussion

The phase diagram in Figure 2, together with the evolution of the electronic occupations in Figure 3, show how the anharmonic dynamics associated with the motions of the apical oxygens that are the out-of-plane phenomena can influence the in-plane parameters, namely, the occupation in the central planar site. In the weak electron-lattice coupling regime, phase I), the hole is already delocalized, thus one may ask why it is important to look at the other regions of the diagram. For $\lambda <{\lambda }_c$ and $K<K_c$, the energy scale that dominates is the Hubbard term U; therefore, since no polarons are formed, neither molecular vibrations synchronizing the apical oxygens nor forming IQTPs, it is expected that the renormalization of the onsite energy, $\epsilon \approx {\epsilon }^0+\delta \epsilon$, associated with the central site gives a larger shift when compared to the lateral ones (left/right). This happens because, since the planar site has two neighboring Cu atoms, while the apical ones are only charge-bridged through one Cu, the energy renormalization at the central site is $\delta {\epsilon }_C\sim 2t^2/U$, while for the lateral ones it is $\delta {\epsilon }_{L,R}\sim t^2/U$. The factor of two for the planar site originates from the two distinct Cu neighbor atoms. Therefore, it is expected that for large Coulomb repulsion U, in the weak electron-lattice coupling regime, the Hubbard term dominates and saturates the central site electronic occupation at larger values than the lateral ones. However, since this electronic degree of freedom is decoupled from the rest of the cluster, no dynamic is associated with this phase. The numerical results are in complete agreement with this physical picture.

Within the Kuramoto analysis, in Figure 4, we show the solution to the Kuramoto’s order parameter Equation (9), where the first order phase transition is clear whenever $K/K_c>1$. In the derivation of Equation (8), we used the mean-field approximation in the anharmonic term of the full Hamiltonian and defined $R=r{e}^{i\psi }$ to be the order parameter. As a result, the mean field version of Kuramoto’s equation includes a $r^2$ term, which gives rise to a first order phase transition in the model. In Figure 4, we also show the evolution, across the path $2\to 3$, of the phonon occupation at the central site, related to the projection of molecular phonons in this site, $n_{\beta }\left(C\right)$, together with the tunneling frequency associated to the IQTP’s, ${\omega }_T=E_1-E_0$, which is the difference in energy between the first excited and ground states. Across the synchronization transition, the formation of IQTP’s is identified by the sudden emergence of a finite tunneling frequency, consistent with the first order character of the synchronization transition at $K_c$. Moreover, the transition project’s finite phonon occupation is at the central site of the cluster.

The pumping of charge onto the plane and the finite projection of molecular phonons at the central site suggest that, after the synchronization transition, there should be an induced, novel interaction between these two degrees of freedom. In order to obtain an effective interaction between the hole occupation at the planar oxygen, $n_3$, and the molecular vibrations described by $\beta$ and ${\beta }^{†}$, we perform a unitary transformation on the full Hamiltonian:

$$\begin{array}{c}\hfill H^{\prime }=e^{-S}H{e}^S,\end{array}$$

(10)

where the transformation matrix S is chosen to eliminate the linear terms associated with the electron–phonon interaction that is proportional to $\lambda$ and also the anharmonic K term in the transformed Hamiltonian, $H^{\prime }$. The matrix that satisfies this constraint is found to be:

$$\begin{array}{c}\hfill S=\frac{\lambda n_1}{\hslash \omega }(b_1-b_1^{†})+\frac{\lambda n_5}{\hslash \omega }(b_5-b_5^{†})+\frac{K}{2\hslash \omega -\hslash \Omega }\left({\beta }^{†}{b}_1{b}_5-\beta b_1^{†}{b}_5^{†}\right).\end{array}$$

(11)

After expanding the transformed Hamiltonian up to third order in the interactions and using the constraint of the electronic occupations across the cluster, $n_1+n_3+n_5=1$, we were able to write an effective interaction between the electron in the plane and the molecular vibration, which is nonzero only after the synchronization transition

$$\begin{array}{c}\hfill H^{\prime }\sim -{\lambda }^{2}K{n}_3(\beta +{\beta }^{†}).\end{array}$$

(12)

The contribution of this effective term is evident when we examine its diagramatic representation in terms of Feynman diagrams, as shown in Figure 4. There, the electrons, solid lines $\left(c\right)$, are coupled to the independent left/right apical oxygen deformations ($b_L$ and $b_R$) with amplitude $\lambda$ at the same time as these phonons (wavy lines) are coupled to the molecular vibration $\beta$ (curly line) with amplitude K. Therefore, the vertex renormalization gives rise to an effective interaction between the electrons and the triatomic molecular phonons with amplitude proportional to ${\lambda }^{2}K$. Moreover, we can argue that since this coupling is only finite after the synchronization transition, as suggested by the numerical calculations, the amplitude of this interaction should also be proportional to the Kuramoto order parameter. This results in an effective interaction of the order ${\lambda }^{2}Kr$, which is zero outside the IQTP phase. This is indeed the case once we recall that, during the calculation of the associated Feynman diagram, a sum over a product of left/right phonon phases is made in the loop, ${\sum }_{i=L,R}{e}^{i{\theta }_i}\equiv r{e}^{i\psi }$. This induced coupling suggests the formation of a novel, planar IQTP, allowing for the tunneling between the two lateral sites, and pointing towards the importance of the synchronization of lattice deformations in different parameters that might be coupled to the superconducting transition in cuprate superconductors.

5. Conclusions

In summary, although we restricted our analysis to a reduced cluster built to represent a much richer crystallographic environment of cuprate superconductors, in this work, we have included two novel features: (i) the planar $O_{pl}$ bridging the two $Cu-O_{ap}$ pairs and (ii) the structural anharmonicity related to the coupling within the dielectric layer, explored by the dynamical analysis based on the Kuramoto model. We have been able to describe the charge-lattice dynamics related to the motion of the apical oxygens in four different parameter regimes of the phase diagram and demonstrated how the lattice vibration couples to the excess charge in the ${\mathrm{CuO}}_2$ plane in the synchronized regime. We emphasize that this cluster is general in the sense that it only needs the above ingredients for the IQTPs and the anharmonicity from the dielectric layer of the compound, widening its possible application to other copper-based [11] as well as non-copper-based systems [37].

The obtained phase diagram points towards an important feature of the excess charge added to cuprates, since the localization/delocalization transition occurs through small variations in the interaction parameters, and depending on how these parameters are energetically related, one can have the electrons coupled to trapped small polarons or to global, synchronized molecular vibrations that give rise to the dynamics associated to IQTPs. Therefore, the addition of extra, out-of-plane dynamical properties to planar and strongly correlated models for cuprates is indeed able to modify the electronic and lattice characteristics of the conducting layer, unveiling its potential relevance to the superconducting phase in these systems.

Our results confirm the prominent coupling of nonadiabatic lattice dynamics and anharmonic oscillations of atoms and charge. As the next step, quantitative implications are delegated to future research. We anticipate that collective patterns of spatially correlated IQTPs (coupled by, e.g., Coulomb or strain fields) could play a role in the pairing of charges and their condensation and the synchronization of the apical oxygen oscillations across distances comparable to the superconducting coherence length. Future research will include coupling (via spin, charge or lattice fluctuations) of the IQTPs to the material medium in which they reside. We also note that the dynamical process of charge coupled to structural instabilities does not induce static deformations of the structure: This is a purely structural process that is dependent on external parameters, such as the ones introduced via the Kuramoto model. We emphasize that we do not have the microscopic proof that the synchronization transition we have studied in this work drives a superconducting one, although this is an intriguing possibility which we will continue to explore. However, our results do elucidate the importance of considering how lattice vibrations and their nonlinear extensions can play a role in high-temperature superconducting cuprates while providing us with an elegant example where a first-order synchronization transition between phonon vibrations, predicted by a mean-field Kuramoto analysis, is confirmed by exact diagonalization.