.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/example_06_Hubbard_dimer.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_example_06_Hubbard_dimer.py: Hubbard Dimer ===================================== This exercise will demonstrate how to handle hopping and multi-site problems within edrixs using the example of a Hubbard dimer. We want to solve the equation .. math:: \hat{H} = \sum_{i,j} \sum_{\sigma} t_{i,j} \hat{f}^{\dagger}_{i,\sigma} \hat{f}_{j, \sigma} + U \sum_{i} \hat{n}_{i,\uparrow}\hat{n}_{i,\downarrow}, which involves two sites labeled with indices :math:`i` or :math:`j` with two electrons of spin :math:`\sigma\in{\uparrow,\downarrow}`. :math:`t_{i,j}` is the hopping between sites, :math:`\hat{f}^{\dagger}_{i,\sigma}` is the creation operators, and :math:`\hat{n}^{\dagger}_{i,\sigma}=\hat{f}^{\dagger}_{i,\sigma}\hat{f}_{i,\sigma}` is the number operator. The main task is to represent this Hamiltonian and the related spin operator using the EDRIXS two-fermion and four-fermion form where :math:`\alpha,\beta,\delta,\gamma` are the indices of the single particle basis. .. math:: \hat{H} = \sum_{\alpha,\beta} t_{\alpha,\beta} \hat{f}^{\dagger}_{\alpha} \hat{f}_{\beta} + \sum_{\alpha,\beta,\gamma,\delta} U_{\alpha,\beta,\gamma,\delta} \hat{f}^{\dagger}_{\alpha}\hat{f}^{\dagger}_{\beta}\hat{f}_{\gamma}\hat{f}_{\delta}. .. GENERATED FROM PYTHON SOURCE LINES 30-35 Initialize matrices ------------------------------------------------------------------------------ We start by noting that each of the two sites is like an :math:`l=0` :math:`s`-orbital with two spin-orbitals each. We will include two electron occupation and build the Fock basis. .. GENERATED FROM PYTHON SOURCE LINES 35-47 .. code-block:: Python import numpy as np import matplotlib.pyplot as plt import edrixs np.set_printoptions(precision=4) ll = 0 case = 's' norb = 4 noccu = 2 basis = edrixs.get_fock_basis_int(norb, noccu) .. GENERATED FROM PYTHON SOURCE LINES 48-67 Create function to populate and diagonalize matrices ------------------------------------------------------------------------------ The Coulomb and hopping matrices :code:`umat` and :code:`emat` will be represented by :math:`4\times4\times4\times4` and :math:`4\times4` matrices, respectively. Note that we needed to specify that these are, in general, complex, although they happen to contain only real numbers in this case. We follow the convention that these are ordered first by site and then by spin: :math:`|0\uparrow>, |0\downarrow>, |1\uparrow>, |1\downarrow>`. Consequently the :math:`2\times2` and :math:`2\times2\times2\times2` block diagonal structures of the matrices will contain the on-site interactions. The converse is true for the hopping between the sites. From here let us generate a function to build and diagonalize the Hamiltonian. We need to generate the Coulomb matrix for the on-site interactions and apply it to the block diagonal. The hopping connects off-site indices with the same spin. :func:`~edrixs.solvers.build_op` transforms the one-body matrix :code:`emat` and the two-body tensor :code:`umat` into the many-body Hamiltonian in our Fock basis and :func:`~edrixs.solvers.ed` diagonalizes it. This problem is tiny, so we use the :code:`'dense'` backend. .. GENERATED FROM PYTHON SOURCE LINES 67-82 .. code-block:: Python def diagonalize(U, t, extra_emat=None): """Diagonalize 2 site Hubbard Hamiltonian""" umat = np.zeros((norb, norb, norb, norb), dtype=np.complex128) emat = np.zeros((norb, norb), dtype=np.complex128) U_mat_1site = edrixs.get_umat_slater('s', U) umat[:2, :2, :2, :2,] = umat[2:, 2:, 2:, 2:] = U_mat_1site emat[2, 0] = emat[3, 1] = emat[0, 2] = emat[1, 3] = t if extra_emat is not None: emat = emat + extra_emat H = edrixs.build_op(emat, umat, basis, backend='dense') e, v = edrixs.ed(H, num_evals=len(basis)) return e, v .. GENERATED FROM PYTHON SOURCE LINES 83-86 The large :math:`U` limit ------------------------------------------------------------------------------ Let us see what happens with :math:`U \gg t`. .. GENERATED FROM PYTHON SOURCE LINES 86-90 .. code-block:: Python e, v = diagonalize(1000, 1) print("Energies are") print(e) .. rst-class:: sphx-glr-script-out .. code-block:: none Energies are [ -0.004 0. 0. 0. 1000. 1000.004] .. GENERATED FROM PYTHON SOURCE LINES 91-94 To analyze what is going on we can determine the spin expectation values of the cluster. Building the operators follows the same form as the Hamiltonian and the previous example. .. GENERATED FROM PYTHON SOURCE LINES 94-103 .. code-block:: Python spin_mom_one_site = edrixs.get_spin_momentum(ll) spin_mom = np.zeros((3, norb, norb), dtype=np.complex128) spin_mom[:, :2, :2] = spin_mom[:, 2:, 2:] = spin_mom_one_site opS = [edrixs.build_op(component, None, basis, backend='dense')() for component in spin_mom] opS_squared = (np.dot(opS[0], opS[0]) + np.dot(opS[1], opS[1]) + np.dot(opS[2], opS[2])) .. GENERATED FROM PYTHON SOURCE LINES 104-107 This time let us include a tiny magnetic field along the :math:`z`-axis, so that we have a well-defined measurement axis and print out the expectation values. .. GENERATED FROM PYTHON SOURCE LINES 107-119 .. code-block:: Python zeeman = np.zeros((norb, norb), dtype=np.complex128) zeeman[:2, :2] = zeeman[2:, 2:] = 1e-8*spin_mom_one_site[2] e, v = diagonalize(1000, 1, extra_emat=zeeman) Ssq_exp = edrixs.cb_op(opS_squared, v).diagonal().real Sz_exp = edrixs.cb_op(opS[2], v).diagonal().real header = "{:<10s}\t{:<6s}\t{:<6s}" print(header.format("E", "S(S+1)", "")) for i in range(len(e)): print("{:<2f}\t{:.1f}\t{:.1f}".format(e[i], Ssq_exp[i], Sz_exp[i])) .. rst-class:: sphx-glr-script-out .. code-block:: none E S(S+1) -0.004000 0.0 0.0 -0.000000 2.0 -1.0 -0.000000 2.0 0.0 0.000000 2.0 1.0 1000.000000 0.0 0.0 1000.004000 0.0 0.0 .. GENERATED FROM PYTHON SOURCE LINES 120-124 For :math:`U \gg t` the two states with double occupancy acquire an energy of approximately :math:`U`. The low energy states are a :math:`S=0` singlet and an :math:`S=1` triplet, which are split by :math:`4t^2/U`, the magnetic exchange energy. .. GENERATED FROM PYTHON SOURCE LINES 126-129 :math:`U` dependence ------------------------------------------------------------------------------ Let us plot the changes in energy with :math:`U`. .. GENERATED FROM PYTHON SOURCE LINES 129-140 .. code-block:: Python plt.figure() t = 1 Us = np.linspace(0.01, 10, 50) Es = np.array([diagonalize(U, t, extra_emat=zeeman)[0] for U in Us]) plt.plot(Us/t, Es/t) plt.xlabel('U/t') plt.ylabel('Eigenstate energies/t') plt.show() .. image-sg:: /auto_examples/images/sphx_glr_example_06_Hubbard_dimer_001.png :alt: example 06 Hubbard dimer :srcset: /auto_examples/images/sphx_glr_example_06_Hubbard_dimer_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 141-143 To help interpret this, we can represent the eigenvectors in terms of a sum of the single particle states. .. GENERATED FROM PYTHON SOURCE LINES 143-163 .. code-block:: Python def get_single_particle_repesentations(v): # binary occupation vectors, in the same order as ``basis`` fock_vectors = np.array(edrixs.get_fock_bin_by_N(norb, noccu)) reps = [] for i in range(6): rep = sum([vec*weight for weight, vec in zip(v[:, i], fock_vectors)]) reps.append(rep) return np.array(reps) t = 1 for U in [10000, 0.0001]: e, v = diagonalize(U, t, extra_emat=zeeman) repesentations = get_single_particle_repesentations(v) print("For U={} t={} states are".format(U, t)) print(repesentations.round(3).real) print("\n") .. rst-class:: sphx-glr-script-out .. code-block:: none For U=10000 t=1 states are [[-0.707 0.707 0.707 -0.707] [ 0. 1. 0. 1. ] [ 0.707 0.707 0.707 0.707] [-1. 0. -1. 0. ] [-0.707 -0.707 0.707 0.707] [ 0.707 0.707 0.707 0.707]] For U=0.0001 t=1 states are [[-0. 1. 1. -0. ] [ 0. 1. 0. 1. ] [ 0.707 0.707 0.707 0.707] [-1. 0. -1. 0. ] [-0.707 -0.707 0.707 0.707] [-1. -0. -0. -1. ]] .. GENERATED FROM PYTHON SOURCE LINES 164-169 For :math:`U \gg t` the ground state maximizes its magnetic exchange energy saving. In the :math:`U \ll t` condition the ground state maximizes its kinetic energy saving. Since both states share the same parity, the cross-over between them is smooth. This type of physics is at play in current research on quantum materials [1]_ [2]_. .. GENERATED FROM PYTHON SOURCE LINES 171-175 .. rubric:: Footnotes .. [1] Y. Wang et al., `Phys. Rev. Lett. 122, 106401 (2019) `_. .. [2] A. Revelli et al., `Science Advances 5, eaav4020 (2019) `_. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.064 seconds) .. _sphx_glr_download_auto_examples_example_06_Hubbard_dimer.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_06_Hubbard_dimer.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_06_Hubbard_dimer.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example_06_Hubbard_dimer.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_