Formation of CuO 2 sublattices by suppression of interlattice correlations in tetragonal CuO

We investigate the tetragonal phase of the binary transition metal oxide CuO (t-CuO) within the context of cellular dynamical mean-field theory. Due to its strong antiferromagnetic correlations and simple structure, analysing the physics of t-CuO is of high interest as it may pave the way towards a more complete understanding of high-temperature superconductivity in hole-doped antiferromagnets. In this work we give a formal justification for the weak-coupling assumption that has previously been made for the interconnected sublattices within a single layer of t-CuO by studying the non-local self-energies of the system. We compute momentum-resolved spectral functions using a Matrix Product State (MPS)-based impurity solver directly on the real axis, which does not require any numerically ill-conditioned analytic continuation. The agreement with photoemission spectroscopy indicates that a single-band Hubbard model is sufficient to capture the material’s low energy physics. We perform calculations on a range of different temperatures, finding two magnetic regimes, for which we identify the driving mechanism behind their respective insulating state. Finally, we show that in the hole-doped regime the sublattice structure of t-CuO has interesting consequences on the symmetry of the superconducting state.

In this paper, we investigate the dynamical influence of the inter-sublattice hopping t d by the means of cellular dynamical mean field theory (CDMFT) [31] :::::::: [31-34] and motivate an efficient block-construction scheme for our cluster calculations. Our key finding is that the inter-sublattice correlations are heavily suppressed as compared to local and short-range intrasublattice correlations, which formally justifies to regard t-CuO as weakly-coupled interlaced CuO 2 lattices. Using a matrix product state [35,36] (MPS)-based impurity solver working directly on the real axis [37][38][39][40][41] and at effectively zero temperature we can reproduce equal energy maps and momentum resolved spectral functions in remarkable agreement with ARPES measurements without the need for analytic continuation. Further we analyse the magnetic ordering in t-CuO as a function of temperature and identify two driving mechanisms for the insulating phase. Finally, we predict the presence of superconductivity (SC) upon hole-doping by applying a complementary cluster technique, the variational cluster approximation (VCA) [42]. As a direct consequence of the sublattice decoupling, we find coexistence of magnetic stripe order and superconductivity of d x y -symmetry, whereas the usual cuprate d x 2 − y 2 order is strongly suppressed.

Model Hamiltonian
Each CuO plane of t-CuO is made up of edge-sharing CuO 4 plaquettes, which can be viewed as consisting of two interpenetrating CuO 2 square lattices. Following this logic, we consider one slab within the a − b plane as shown in Fig. 1(a). We consider a single band Hubbard model [43]: with i, j being site indices and σ ∈ {↑, ↓}. The single particle terms (t d = −0.1 eV, t = 0.44 eV, t = −0.2 eV, t = 0.075 eV) were obtained as a result of fitting the magnon dispersion, measured by RIXS, with a t-J model in Ref.
[24]. Contrary to the usual CuO 2 planes found in cuprate superconductors, the interstitial O atoms within one slab favor next-nearest neighbour (NNN) hopping t between Cu sites rather than nearest-neighbour (NN) hopping t d .
We use an Hubbard interaction strength of U = 7 eV, a significantly higher value than the one from Ref.
[24] but necessary for obtaining a gap that is larger than the experimental lower bound 2.35 eV [17]. Similar values of U have been used in LDA+U calculations [22,23].

