Metal insulator transition

The half-filled Hubbard model undergoes a transition from a metallic state to an (Mott) insulating state when increasing the interaction strength \(U\). One model where this can be solved exactly is the Hubbard model on the Bethe graph in the limit of infinite connectivity. In this limit Dynamical Mean Field Theory (DMFT) is exact and can be used to study the transition.

The DMFT lattice self-consistency relation on the Bethe graph has the form

\[\Delta_\sigma(\tau) = t^2 G_\sigma(\tau)\]

where \(\Delta_\sigma(\tau)\) is the hybridization function, \(G_\sigma(\tau)\) is the local single particle Green’s function with spin \(\sigma \in \{\uparrow, \downarrow\}\), and we set the nearest neighbour hoppint \(t\) to unity (\(t=1\)).

At low temperature the transition is first order, but at higher temperatures it turns into a cross over. Here we will study this crossover at inverse temperature \(\beta = 5\).

Solver setup

Since we will be performing self-consistent DMFT calculations we modularize the solver setup using a helper function setup_solver().

[2]:

from triqs_soehyb import Solver

from triqs.gf import make_gf_dlr_imtime, make_gf_dlr_imfreq, SemiCircular

def setup_solver():

    S = Solver(
        beta=5.0, # Inverse temperature
        gf_struct=[('up', 1), ('do', 1)], # Green's function structure 1st index: name, 2nd index: dimension of subspace
        eps=1e-12, # Accuracy of Discrete Lehmann Representation (DLR) used for imaginary time response functions
        w_max=10.0, # DLR freqiency cut-off (the spectrum of the model must be in the range [-w_max, +w_max]
        )

    for bidx, delta_tau in S.Delta_tau:
        delta_w = make_gf_dlr_imfreq(delta_tau)
        delta_w << SemiCircular(2.0)
        delta_tau[:] = make_gf_dlr_imtime(delta_w)

    return S

Warning: could not identify MPI environment!

Starting serial run at: 2025-08-19 17:41:44.159873

Local Hubbard interaction

The local interaction in the Hubbard model has the form

\[H_{int} = U n_\uparrow n_\downarrow - \mu ( n_\uparrow + n_\downarrow )\]

where \(U\) is the so called Hubbard-\(U\) and \(\mu\) is the chemical potential, which at half filling is given by \(\mu = U/2\). The helper function get_h_int(U) constructs the local Hamiltonian using Triqs second-quantization operators.

[3]:

from triqs.operators import n

def get_h_int(U, mu=None):
    if mu is None:
        mu = U / 2
    h_int = U * n('up',0) * n('do', 0) - mu * ( n('up', 0) + n('do', 0) )
    return h_int

h_int = get_h_int(U=5.0)

print(f'h_int = {h_int}')

h_int = -2.5*c_dag('do',0)*c('do',0) + -2.5*c_dag('up',0)*c('up',0) + 5*c_dag('do',0)*c_dag('up',0)*c('up',0)*c('do',0)

DMFT self-consistency

Finally we design a function solve_dmft_self_consistency(S, h_int, order, ...) that takes a solver S, a local Hamiltonian h_int, the wanted expansion order, and performs the DMFT self-consistency iterations.

[4]:

def solve_dmft_self_consistency(S, h_int, order,
                                dmft_maxiter=100, dmft_tol=1e-6):

    for dmft_iter in range(1, dmft_maxiter+1):

        S.solve(h_int=h_int, order=order, maxiter=100, verbose=False)

        diff = np.max(np.abs(S.Delta_tau['up'].data - S.G_tau['up'].data))

        print(f'DMFT: dmft_iter = {dmft_iter}, diff = {diff:2.2E}')

        if diff < dmft_tol: break

        S.Delta_tau[:] = S.G_tau # DMFT Bethe lattice self consistency

    S.d = S.S.get_expectation_value(n('up', 0) * n('do', 0)).real
    S.n_up = S.S.get_expectation_value(n('up', 0)).real
    S.n_do = S.S.get_expectation_value(n('do', 0)).real

    return S

Solve NCA and OCA for a sweep in \(U\)

To study the metal-insulator cross over we solve the Hubbard model on the Bethe graph for a range of \(U\) values at first and second expansion order (a.k.a. the non-crossing and one-crossing approximation NCA/OCA).

[5]:

orders = [1, 2]
Us = np.arange(3.0, 7.0, 1.0)

