Cluster dynamical mean-field calculations for TiOCl

Based on a combination of cluster dynamical mean field theory (DMFT) and density functional calculations, we calculated the angle-integrated spectral density in the layered s=1/2 quantum magnet TiOCl. The agreement with recent photoemission and oxygen K-edge x-ray absorption spectroscopy experiments is found to be good. The improvement achieved with this calculation with respect to previous single-site DMFT calculations is an indication of the correlated nature and low-dimensionality of TiOCl.


Introduction
The low-dimensional quantum spin system TiOCl has received a lot of attention in recent years due to its anomalous behavior in a wide range of temperatures. It consists of Ti-O bilayers in the ab-plane separated by layers of Cl − ions stacked along the crystallographic c-axis (see figure 1). Ab initio density functional calculations [1]- [3] showed that the system at room temperature can be described in terms of spin-1/2 (Ti 3+ ) Heisenberg chains running along the crystallographic b-axis, along with small but non-negligible inter-chain couplings. Upon varying the temperature, TiOCl undergoes two successive phase transitions, one of second-order nature at T c2 = 91 K and one of first-order nature at T c1 = 67 K to a spin-Peierls dimerized state [1,4]. The nature of the transition at T c2 = 91 K has been discussed within various scenarios of orbital fluctuations [2,5] 5 , frustration of interchain interactions [7,8], the role of phonons [9]- [11], spin and correlation effects [12,13] and is still debatable. Even the behavior of TiOCl in the hightemperature phase has not been understood in a satisfactory way [14,15].
TiOCl can be viewed as a Mott insulator where the insulator gap is driven by electron correlation. Attempts to describe the influence of correlation on the properties of TiOCl have been performed by considering the local density approximation (LDA) + U approach [1]- [3] as well as the more elaborate single-site LDA + dynamical mean field theory (DMFT) approach [12,13] where dynamical fluctuations-absent within the LDA + U approachare considered. Nevertheless, comparison of the spectral function obtained from single-site LDA + DMFT calculations with photoemission spectroscopy (PES) experiments shows only a moderate agreement [14]. While the width of the computed spectral function is in accordance with the PES results, the magnitude of the optical correlation gap was highly underestimated. Also the agreement of the line shape is not very satisfactory.
The effective low-dimensionality of the material suggests that correlation-driven inter-site fluctuations may be important for the proper description of TiOCl. In fact, it was shown in a recent work [16] on a one-dimensional extended Hubbard model that consideration of two-site clusters within the cluster extension of DMFT (cluster-DMFT) leads to a satisfactory description of the charge gap. One may therefore expect for TiOCl-which is effectively a low-dimensional system-that an improvement upon previous single-site LDA + DMFT results may be achieved by including nonlocal correlation effects. Demonstration of such a property in a real material of such complexity as TiOCl has not been done previously. In order to analyze this proposal, we We use the LDA Wannier functions obtained by an N th order muffin-tin orbital-based (NMTO) downfolding method ( [17] and references therein) to construct a Hubbard Hamiltonian which we solve within the cluster-DMFT approximation.

Crystal structure and energy levels
The octahedral environment of Ti [TiO 4 Cl 2 ] in TiOCl splits the five degenerate d levels into t 2g and e g blocks. Since the Ti 3+ ion is in a 3d 1 configuration, the t 2g states are 1/6th filled [1,2] and the TiOCl system can be described by a low-energy multiband, t 2g Hubbard Hamiltonian [12]. The [TiO 4 Cl 2 ] octahedra are quite distorted, the t 2g states therefore show further splittings into lower d x y and higher d − = 1 √ 2 (d x z − d yz ) and d + = 1 √ 2 (d x z + d yz ) orbitals, with theẑ-axis pointing along the crystallographic a-axis andx-andŷ-axes rotated 45 • with respect to the crystallographic b-and c-axes. With this choice of axes, one has the usual convention of d x 2 −y 2 , and d z 2 pointing towards O and Cl neighbors. This however leads to a non-diagonal form of the Hamiltonian. The LDA density matrix (M) for TiOCl has then matrix elements between the d x z and d yz orbitals which are nonzero d x z |M|d yz = d yz |M|d x z = 0. A representation of the density of states in this basis reflects only the information of the diagonal matrix elements d x z |M|d x z = d yz |M|d yz (see figure 2 of [2]). Diagonalization of this matrix provides the corresponding representation into d − and d + . If one chooses instead a coordinate system with z = a,x = b,ŷ = c then the Hamiltonian is diagonal with d x 2 −y 2 , d yz and d x z forming the t 2g block and d x y , d z 2 forming the e g block (for a detailed account see [18]).
In the following, we will work in the coordinate system with d x y , d − and d + as t 2g states. The t 2g crystal field splittings calculated within NMTO are 1 (d x y ↔ d − ) = 0.29 eV and 2 (d x y ↔ d + ) = 0.56 eV which are in reasonable agreement with infrared absorption spectroscopy measurements [19] where 2 = 0.65 eV, and also comparable to recent cluster calculations for the TiO 4 Cl 2 octahedron [19] 1 = 0.25 eV and 2 = 0.69 eV.

