Model and Method

1Models and method¶

In this chapter we discuss the tools we have used for solving the two models that are relevant for us. These are, broadly, Kondo lattice type models relevant for the double perovskites, and the Hubbard model for the Mott transition. The first section introduces the models. This is followed by a discussion of the approximations we use to set up a tractable numerical problem. The section after discusses the numerical methods used. In the last section we describe the benchmarks against which we have tested our results.

1.1Models¶

We consider two kinds of models:

In both cases the first term (say $H_0$) defines the non-interacting (band) problem, ${\vec \sigma}$ is the electron spin operator. $H_{KLM}$ refers to a Kondo lattice model where the local moments ${\bf S}_i$ live on a frustrated lattice, while $H_{Hubb}$ refers to the Hubbard model with electron repulsion $U$. The KLM is a ‘two species’ model, involving spins and fermions, while the Hubbard model involves interacting electrons (at half filling in our case). For the KLM, the usual approach is to use a variant of classical Monte Carlo (assuming that the spins can be treated as classical) while the Hubbard model is solved via quantum Monte Carlo, cluster DMFT, or some form of cluster perturbation theory. We discuss the KLM problem first in the next section, and then describe our approach to the Hubbard problem.

1.2Kondo lattice systems¶

The Kondo lattice is a ‘two-species’ model, in which electrons delocalise on an underlying lattice, via hoppings $t_{ij}$, and couple locally to core spins ${\bf S}_i$. The core spins are supposed to arise from high spin states of electrons well below the Fermi level. The coupling between the electrons and spins, $J$, would be large in our system.

The electronic properties emerge as a function of spatial arrangements of spins. The probability $P\{{\bf S}_i\}$ of a spin configuration is decided by the electron free energy in that background.

To obtain $P\{{\bf S}_i\}$, one has to evaluate the trace over the fermions. This is analytically computable only in limiting cases, and in general one has to resort to some numerical scheme. We use an exact diagonalization based Monte Carlo (ED-MC) to sample the distribution $P\{{\bf S}_i \}$. The electronic properties are computed by diagonalizing $H_{KLM}$ in the equilibrium configurations of the sampled $P\{{\bf S}_i \}$. The detailed procedure is discussed in Section Real space Monte Carlo .

The Hubbard model¶

The Hubbard model describes interacting electrons on a lattice, and is defined by

$t_{ij}$ is the hopping amplitude, usually taken to be non-zero only between nearest neighbours $\langle ij\rangle$ on a given lattice. It can be non-zero for further neighbours as well and for a specific material should be determined via ab initio calculation.

The Hubbard problem looks very different from the classical spin Kondo lattice. Beyond weak coupling it has been traditionally handled via quantum Monte Carlo and exact diagonalization tools. Approximate schemes like DMFT also ultimately resort to these tools. These methods have a serious size limitation, despite the enormous increase in computing power over the last two decades, and are still not able to access complex magnetic states.

We use a new approach, discussed in the next section, which focuses on spatial order and its fluctuations and allows us to access large lattice sizes. It is approximate in nature, and maps the Hubbard problem at half filling to an effective Kondo lattice like problem.

1.3Approximations¶

‘Classical’ approximation in the Kondo lattice¶

The model for that we will use for the double perovskites involves two sets of spatially alternating electronic levels connected by a hopping amplitude, with the electrons on one of these sites coupling to a core spin. This model involves two simplifications with respect to the real material.

• a The conduction states in reality are orbitally degenerate.

• b The core spin is a quantum degree of freedom.

In a material like SFMO the bands near the Fermi level are formed by $3~t_{2g}$ ($xy$, $yz$, $zx$) orbitals which delocalize in $x-y$, $y-z$, $z-x$ planes respectively. Multiple orbitals could, in principle, have a qualitative effect on the physics, as it does in the manganites via the Jahn-Teller distortion. The double perovskites do not seem to involve strong lattice effects and the orbital degeneracy survives. In this situation a single band model can describe the qualitative physics. In addition, a comparison to ab initio results, which include the full electronic structure, confirms the usefulness of the single band approach.