for order in orders:

    S = setup_solver()

    for U in Us:

        print(f'--> order = {order}, U = {U}')

        S.U = U
        h_int = get_h_int(U)

        S = solve_dmft_self_consistency(S, h_int, order)

        filename = f'data_order_{order}_U_{U:011.8f}.h5'
        print(f'--> Storing: {filename}')
        with HDFArchive(filename, 'w') as A:
            A['S'] = S

  ___      ___    _  ___   _____
 / __| ___| __|__| || \ \ / / _ )
 \__ \/ _ \ _|___| __ |\ V /| _ \
 |___/\___/___|  |_||_| |_| |___/  [github.com/TRIQS/soehyb]

beta = 5.0, lamb = 5.00E+01, eps = 1.00E-12, N_DLR = 27
fundamental_operators = [1*c('up',0), 1*c('do',0)]
H_mat.shape = (4, 4)
H_loc = 0

--> order = 1, U = 3.0
DMFT: dmft_iter = 1, diff = 1.27E-01
DMFT: dmft_iter = 2, diff = 1.24E-02
DMFT: dmft_iter = 3, diff = 1.58E-03
DMFT: dmft_iter = 4, diff = 3.43E-04
DMFT: dmft_iter = 5, diff = 6.34E-05
DMFT: dmft_iter = 6, diff = 1.24E-05
DMFT: dmft_iter = 7, diff = 2.39E-06
DMFT: dmft_iter = 8, diff = 4.62E-07
--> Storing: data_order_1_U_03.00000000.h5
--> order = 1, U = 4.0
DMFT: dmft_iter = 1, diff = 4.31E-02
DMFT: dmft_iter = 2, diff = 1.11E-02
DMFT: dmft_iter = 3, diff = 3.24E-03
DMFT: dmft_iter = 4, diff = 9.52E-04
DMFT: dmft_iter = 5, diff = 2.83E-04
DMFT: dmft_iter = 6, diff = 8.44E-05
DMFT: dmft_iter = 7, diff = 2.52E-05
DMFT: dmft_iter = 8, diff = 7.50E-06
DMFT: dmft_iter = 9, diff = 2.24E-06
DMFT: dmft_iter = 10, diff = 6.68E-07
--> Storing: data_order_1_U_04.00000000.h5
--> order = 1, U = 5.0
DMFT: dmft_iter = 1, diff = 3.63E-02
DMFT: dmft_iter = 2, diff = 5.84E-03
DMFT: dmft_iter = 3, diff = 1.63E-03
DMFT: dmft_iter = 4, diff = 4.16E-04
DMFT: dmft_iter = 5, diff = 1.06E-04
DMFT: dmft_iter = 6, diff = 2.69E-05
DMFT: dmft_iter = 7, diff = 6.85E-06
DMFT: dmft_iter = 8, diff = 1.74E-06
DMFT: dmft_iter = 9, diff = 4.42E-07
--> Storing: data_order_1_U_05.00000000.h5
--> order = 1, U = 6.0
DMFT: dmft_iter = 1, diff = 3.04E-02
DMFT: dmft_iter = 2, diff = 2.63E-03
DMFT: dmft_iter = 3, diff = 5.08E-04
DMFT: dmft_iter = 4, diff = 1.05E-04
DMFT: dmft_iter = 5, diff = 2.09E-05
DMFT: dmft_iter = 6, diff = 4.09E-06
DMFT: dmft_iter = 7, diff = 8.00E-07
--> Storing: data_order_1_U_06.00000000.h5
  ___      ___    _  ___   _____
 / __| ___| __|__| || \ \ / / _ )
 \__ \/ _ \ _|___| __ |\ V /| _ \
 |___/\___/___|  |_||_| |_| |___/  [github.com/TRIQS/soehyb]

beta = 5.0, lamb = 5.00E+01, eps = 1.00E-12, N_DLR = 27
fundamental_operators = [1*c('up',0), 1*c('do',0)]
H_mat.shape = (4, 4)
H_loc = 0

