# Multi-orbital models¶

In pybinding, if an onsite or hopping energy term is defined as a matrix (instead of a scalar), we refer to the resulting model as multi-orbital. The elements of the matrix term may correspond to different spins, electrons and holes, or any other degrees of freedom. These can have different physical meaning depending on the intend of the model. Because we’re talking in generic terms here, we’ll use orbital as a blanket term to refer to any degree of freedom, i.e. matrix element of an onsite or hopping term.

This section describes how these models can be defined and how the presence of multiple orbitals affects modifier functions and the results obtained from solvers. In general, it is as simple as replacing a scalar value with a matrix while all of the principals described in the Tutorial still apply.

Download this page as a Jupyter notebook

## Onsite and hopping matrices¶

Starting from the very beginning, the orbital count of a site is determined by the shape of the onsite energy matrix. Let’s take a look at a few possibilities:

lat = pb.Lattice([1, 0], [0, 1])
("A", [0.0, 0.0], 0.5),          # single-orbital: scalar
("B", [0.0, 0.2], [[1.5,  2j],   # two-orbital: 2x2 Hermitian matrix
[-2j, 1.5]]),
("C", [0.3, 0.1], np.zeros(2)),  # two-orbital: zero onsite term
("D", [0.1, 0.0], [[4, 0, 0],    # three-orbital: only diagonal
[0, 5, 0],
[0, 0, 6]]),
("E", [0.2, 0.2], [4, 5, 6])     # three-orbital: only diagonal, terse notation
)


The onsite term is required to be a square Hermitian matrix. If a 1D array is given instead of a matrix, it will be interpreted as the main diagonal of a square matrix (see sublattices D and E which have identical onsite term specified with different notations).

As seen above, sublattices don’t need to all have the same orbital count. The only thing to keep in mind is that the hopping matrix which connect a pair of sublattice sites must have the appropriate shape: the number of rows must match the orbital count of the source sublattice and the number of columns must match the destination sublattice.

lat.add_hoppings(
([0, 1], "A", "A", 1.2),          # scalar
([0, 1], "B", "B", [[1, 2],       # 2x2
[3, 4]]),
([0, 0], "B", "C", [[2j, 0],      # 2x2
[1j, 0]]),
([0, 0], "A", "D", [[1, 2, 3]]),  # 1x3
([0, 1], "D", "A", [[7],          # 3x1
[8],
[9]]),
([0, 0], "B", "D", [[1j, 0,  0],  # 2x3
[2,  0, 3j]])
)


If a matrix of the wrong shape is given, an informative error is raised:

>>> lat.add_one_hopping([0, 0], "A", "B", 0.6)
RuntimeError: Hopping size mismatch: from 'A' (1) to 'B' (2) with matrix (1, 1)
>>> lat.add_one_hopping([0, 1], "D", "D", [[1, 2, 3],
...                                        [4, 5, 6]])
RuntimeError: Hopping size mismatch: from 'D' (3) to 'D' (3) with matrix (2, 3)


After the Lattice is complete, a Model can be built as usual:

>>> model = pb.Model(lat, pb.primitive(2, 2))
>>> model.system.num_sites
20  # <-- 5 sites per unit cell and 2x2 cells: 5*2*2 == 20
>>> model.hamiltonian.shape
(44, 44)  # <-- 11 (1+2+2+3+3) orbitals per unit cell and 2x2 cells: 11*2*2 = 44


Sites refer to physical locations so their total count corresponds to the number of sublattices (A to E) multiplied by the number of times the unit cell is repeated. The Hamiltonian matrix is larger than num_sites due to the extra orbitals.

## Effect on modifier functions¶

The @onsite_energy_modifier and @hopping_energy_modifier functions work equally well for single- and multi-orbital models. In case of the latter, the energy argument of the modifiers will have a shape matching the onsite/hopping matrix term.

@pb.onsite_energy_modifier
def potential(energy, x):
"""Linear onsite potential as a function of x for a 2-orbital model"""
return energy + np.eye(2) * x


Note the np.eye(2) in the code above. The number 2 matches the 2-orbital structure of a specific model. Without this, energy + x would also add the value to the off-diagonal elements of the onsite matrix which is not desirable in this case.