The local moment in SFMO has magnitude $S=5/2$. So, while the original Kondo lattice model of $S=1/2$ moments coupled to electrons essentially involves quantum dynamics of the spins, the double perovskite case can be approximated by treating the core spins as classical. The quantum effects can be built in later via a $1/S$ expansion (Das et al. (2013)).

Static approximation in the Hubbard model¶

The presence of two body interaction requires the use of approximations, or numerical tools, when studying the large $U$ behaviour of the Hubbard model. We opt for an auxiliary field decomposition that captures the mean field ground state but retains all thermal fluctuations. We use a Hubbard-Stratonovich (HS) transformation, below, to decouple the interaction in terms of auxiliary fields $\phi_i$ and ${\bf m}_i$ which couple to charge and spin densities respectively. The specific form of the decomposition was suggested by Hubbard (Hubbard (1959), Schulz (1990)) and explored further by Schulz (Schulz (1990)). It is motivated by the need to capture Hartree-Fock theory at saddle point, and retain the correct spin rotation invariance (and Goldstone modes). In terms of the Grassmann field $\psi_{i \sigma}(\tau)$, we have:

where charge and spin densities are defined as,

The partition function is now written as

where, the Lagrangian ${\cal L}(\tau)$ is (all operators have $\tau$ dependence)

At half filling $(n=1)$, we make the following two approximations:

• ($a$) We freeze the $\phi_i$ to its saddle point value $(U/2) \langle n_i \rangle = U/2$. This is motivated by the large $U$ which in the half-filled situation would penalise charge fluctuations. It is worth noting that away from half filling there would be $\phi_i$ fluctuations even at large $U$.

• ($b$) We neglect the $\tau$ dependence of ${\bf m}_i(\tau)$ field but retain its full spatial dependence. This allows unbiased choice of the $T=0$ magnetic state, and retains the thermal spin fluctuations on that state.

These two approximations lead to a model of electrons coupled to a the spatially fluctuating classical field ${\bf m}_i$. After making $\bf m_{i}$ dimensionless, by simple scaling ${\bf m}_{i} \rightarrow \frac{U}{2}{\bf m}_{i}$, the half filled Hubbard model is mapped to the following effective Hamiltonian.

where $\vec{\sigma}_i=\sum_{a,b}c^{\dagger}_{ia}\vec{\sigma}_{ab}c_{ib}=2{\bf s}_i$. This model ((9)) has a classical part $H_{cl}$, and electronic part $H_{el}$ describing electrons coupled to ${\bf m}_i$ fields, where -

This resembles the KLM, with the crucial difference that the ${\bf m}_i$ are not of fixed magnitude. Nevertheless, the Kondo lattice like model for the double perovskites, and the approximate model derived from Hubbard, can be solved by the same tools.

Neglecting time dependence in the Hubbard case has several physical consequences that we will discuss later. In terms of advantages it allows:

• The access to non trivial magnetic states and captures their $T_c$ scales accurately.

• It can access spectral and transport properties without any need for analytic continuation.

• It is not tied to long range order, and can describe a spin disordered Mott state and correlated metal as well.

Consequences of the approximations¶

The results of the KLM with $S=5/2$ are well approximated by treating the spins as classical. Quantum effects, for example spin waves, can be extracted via a $1/S$ expansion (Das et al. (2013)).

For the Hubbard model, as we mentioned above, and would see in detail in Chapter Mott transition on the anisotropic triangular lattice, our approximate scheme allows access to the finite temperature physics rather well. However, the neglect of quantum dynamics of the auxiliary fields ${\bf m}_i$ results in: (a) missing the correlation effects in the ground state of the metal, overestimating its energy and underestimating the critical $U/t$ for metal-insulator transition, and (b) missing a possible quantum spin liquid insulator (Sahebsara & Senechal (2008), Yang et al. (2010)) at intermediate $U/t$. We will discuss these issues in the relevant chapter.

1.4Numerical Methods¶