--> order = 2, U = 3.0
DMFT: dmft_iter = 1, diff = 8.02E-02
DMFT: dmft_iter = 2, diff = 3.40E-03
DMFT: dmft_iter = 3, diff = 3.65E-04
DMFT: dmft_iter = 4, diff = 5.09E-05
DMFT: dmft_iter = 5, diff = 7.03E-06
DMFT: dmft_iter = 6, diff = 9.69E-07
--> Storing: data_order_2_U_03.00000000.h5
--> order = 2, U = 4.0
DMFT: dmft_iter = 1, diff = 4.58E-02
DMFT: dmft_iter = 2, diff = 1.46E-02
DMFT: dmft_iter = 3, diff = 6.23E-03
DMFT: dmft_iter = 4, diff = 2.71E-03
DMFT: dmft_iter = 5, diff = 1.20E-03
DMFT: dmft_iter = 6, diff = 5.39E-04
DMFT: dmft_iter = 7, diff = 2.42E-04
DMFT: dmft_iter = 8, diff = 1.09E-04
DMFT: dmft_iter = 9, diff = 4.91E-05
DMFT: dmft_iter = 10, diff = 2.22E-05
DMFT: dmft_iter = 11, diff = 9.99E-06
DMFT: dmft_iter = 12, diff = 4.50E-06
DMFT: dmft_iter = 13, diff = 2.03E-06
DMFT: dmft_iter = 14, diff = 9.16E-07
--> Storing: data_order_2_U_04.00000000.h5
--> order = 2, U = 5.0
DMFT: dmft_iter = 1, diff = 4.13E-02
DMFT: dmft_iter = 2, diff = 1.34E-02
DMFT: dmft_iter = 3, diff = 6.27E-03
DMFT: dmft_iter = 4, diff = 2.80E-03
DMFT: dmft_iter = 5, diff = 1.23E-03
DMFT: dmft_iter = 6, diff = 5.41E-04
DMFT: dmft_iter = 7, diff = 2.37E-04
DMFT: dmft_iter = 8, diff = 1.04E-04
DMFT: dmft_iter = 9, diff = 4.56E-05
DMFT: dmft_iter = 10, diff = 2.00E-05
DMFT: dmft_iter = 11, diff = 8.75E-06
DMFT: dmft_iter = 12, diff = 3.84E-06
DMFT: dmft_iter = 13, diff = 1.68E-06
DMFT: dmft_iter = 14, diff = 7.37E-07
--> Storing: data_order_2_U_05.00000000.h5
--> order = 2, U = 6.0
DMFT: dmft_iter = 1, diff = 3.30E-02
DMFT: dmft_iter = 2, diff = 5.61E-03
DMFT: dmft_iter = 3, diff = 1.83E-03
DMFT: dmft_iter = 4, diff = 6.47E-04
DMFT: dmft_iter = 5, diff = 2.18E-04
DMFT: dmft_iter = 6, diff = 7.20E-05
DMFT: dmft_iter = 7, diff = 2.37E-05
DMFT: dmft_iter = 8, diff = 7.77E-06
DMFT: dmft_iter = 9, diff = 2.55E-06
DMFT: dmft_iter = 10, diff = 8.34E-07
--> Storing: data_order_2_U_06.00000000.h5

Load result

To visualize the crossover we load the results from file using a helper function and class.

[6]:

import glob

class AttrList:
    def __init__(self, l): self.l = l
    def __getattr__(self, attr): return [ getattr(x, attr) for x in self.l ]


def load_order(order):

    filenames = np.sort(glob.glob(f'data_order_{order}_U_*.h5'))

    results = []
    for filename in filenames:
        print(f'--> Loading: {filename}')
        with HDFArchive(filename, 'r') as A:
            results.append(A['S'])

    return AttrList(results)

nca = load_order(1)
oca = load_order(2)

--> Loading: data_order_1_U_03.00000000.h5
--> Loading: data_order_1_U_04.00000000.h5
--> Loading: data_order_1_U_05.00000000.h5
--> Loading: data_order_1_U_06.00000000.h5
--> Loading: data_order_2_U_03.00000000.h5
--> Loading: data_order_2_U_04.00000000.h5
--> Loading: data_order_2_U_05.00000000.h5
--> Loading: data_order_2_U_06.00000000.h5

Visualization

To visualize the metal insulator cross over we plot the double occupancy \(\langle n_\uparrow n_\downarrow \rangle\) as a function of \(U\) for NCA and OCA.

[7]:

plt.plot(nca.U, nca.d, 'o-', label='NCA')
plt.plot(oca.U, oca.d, 'o-', label='OCA')

plt.ylabel(r'$\langle n_{\uparrow} n_{\downarrow} \rangle$')
plt.xlabel('U')
plt.ylim(bottom=0, top=0.12)
plt.legend(loc='best');

../_images/tutorials_Metal_insulator_transition_13_0.svg

This reproduces the results in Fig. 4a in M. Eckstein, P. Werner, Phys. Rev. B 82, 115115 (2010).