LDA + DMFT calculations
Following our previous work in [12], we consider as the starting point for the description of TiOCl, the LDA-NMTO Hamiltonian, H LDA mm , to which direct and exchange terms of the screened onsite Coulomb interaction U mm and J mm of Hubbard-Hund type are added 6 , Here,ĉ σ † Rm (ĉ σ Rm ) denotes the creation (annihilation) operator for an electron at site R in orbital m with spin σ, andn σ Rm =ĉ σ † Rmĉ σ Rm is the particle number operator. U mm and J mm are parametrized as follows: U mm = U is the Coulomb repulsion between electrons in the same orbital, U m =m = U − 2J is the average repulsion, and J m =m = J is the Hund's rule coupling. We assume the double counting correction to be orbital independent within the t 2g states, thus resulting in a simple shift of the chemical potential. The choice of the local orbitals (Ti-t 2g Wannier functions) is done via the NMTO-downfolding technique [17], where all the partial waves other than Ti − d x y , Ti − d x z and Ti − d yz are downfolded (these notations refer to the choice of coordinate system as explained previously). As shown in [12] and earlier communications [20], the recent implementation of LDA + DMFT [21] allows the many-body Hamiltonian to be solved, including all off-diagonal elements in the orbital space of the local self-energy, mm . This has been crucial for TiOCl since the choice of the coordinate system as discussed above does not result in a diagonal form of the onsite matrix.
Our starting point is the high-temperature crystal structure. Since the low-temperature crystal structure of TiOCl shows a doubling of the cell along the crystallographic b-axis, a natural choice of Ti pairs are those along the b-axis. For our cluster calculation, we have therefore carried out a supercell calculation with the unit cell doubled along the b-axis. This results in four Ti atoms in the unit cell, marked as 1, 2, 3 and 4 in figure 1, with two Ti-Ti pairs located in the upper and lower layer of a given bilayer. Our self-energy will have accordingly off-diagonal elements in the orbital space as well as in the site space. The subscript M in the self-energy M M is defined as a composite site-orbital index. Such formulation has been already successfully applied in the case of VO 2 and Ti 2 O 3 [22,23].
For the present problem, a 12 × 12 block-diagonal self-energy matrix, M M is, therefore constructed where M is the composite index (m, i c ) with m denoting the orbital index, 1, 2, 3 5 for a t 2g basis and i c = (1, 2) and (3,4) denote the intradimer site indices for two Ti-Ti pairs in the unit cell. We neglect the interdimer correlation connecting the two pairs (1,2) and (3,4), and set the interdimer components of the self-energy in the present calculation to zero resulting in a block diagonal form of the self-energy. This may be a reasonable approximation considering the fact that the effective intradimer interaction between 1 and 2, or 3 and 4 is about an order of magnitude larger compared to interdimer interactions connecting two pairs belonging to different layers in this system [2]. With this approximation, the self-energy takes the form: whereˆ 11 andˆ 22 denote the on-site self-energy corresponding to sites 1 and 2 within a pair. Each of these matrices is therefore a 3 × 3 matrix.ˆ 12 (ˆ 21 ) gives the intersite, intradimer component of the self-energy. Note that the presence of the intersite component of the selfenergy,ˆ 12 (ˆ 21 ) gives rise to some 7 k-dependence of , as expected for a cluster-DMFT calculation.
We further assume that the dimers belonging to two different layers of the bilayer are similar, therefore the upper and lower sub-blocks of the matrix given by (1) are identical, i.e. 33 =ˆ 11 ,ˆ 44 =ˆ 22 and so on. It is therefore enough to work with a 6 × 6 block of this self-energy in the DMFT self-consistency condition and then construct the full 12 × 12 block from the 6 × 6 block to which H LDA mm is added to generate the Green's function. One may note that, with such a choice, the bond ordered ground state is allowed in the calculations, although one needs to include the electron-phonon interaction in the Hamiltonian in order to drive that state. For a single chain of Ti atoms running along the b-axis, the dimers may be formed in even or odd bonds. For a bilayer, as is the present case, there are, however, four possible arrangements [8], with even/even, odd/odd, odd/even and even/odd bonds in upper/lower layers. With our present choice of supercell, it is not possible to distinguish between these cases. Consideration of such scenarios would require to have a supercell which is four times enlarged along the b-axis. A second possible dimer pattern scenario, which is also not included in the present calculation, is the in-phase/out-of-phase dimer arrangements of neighboring chains within a layer. Such arrangements can be distinguished in a supercell calculation, which in addition to doubling along the b-axis is twice enlarged along the a-axis. Both calculations are presently computationally too expensive.
The 6 × 6 impurity problem is solved by a numerically exact quantum Monte Carlo (QMC) scheme. Within the given computational resources, we could reach a temperature of 1400 K with 10 6 QMC sweeps and 68 time slices. U and J values were chosen to be 4 and 0.7 eV respectively. Whereas J = 0.7 eV is a generally accepted value for an early transition-metal oxide (hardly changed from its atomic value), the precise determination of U is a delicate issue. Ab initio techniques such as constrained LDA give only rough estimates. For an early transition metal like Ti, U is expected to be between 3 and 5 eV [2]. We experimented in the past calculations [12] with different choices of U in terms of a single-site DMFT calculation and fixed U to 4.0 eV, since it provided the best possible spectra in terms of bandgap. We consider here the same U and J values as above so as to compare the single-site and cluster-DMFT results. The maximum entropy method [24] has been employed for the analytic continuation of the Green's function from the imaginary to the real axis for the calculation of the spectral functions.
Experimentally, the electron removal and addition spectra were probed by angle-integrated PES and x-ray absorption spectroscopy (XAS), respectively. Details of the PES experiment are described in Hoinkis et al [14]. XAS was measured at the PM-3 beamline of BESSY (Berlin, Germany) using its MUSTANG endstation. Due to finite Ti 3d-O 2p hybridization the oxygen K-edge spectrum can be taken as an approximate measure of the Ti 3d electron addition spectrum. The energy resolution for PES and XAS amounts to 100 and 50 meV, respectively. Figure 2 shows the Ti t 2g -dominated total spectral function computed with the NMTO + cluster-DMFT method described above (full line) in comparison with the angle integrated photoemission data (blue solid dots) and oxygen K-edge XAS measurements (red open dots) in the range of energies between −5 and 7 eV. Note that the experimental XAS spectrum shows both the t 2g and e g manifolds. The insulating behavior of TiOCl is correctly described by this calculation with a charge gap of about 1.1 eV, in rough agreement with the observed 2 eV optical gap [8,25]. The important point to note is the enhancement of the gap value compared to the single-site DMFT result which was about 0.3 eV [12]. This inevitably points towards the importance of the inter-site fluctuations for a good description of a low-dimensional system like TiOCl.  [12] with photoemission. To facilitate better comparison between experimental and theoretical lineshapes an inelastic background has been subtracted from the photoemission spectrum. All spectra are aligned to the same first order moment and scaled to the same integrated weight. The absolute energy scale refers to that of experiment.