Both the models, the KLM (equation (\ref{eq:klm})) and effective Hubbard model (equation (9)), have one feature in common. They both describe the motion of electrons in the presence of spatially fluctuating classical fields, which are localized spins for KLM, and the auxiliary fields for the effective Hubbard model. Both problems are to be solved at strong coupling ($J/t \gg 1$ for KLM, or $U/t \gg 1$ for Hubbard model), which is well beyond the perturbative regime. In the absence of analytical tools one resorts to numerical methods. The numerical methods we use are broadly the following:

• (a) Real space Monte Carlo, where one ‘anneals’ the classical fields from high to low temperature and samples the distribution of the classical field using Metropolis sampling.

• (b) Variational minimization, in which we use a family of spatially periodic configurations of the classical variable and minimize the electronic energy over them.

As we discussed, the KLM (equation (2)) and effective Hubbard model (equation (9)), both describe electron propagation in a classical background of spins $\{{\bf S}_i\}$ of fixed magnitude (KLM), or the auxiliary fields $\{{\bf m}_i\}$ (effective Hubbard model). For a given configuration $\{{\bf m}_i\}$ (or $\{{\bf S_i}\}$), the electronic Hamiltonian $H_{el}$ is quadratic in electronic operators, and one just needs to diagonalize it to get the single particle eigenvalues. The total energy of the system is sum of (i) the classical part and (ii) the electronic part which is sum of eigenvalues up to the chemical potential at $T=0$, or the free energy at finite $T$.

However, the configurations $\{{\bf m}_i\}$ (or $\{{\bf S_i}\}$), themselves have to be determined from a certain distribution. At finite temperature, the thermal distribution of the fields is given by $P\{{\bf m}_i\}$ (or $P\{{\bf S}_i\}$):

The fermion trace is not analytically calculable when the magnetic coupling term is comparable to the kinetic energy and we need to converge towards the $P\{{\bf m}_i\}$ via an iterative combination of the Metropolis algorithm and exact diagonalisation. This is discussed next.

Real space Monte Carlo¶