The modifier defined above will only work for 2-orbital models. In general, we might want to create modifiers which work with any n-orbital model or with a mixed number of orbitals. For this we can use the sub_id modifier argument and its .eye attribute which supplies the correct matrix shape for any sublattice:

@pb.onsite_energy_modifier
def potential(energy, x, sub_id):
"""Same as above, but works for any n-orbital model"""
return energy + sub_id.eye * x


Even more generally, if we wish to apply completely different functions to the various sublattices, the sub_id argument can be used to create different branches in the modifier:

@pb.onsite_energy_modifier
def potential(energy, x, sub_id):
"""Applies different functions to different sublattices"""
if sub_id == "A":
return energy + x  # we know sublattice A is single-orbital
elif sub_id == "D":
energy[x > 0] += sub_id.eye * x  # the notation can be mixed with numpy indexing
return energy                    # apply only to sites where x > 0
elif sub_id == "B":
sigma_y = np.array([[0, -1j],
[1j,  0]])
return energy + sigma_y * 1.3 - np.eye(2) * 0.6  # add multiple 2x2 matrices
else:
return energy  # leave the other sublattices unchanged


This branching behavior is only supported by the sub_id and hop_id arguments. Do not try to create branches like this using any of the other modifier arguments:

"""Creating a position-dependent potential"""
# This is an error with anything except sub_id or hop_id
if x > 0:
return energy + 1
else:
return energy - 1

# Use this notation instead
energy[x >  0] += 1
energy[x <= 0] -= 1


On the other hand, sub_id and hop_id can be used with either of these variants with just a single caveat:

"""Sublattice-dependent potential"""
# This always works with sub_id and hop_id
if sub_id == "A":
return energy + 1
else:
return energy - 1

# This only works when all sublattices have the same number of orbitals,
# but it will raise an error for mixed orbital counts.
energy[sub_id == "A"] += 1
energy[sub_id == "B"] -= 1


## Local properties and plotting¶

When examining the local properties of a multi-orbital model, it is important to make the distinction between system indices which correspond to sites (unique positions) and Hamiltonian indices which correspond to the onsite or hopping terms in the Hamiltonian.

As shown in one of the previous examples, the number of sites in a system does not have to be equal to the size of the Hamiltonian matrix (hamiltonian.shape[0] >= num_sites). This affects how the system and Hamiltonian are indexed. System indices are always scalars and point to a single site position. For single-orbital models there is a 1:1 correspondence between system and Hamiltonian indices. However, for multi-orbital models the Hamiltonian indices are 1D arrays with a size corresponding to the number of orbitals on the target site.

>>> model = pb.Model(lat, pb.primitive(2, 2))
>>> sys_idx = model.system.find_nearest(position=[0, 0], sublattice="D")
>>> sys_idx  # <-- Points to a site on sublattice D which is closest to the target position.
15           #     It's always a scalar.
>>> model.system.x[sys_idx]
0.1  # <-- Not exactly 0 as requested, but the closest site to it.
>>> model.system.y[sys_idx]
0.0
>>> ham_idx = model.system.to_hamiltonian_indices(sys_idx)
>>> ham_idx   # <-- Array of integers which can be used to index the Hamiltonian matrix.
[29, 30, 31]  #     Size 3 because the selected site is on the 3-orbital sublattice D.
>>> ham = model.hamiltonian.todense()
>>> ham[np.ix_(ham_idx, ham_idx)]  # Returns the onsite hopping term of sublattice D.
[[4, 0, 0],
[0, 5, 0],
[0, 0, 6]]


Functions which compute various local properties take into account the presence of multiple orbitals on a single site. For example, when calculating the local density of states, one of the input parameters is the target site position. By default, the resulting LDOS is calculated as the sum of all orbitals but this is optional as shown in the following example:

"""Calculate the LDOS in the center of a MoS2 quantum dot"""
from pybinding.repository import group6_tmd

model = pb.Model(group6_tmd.monolayer_3band("MoS2"),
pb.regular_polygon(6, 20))

kpm = pb.kpm(model)
energy = np.linspace(-1, 3.8, 500)