Results
The calculated spectral weight distribution below and above the chemical potential µ, i.e. the lower and upper Hubbard bands compare reasonably well with the PES and oxygen K-edge absorption data, respectively, though one should keep in mind that the temperature used in the theoretical calculation was rather high compared to that of the experiments. We observe that the spectral function below µ is dominated totally by the d x y -like contribution with a 99% occupancy, consistent with polarization-dependent optical spectroscopy [8] and ARPES [14] experiments. Note that our calculations do not show any shape for the observed empty e g states since these bands were not explicitly considered in the cluster-DMFT calculations. They were downfolded and included only as tails of the t 2g NMTO Wannier orbitals and may be responsible for the slightly larger width of the calculated upper Hubbard band compared to the t 2g peak of the XAS spectrum.
In figure 3, we present a comparison of the cluster-DMFT with the single-site DMFT results [12]. Both calculations were performed within the same NMTO basis and considering the same QMC impurity solver. We notice improvement of the present results with respect to the single-site DMFT results. This may emphasize once more the importance of including dynamical Ti-Ti intersite correlations for the description of the spectral properties of TiOCl. For comparison, we also show the spectra obtained from LDA + U calculations which has a much narrower width as has been already stressed in [14].
While the accordance between theory and experiment turn out to be reasonably good as we have shown above, there are a few sources of improvement, which are all related to the computer feasibility of the calculations and are out of the scope of the present work. (i) The presented results were performed at quite high temperatures and with moderate numbers of QMC sweeps and time-slice steps regarding the QMC impurity solver. Calculations at lower temperatures are expected to provide a better resolution of the results. In particular, it will be necessary to check the temperature dependence of the spectra which is at the moment beyond our capability 8 . (ii) Consideration of larger supercells is expected to distinguish between the odd/even and even/even arrangements as well as in-phase/out-of-phase patterns and may provide a clue of the possible bilayer frustration present in this system. (iii) In order to be both computer realistic and manageable, we have used the NMTO t 2g Wannier functions for a low-energy subset of LDA bands, while in the actual case there are also O and Cl-p dominated bands and also higher lying Ti-e g states. The influence of the Cl and O-p-like bands (and also Ti-e g bands) has been considered implicitly in the construction of the t 2g Wannier functions which take care of the hybridization effects coming from Cl and O-p-like bands (see figure 3 in [2]), but they have not been considered explicitly. In principle, one should therefore carry out calculations involving O and Cl-p bands too, which is an extremely computer intensive job and not possible to carry out in the present framework.

Conclusions
By means of cluster-DMFT calculations implemented in the NMTO Wannier function basis, we have shown that the Ti-Ti intersite correlations in TiOCl play an important role for the proper description of photoemission and O K-edge XAS spectra. Although our calculations were limited by the computing feasibility, the improvement of the value of the optical gap and the overall comparison of the theoretical spectra with the experimental one indicates the significant role of the off-site, intra-dimer correlations. The present results suggest that the nature of fluctuations observed in a large variety of experiments [5,10] may be governed by correlated Ti-Ti dimers, existing up to rather high temperatures.