We set up the electronic Hamiltonian in real space basis for a given lattice size for an initial configuration of the classical field and attempt to ‘update’ this configuration. The energy cost of this update has to be computed through the free energy change of the entire system. The transition probability between configurations is $\propto \exp{\left(-\beta(F'-F)\right)}$, where $F'$ and $F$ are free energy in the attempted and initial configuration respectively. Computing $F'-F$ requires diagonalizing $H_{el}$. This involves an ${\cal O}(N^3)$ computational cost per update, i.e., the cost per MC sweep would be $N^4$. This limits the accessible lattice size to $N\sim 100$.

To access larger sizes within reasonable time we use a cluster algorithm (TCA) (Kumar & Majumdar (2006)) for estimating the update cost. Here, rather than diagonalize the full $H_{el}$ for every attempted update, we compute the update cost by diagonalizing a cluster (of size $N_c \sim 64$, say) around the reference site. This scheme has been extensively benchmarked (Kumar & Majumdar (2006)) and works reasonably unless we are attempting to recover states with a large spatial period. A cartoon of the cluster scheme is shown in Figure 1. The computational cost now scales as $N{\cal O}(N_c^3)$ which is linear in lattice size $N$. This lets us access large lattice sizes ($24\times 24$ in 2D (with $N_c=8^2$), and $12\times 12\times 12$ (with $N_c=4^3$) in 3D.

Once equilibrium is attained within the Monte Carlo scheme at a given temperature we use the MC configurations to calculate various properties of interest. Below we list some of the key indicators that we track in the KLM and Hubbard model.

Auxiliary field properties¶

• $P(m)$: This describes the magnitude distribution of the ${\bf m}_i$ fields:

  $$P(m) = \frac 1N \sum_{i}\langle\delta(m - |{\bf m}_i|)\rangle$$

  (13)

  The sum is over all the lattice sites, and $\langle .. \rangle$ denote the average over thermally sampled field configurations. This quantity is irrelevant for KLM model, as the spins are fixed magnitude fields. For the Hubbard model it provides crucial insight into Mott physics.

• Structure factor $S({\bf q})$: The angular correlation is accessed via the thermally averaged structure factor $-$

  $$S({\bf q}) = \frac{1}{N^2}\sum_{ij}\langle{\bf m_{i}} \cdot{\bf m_{j}}\rangle e^{i{\bf q}\cdot({\bf X}_i-{\bf X}_j)}$$

  (14)

  at each temperature. The onset of rapid growth in $S({\bf q})$ at some ${\bf q} = {\bf Q}$, say, indicates a magnetic transition. The electronic properties are calculated by diagonalizing $H_{el}$ on the full lattice for equilibrium $\{{\bf m}_i\}$ configurations. The nature of these correlations is characterized by the ${\bf q}$ at which the maximum.

Electronic properties¶

• Density of states (DOS) $N(\omega)$ at a given temperature is obtained by taking thermal average of the single configuration DOS over equilibrium ${\bf m}_{i}$ samples~$-$

  $$N(\omega) = (1/N) \sum_{\alpha}\langle\delta(\omega - \epsilon_\alpha)\rangle$$

  (15)

  The $\epsilon_{\alpha}$ above are the eigenvalues in a single equilibrium configuration.

• Optical Conductivity $\sigma(\omega)$: The conductivity is calculated using the Kubo formula (ref. Allen (2006)), which takes the form for a non-interacting disordered system

  $$\sigma^{xx}(\omega) = \frac{\sigma_{0}}{N}\sum_{\alpha,\beta} { {n_{\alpha}-n_{\beta}} \over {\epsilon_{\beta}-\epsilon_{\alpha}} } |\langle \alpha|J_x|\beta\rangle|^2 \delta(\omega-(\epsilon_{\beta}-\epsilon_{\alpha}))$$

  (16)

  where the current operator $\vec{\mathbf J}$ is

  $$\begin{align*} \vec{\mathbf J} = -i\sum_{i,\sigma,\vec{\delta}}\left[\vec{\mathbf \delta}t_{\vec{\mathbf\delta}} ~c^{\dagger}_{{\mathbf X}_i,\sigma}c_{{\mathbf X}_i+\vec{\mathbf\delta},\sigma}-\textrm{h.c.})\right] \end{align*}$$

  (17)

  The d.c conductivity is the $\omega \rightarrow 0$ limit of the result above. $\sigma_{0}$=$\frac{\pi e^2}{\hbar}$ in 2D and $\frac{\pi e^2}{\hbar a_0}$ in 3D ($a_0$ is the lattice spacing in the relevant direction). $n_{\alpha}=f(\epsilon_{\alpha})$ is the Fermi function, and $\epsilon_{\alpha}$ and $|\alpha\rangle$ are respectively the single particle eigenvalues and eigenstates of $H_{eff}$ in a given background {${\bf m}_i$}.

• Spectral function $A({\bf k},\omega)$: We extract the thermal and spin averaged spectral function $A({\bf k},\omega)$ as follows. The retarded Greens function

  $$G_{\sigma}({\bf k},t)=-i\theta(t)\langle \{c_{{\bf k}\sigma}(t),c^{\dagger}_{{\bf k}\sigma}(0)\}\rangle$$

  (18)

  can be simplified to the following

  $$G_{\sigma}({\bf k},t)=-i\theta(t)\sum_{\alpha} |\langle {\bf k}\sigma|\alpha \rangle|^2e^{-i\epsilon_{\alpha}t}$$

  (19)

  where $\{|\alpha\rangle\}$ are the single particle eigenstates and $\epsilon_{\alpha}$ are eigenvalues in a given $\{ {\bf m}_i \}$ background. Doing the Fourier transform to frequency, one gets

  $$G_{\sigma}({\bf k},\omega)= \sum_{\alpha}\frac{|\langle {\bf k}\sigma|\alpha \rangle|^2}{ \omega-\epsilon_{\alpha}+i0^+}$$

  (20)

  Then, we have the spectral functions as

  $$A_{\sigma}({\bf k},\omega)=-{1 \over \pi}\text{Im}G_{\sigma}({\bf k},\omega) = \sum_{\alpha}|\langle {\bf k}\sigma|\alpha \rangle|^2\delta(\omega-\epsilon_{\alpha})$$

  (21)

  We average this over the sampled thermal configurations. The low frequency part of the spectral function $A_{\sigma}({\bf k}, \omega=0)$ serves to identify the shape of the Fermi surface, and the hot and cold spot locations.

Variational minimization at $T=0$¶

When the Monte Carlo suggests that the system is heading towards a state with periodic spatial behaviour, a simple way to access the possible ground state is to try a family of variational configurations. For simple enough configurations (discussed further on) the electronic model may be analytically solvable. Even when that is not possible, the numerical cost of diagonalising a system of size $N$ for $N_{trial}$ configurations is $\sim N_{trial} N^3$, compared to the $N_{tempr} N_{sweep} N^4$ cost involved in the MC. $N_{tempr}$ is the number of temperature points and $N_{sweep}$ is the number of MC sweeps per temperature.

Taking our representation of the Hubbard model (equation (9)) and assuming a set of periodic configurations ${\bf m}_i$ defined by a magnitude $m$ and period ${\bf q}$:

we have

This can be simplified by doing Fourier transformation to

Where $\epsilon_{\bf k}$ is the tight binding dispersion for the lattice. For an anisotropic triangular lattice this is given by $\epsilon_{\bf k} = -2t(\cos k_x+\cos k_y) - 2t'\cos (k_x+k_y)$. The eigenvalues can be readily obtained, and one gets

Once we get the dispersion, the total energy is

To get the ground state, we minimize the total energy $E(m,{\bf q})$ with respect to the magnitude ‘$m$’ and vector ${\bf q}$. On a two dimensional $L\times L$ lattice, ${\bf q}$ can take $L^2$ values, as ${\bf q} = \frac{2\pi}{L} (q_x,q_y)$, where $q_x$,$q_y$ are integers ranging from 0 to $L-1$. One can discretise $m$ (bounded between $[0,1]$) into $\sim 100$ intervals. As we will see in the solution of the 2D Hubbard model this provides some intuition into the nature of possible magnetic order.

1.5Benchmarks¶

Here we discuss selected benchmarks for our methods. We first provide a quick comparison of our results with results from other methods (to be elaborated in individual chapters later), and then discuss issues of size dependence.

Comparison with other methods¶

For the double perovskite model in three dimensions there is no other theoretical effort probing the regime that we explore in this thesis. We will discuss the issue of size dependence of our results in the appropriate chapter.

For the Hubbard model in 2D, on the anisotropic triangular lattice, there is a significant body of work. We had summarised these in the first chapter. Here we provide a quick comparison of these results with the ground state obtained by our approach. This is to establish the basic usefulness of the auxiliary field approach. Detailed indicators at finite temperature will be taken up in Chapter 4.

The principal methods for the ground state include variational Monte Carlo (VMC) used by Luca et al (Tocchio et al. (2009), Tocchio et al. (2013)), exact diagonalization (ED) by Clay et al (Clay et al. (2008)), and variational cluster perturbation theory (VCPT) by Sahebsara et al (Sahebsara & Senechal (2008), Sahebsara & Senechal (2006)).

Figure 19 shows the comparison of results from these methods with that from the auxiliary field based minimisation (equivalent at $T=0$ to unrestricted Hartree-Fock (UHF) theory). The qualitative trends in the UHF result match with that in the more sophisticated approaches. At finite $T$, where the wavefunction based approaches have limited use, the auxiliary field scheme would have a distinct advantage.

Most of the finite temperature studies have employed dynamical mean field theory (DMFT) and its cluster variant (C-DMFT). DMFT recovers some of the generic features in transport (Limelette & others (2003)) and optics (Merino & others (2008)) but does not capture the impact of magnetic correlations. C-DMFT has not been employed to explore response functions yet. We will compare our results to these wherever possible.

On the FCC lattice, there is hardly any theoretical work, except for a very early attempt using conventional static mean field theory (Grensing et al. (1978)). A recent experimental paper (Phuoc & others (2013)) does present a DMFT based phase diagram but that again does not retain any of the magnetic correlations that would be important close to the Mott transition.

Size dependence within the MC¶

In this section we provide a discussion of the size dependence of our results in the 2D anisotropic triangular lattice Hubbard model, and establish that the $24 \times 24$ lattice that we use in Chapter. 4 is adequate for most physical properties. This will allow us to focus on the physics issues there without digressing into questions of finite size.

In what follows we show results on the anisotropic triangular lattice at $t'/t=0.8$ and focus on the magnetic indicators. Electronic properties are controlled by these.

Figure 3 shows the magnitude distribution $P(m)$ of the auxiliary fields ${\bf m}_i$, at two temperatures ($T/t=0.1,0.005$) for four representative values, $U/t = 4,5,6,8$ across the Mott transition, for lattice sizes $L \times L$ with $L$ ranging from 8 to 24. The distribution is broad at high temperature and sharply peaked at low temperature. The $P(m)$ is practically insensitive to lattice size $L$ at all $U$ at high temperature. At low temperature and large $U$ the results are again independent of $L$ but at the lowest $U$, close to the Mott transition value, the $L=8$ data differs from the larger $U$ results. Beyond $L=12$ however even here the results are essentially size independent.

We probed the underlying state more closely in terms of magnetic correlations, rather than the gross size distribution. The allowed values of momentum change with system size. On the $L\times L$ lattice the allowed ${\bf q}$ are given by ${\bf q} = 2\pi(\frac{q_x}{L},\frac{q_y}{L})$ where $q_x,q_y$ are integers in the range $0,1,2,...L-1$. In Figure 4 we show the $q_x$, say, that are available for varying lattice size. For each choice of $L$ the points on the vertical line indicate this.

The next figure shows the temperature dependence of the peak in the magnetic structure factor for two values of $U$ and a wide range of $L$. The left panel establishes that there is no magnetic order (and in fact no magnetic moment) at $T=0$ in the $U=4$ problem. The $L=8$ case shows a weak correlation, associated with the finite size effect highlighted for $P(m)$. The strength of correlation at finite temperature does depend on $L$ but is reasonably stable by the time $L=24$. For $U=8$, right panel, the lower $L$ values give significantly enhanced values, and strong size dependence, in the $S({\bf q})$ peak, but the results stabilise for $L \ge 20$.

The full ${\bf q}$ dependence of $S({\bf q})$ is quite instructive. Figure 6a shows the result at $U/t=4$ and Figure 6b at $U/t=8$, in each case for two temperatures, $T/t=0.02$ and $T/t=0.1$. The panels from left to right are for increasing $L$ as indicated.

For $U/t=4$, low $T$ suggests ordering at a wavevector $\sim (0.8 \pi, 0.8 \pi)$ but this is systematically diminished with growing size. There are no moments and no order in the large $L$ system. The higher $T$ result suggests short range correlation and that pattern is reasonably stable with $L$. At this $U/t$ we need $L \ge 20$ to get a reliable picture of the ground state.

At larger $U$, Figure 6b, the moments are stable but the ordering pattern depends on $L$. $L=8$ has a peak at $(\pi, \pi)$, and weak satellites, suggesting Neel order. With increasing $L$ however the peak shifts away from Neel order to ${\bf q} \sim (0.8 \pi, 0.8 \pi)$, a spiral state (weak a weaker peak in $S({\bf q})$). Again for $L \ge 20$ the pattern is stable. Accessing the correct ground state even deep in the Mott window requires fairly large $L$. The high temperature data has a ${\bf q}$ dependence that is similar at all $L$, and is quantitatively stable for $L \ge 20$.

Overall, conclusions about the magnetic ground state and fluctuations can depend strongly on system size, but for $L=24$ that we employ in 2D in most of the calculations the results are free of finite size artifacts.

Variational minimisation also involves a size dependence and we discover that for $L \ge 20$ the VC results are reasonably stable and the Monte Carlo and VC based estimates are mutually consistent. Figure 7 compares the $L$ dependence of the ordering wavevector at two $U$ values and suggests that for $L=24$ the conclusions are similar. Most of our MC results are at $L=24$ but we have checked some data points at $L=30$ as well.

The checks above will allow us to focus on the physics that emerges from our calculation without having to check for finite size artifacts frequently.