Sublattice decoupling
At the single particle level it is hard to argue for the decoupling of the two sublattices since the nearest-neighbour hopping t d is of the same order of magnitude as the next-nearest neighbour hopping (t d ∼ − t 4 ). Therefore it is important to also take into account the self-energy which captures the modification of the non-interacting Hamiltonian due to the presence of electronic interactions in the correlated material.
To this aim, we choose the dimer cluster including the next-nearest neighbour hopping t ( Fig. 1 (d)) and the plaquette cluster containing two such dimers connected by the nextneighbour hopping t d (Fig. 1 (e)). The following results have been obtained by a MPS-based impurity solver [57-59] working on the imaginary axis and were computed using CDMFT at effectively zero temperature (T = 0 K). More details on the solver can be found in Sec. 4 and App. A.
In Fig. 3(a) we show an equal energy cut on the top of the valence band, which agrees well with the energy map measured in experiment (Ref. [17], Fig. 1(a)): We recover the strong maxima in the middle of the BZ, which are offset by 90 • . We also reproduce the replica features outside the single-sublattice BZ (dashed black line) that experimentally justified the assumption of only weakly coupled sublattices. To elaborate on this point in more detail, we note that on the one hand the two sublattices would be entirely decoupled only for t d = 0, yielding the spectral function of a single sublattice. In such a case, the features inside the first BZ of a single sublattice would be periodically replicated outside the BZ. On the other hand, the vanishing inter-sublattice self-energy (see Fig. 2(c)) keeps the hopping t d bare, whereas the intra-sublattice self-energy enhances the hopping t (see Fig. 2(b)) by a factor of ∼ 2. This effectively renders the t hopping ∼ 10 times stronger than t d , explaining the close resemblance of the replica features with respect to the ones in the original BZ. Note that unlike in ARPES [69] there are no matrix-element effects present in our calculation, which is why our replicas do not undergo any additional intensity modulations. In order not to favor any direc- tion by using an asymmetric super-cluster we average over possible cluster orientations. This procedure is described in more detail in App. D. The remaining difference between the x and y direction in Fig. 3 is entirely due to the magnetic stripe order. In panel (b), we show the momentum resolved spectral function of the valence band using the block-construction and compare it to the ARPES spectrum along the high symmetry path through the BZ depicted in Fig. 3(a). Comparing our results to ARPES (cf. Fig. 2(a) in Ref.
[17]) we find overall good agreement. In particular the low-energetic Zhang-Rice-like bands, which are separated from the lower Hubbard band at higher binding energy coincide (see inset of Fig. 3(b)). We identify this band to stem from a spin-polaron, i.e. a hole propagating in an antiferromagnetic background. , ::::::::: similarly :: to ::: the :::::::::::::: interpretation :: of ::::: Refs. :::::::::::: [68,[70][71][72] ::: for ::::::::: standard ::::::::: cuprates. : The incoherent and very dispersive spectrum without well-defined structures around the M and Γ points are also consistent with the measured spectrum. Moreover we reproduce the experimentally observed missing spectral weight at the X point, a feature which was not obtained within a self-consistent Born approximation calculation based on a Zhang-Rice singlet (ZRS) [10] spin-model [28]. An obvious feature that the calculations presented in this work can not reproduce are the contributions from a lower lying band marked with β in the experimental data [17], which is not included in our low-energy model. However, apart from these features the agreement between our model and the experiment is striking.

Finite temperature analysis
All results so far presented were computed at T = 0 K, however, there have been multiple predictions about the Néel temperature T N for the antiferromagnetic ordering of t-CuO in the literature [15,18,19,25], which underlines the necessity to better understand the finitetemperature behavior of the system. To this end, we employ CDMFT with a continuous-time Quantum Monte Carlo solver using the dimer and the block-construction clusters (Figs. 1(d),(g)).
In Fig. 4 we show the staggered magnetization as a function of temperature as well as the spinand momentum-resolved self-energy for three characteristic temperatures. On the right, we show the self-energy at the two cluster momenta K 1 = (0, 0) (dashed) and K 2 = 0, π a (solid) respectively.

Superconductivity away from half-filling
In order to address the question of superconductivity upon doping the system, we employed the variational cluster approximation (VCA) method [42,78,79]. This technique is particularly well suited to study the energetics of different symmetry breaking solutions of the model and their competition. It is based on finding stationary points of the self-energy functional Ω(Σ), which approximates the grand potential of the (lattice) system in the space of cluster self-energies [80]. These self-energies are parametrized via suitably chosen one-body parameters of the quantum cluster, potentially augmented by Weiss fields to allow for solutions with long-range order. We checked for different singlet-pairing channels and found stable solutions for superconducting Weiss fields of d x y and d x 2 − y 2 symmetry. These two pairing channels have been discussed already in the context of the t − t − U Hubbard model [81], which would correspond to only take into account t d and t. Whereas the d x 2 − y 2 channel and its competition with Néel antiferromagnetism are important close to half-filling for t/t d < 1, the (π, 0) collinear antiferromagnet and SC of d x y symmetry were identified to be key for t > t d [81,82]. Here, however, we focus on the superconducting channels of d x y and d x 2 − y 2 symmetry away from half-filling. For both d x y and d x 2 − y 2 symmetry, the superconducting solutions are energetically favoured over the normal state (PM) solution for fillings n < 1. The same is true for antiferromagnetic stripe order (AFS) , see Fig. 5, which we find even lower in energy down to n ≈ 0.87. However, when allowing for competition between these symmetry-breaking orders, a coexistence of superconductivity and AFS (AFS+SC) leads to the overall lowest energy solution at zero temperature, see red curve in Fig. 5. Comparing the corresponding antiferromagnetic and superconducting order parameters of these solutions shows that they are reduced in the coexistence solution as compared to the pure AFS and SC solutions. This indicates a competition between magnetic and superconducting orders upon doping. Most interestingly, the order parameter of the d x 2 − y 2 SC is strongly suppressed by the presence of antiferromagnetic stripe order such that the Cooper pairing is mainly of d x y -symmetry. Finally, we note that SC of d x y symmetry actually corresponds to d x 2 − y 2 symmetry within each of the two sublattices. Therefore, in the context of sublattice decoupling, our energetically most favourable solution could be interpreted as the emergence of a d x 2 − y 2 superconducting state with coexisting (Néel) antiferromagnetic order on each sublattice.  Fig. 1(g), i.e. without block-construction and optimized the functional in addition with respect to the (cluster) chemical potential µ (µ ).

Methods
Our method of choice to treat the interacting many-electron problem is a cluster extension of dynamical mean-field theory (CDMFT) [?, 31, 83] ::::::::::::::: [31-33, 44, 83]. In CDMFT, the full lattice problem is mapped to an effective cluster of several sites which is dynamically coupled to an electronic reservoir that represents the rest of the solid. This cluster-bath system is solved numerically for its Green's function and linked to the full lattice problem via a self-consistency loop. Due to the long-range magnetic stripe order, see Fig. 1(c), we perform magnetic calculations choosing cluster tilings, that are inline with the order. We propose several cluster geometries as shown in Fig. 1(d)-(g) for our CDMFT calculations, both to investigate the question of weak coupling ( Fig. 1(d,e)) as well as to obtain further observables like the spectral function ( Fig. 1(d,f,g)). In order to allow for a polarized solution, the CDMFT loop was initialized with a strongly polarized (constant) self-energy. To investigate possible superconducting solutions upon doping, we use the variational cluster approximation (VCA) [42,78,79], which is an established variational quantum cluster technique well suited to check for different symmetry-breaking orders of the lattice system [84,85]. Both techniques can be explained in terms of self energy functional theory [80,86] and are in this sense complementary [87]. They rely on the solution of the embedded cluster problem, for which we used different solvers as detailed below.

Imaginary axis MPS impurity solver
We use the MPS-based impurity solver introduced in Ref.
[57] and successfully applied in the context of DFT+DMFT in Refs. [58,59,88]. Using Matrix Product States (MPS) [35,36] we are able to access effectively zero temperature. MPS need a Hamiltonian formulation of the cluster impurity problem, which we obtain by following the fitting procedure introduced in Ref. [89] in the context of exact diagonalisation (ED). However with MPS we are able to treat a larger if the norm becomes smaller than 10 −8 . Over the course of the last two iterations the self-energies presented in Fig. 2 did not change by more than 7 · 10 −4 eV on the diagonals (Fig. 2 (a)), by not more than 5 · 10 −4 eV on the offdiagonal corresponding to the intra-sublattice hopping t (Fig. 2 (b)), and by not more than 6 · 10 −6 eV on the inter-sublattice component (Fig. 2 (c)). To improve the numerical accuracy for the computation of the high-Matsubara frequency part (tail) of self-energies we use the additional correlator introduced by Bulla et al. [111].

A.1.2 Real axis
For real axis calculations we use a broadening of η = 0.05 eV. In ground state searches we allow for a maximal bond dimension of 1536. We time-evolve using TDVP [109, 110] up until T max = 40 eV −1 (T max = 60 eV −1 ) for the block construction (dimer) calculations. Even though more costly in real time calculations we also use the additional correlator introduced by Bulla et al. [111] to improve the quality of the self-energy as compared to Dyson's equation. For calculations with the dimer cluster we use a mixing factor of 0.7 in the last few iterations.
to t d to be zero, which is a fair approximation as we confirmed when inspecting Fig. 2. Then, the two interlaced sublattices are no longer interconnected and the cluster self-energy is of block structure when regrouping cluster site indices that belong to the same sublattice. This is for instance the case for the differently coloured sites in Fig. 1(g). Therefore, it is sufficient to solve the impurity problem on one of the two sublattices. After closing the self-consistency loop we project down onto one of those blocks and obtain an impurity problem on the unit cell of a single sublattice ( Fig. 1(f) in the main text). Solving this problem yields one of the two blocks of the aforementioned self-energy. This construction amounts to a momentum resolution as obtained by considering the cluster depicted in Fig. 1(g), but with the computational effort of considering the unit cell shown in Fig. 1(f). The approximations described in this paragraph in essence correspond to treating the inter-sublattice hopping included in the diamond on the single-particle level via a feedback from the self-consistency loop [56]. A similar block construction of a supercluster was already used successfully within CDMFT [112,113].

D Averaging dimer/diamond orientations
Most of the results in this work were computed using the t-dimer or the block construction as unit cells. Fig. A2(a,b) shows two different orientations for these unit cells, that are in line with the stripe order proposed in Ref. [24]. Choosing one of them would artificially introduce an asymmetry, which is why the results presented in the paper were averaged over the two possible orientations. This approach goes by the name oriented cluster DMFT and was already introduced in Ref. [114,115] and applied to Sr 2 IrO 4 [114][115][116].
For completeness we show the equal energy maps obtained for every orientation compared to their respective mean in Fig. A2(c,d). We observe that the dimer results are far more sensitive to the orientation, however apart from the minimum in the middle of the BZ their average (d) Figure A3: Comparison of equal energy maps (a,c) and spectral functions on a path through the BZ (b,d) for the two different possible direction of the stripe order. The results shown were obtained with the block construction discussed in the main text. The heat maps were normalised to the maximal value shown. The path through the BZ is the one depicted in Fig. 3(a) of the main paper.
is already very similar to the energy maps computed with the block construction. This has two promising implications. First, it implies that the dimer results already capture very well the physics in t-CuO, indicating that the most important physical content of the extended unit cells is actually the delocalisation along the dominating bonds, as was already argued in the main text. On a second note we interpret the fact that the block construction result is almost independent of the orientation as an indication for convergence in cluster size.

E Stripe orientation
As mentioned in the main text we initialized our solvers such that the stripe order depicted in Fig. 1 is favoured. However there is no reason why the stripes should specifically be oriented in x-direction as opposed to y-direction. Thus in order to deliver a more complete picture we compare in Fig. A3 the results obtained when favouring the order in y-direction to those presented in the main text. In Ref.
[27] the possibility of having multiple domains in the ARPES sample was mentioned. This would yield a spectral function that amounts to some weighted superposition of the two spectra and equal energy maps shown in Fig. A3, where the weight would depend on the portion of the sample that has a certain :::::::: magnetic : stripe orientation. However as both :::::::: magnetic : stripe orders yield qualitatively very similar results we only discuss one of them in the main text.

F Estimation of the critical temperature
In the main text we mention an estimate for the critical temperature and also depict it with error bars in Fig. 4(a). In order to extract this estimate we fitted a function of the form to the staggered magnetization. Here γ, T c and β are fit parameters and θ is the Heaviside step function, that was added in order to make the fits more stable. Note that unlike the rest of the manuscript here β is the critical exponent of the transition, while β c in the following denotes the inverse critical temperature. By inspection of the self-energies in the transition region we find upper and lower boundaries for β c . We set the lower boundary such that the spin splitting vanishes and the upper one such that the imaginary part of the diagonal components of the self-energy tends to 0 as ω n → 0. By this criterion we identify β c = 16 eV −1 (β c = 20 eV −1 ) and β c = 13 eV −1 (β c = 17 eV −1 ) as upper and lower boundary for the dimer and block-construction clusters respectively.
Varying the upper and lower boundaries of the fit interval we obtain a collection of fits, of which we discard those, which either display a deviation bigger than 0.05 from any data point or which do not give β c in the region that was determined by inspection of the self-energies. Thus we end up with a collection of valid fits over which we average the resulting β c . The error bars in Fig. 4(a) correspond to the standard deviation in the set of valid fits.
The average values we obtain for the critical exponent are β = 0.44±0.15 (β = 0.66±0.34) for the block-construction (dimer) respectively. The errors are again determined as the standard deviation in the set of valid fits. Finally, we note that the exponents are in good agreement with the expected mean-field critical exponent of β = 0.5 [117].

G Comparing CTQMC and MPS data at low temperature
At lowest temperatures (e.g. for β = 50eV −1 ), the Green's functions obtained from applying the CTQMC solver and the MPS-based solver at T = 0 K coincide, see Fig. A4. It is important to point out that the only notable deviations occur in the off-diagonal part of the cluster Green's function, where the absolute value of G i, j (iω n ) is small (∼ 10 −3 ). In order to achieve good agreement between the two methods we ensured a sufficiently large sampling in the CTQMC solver: maximal sampling count was larger than 7 · 10 6 for the dimer clusters, and 2 · 10 6 for the block construction.

H Self-energy temperature dependence
While already being in the ordered phase, at finite temperatures (e.g. β = 20 eV −1 ) the real part of the self-energy in Fig. 4(b) in the main text shows an additional dynamic splitting at low Matsubara frequencies. This dynamical effect decreases when temperature is lowered. The static part given by the high-frequency tail however shows a constant increase. To identify the leading mechanism we derive in the following a simple toy-model able to capture this behaviour by including thermal effects only at the single-site level. Following the idea of Stepanov et al.
[74], we consider a single Hubbard site subject to a small external magnetic field representing the spin-exchange coupling between neighbouring spins