# 6. Fields and effects¶

This section will introduce @onsite_energy_modifier and @hopping_energy_modifier which can be used to add various fields to the model. These functions can apply user-defined modifications to the Hamiltonian matrix which is why we shall refer to them as modifier functions.

Download this page as a Jupyter notebook

## 6.1. Electric potential¶

We can define a simple potential function like the following:

@pb.onsite_energy_modifier
def potential(x, y):
return np.sin(x)**2 + np.cos(y)**2


Here potential is just a regular Python function, but we attached a pretty @ decorator to it. The @onsite_energy_modifier decorator gives an ordinary function a few extra properties which we’ll talk about later. For now, just keep in mind that this is required to mark a function as a modifier for use with pybinding models. The x and y arguments are lattice site positions and the return value is the desired potential. Note the use of np.sin instead of math.sin. The x and y coordinates are numpy arrays, not individual numbers. This is true for all modifier arguments in pybinding. When you write modifier functions, make sure to always use numpy operations which work with arrays, unlike regular math.

Note

Modifier arguments are passed as arrays for performance. Working with individual numbers would require calling the potential function individually for each lattice site which would be extremely slow. Arrays are much faster.

To use the potential function, just place it in a Model parameter list:

from pybinding.repository import graphene

model = pb.Model(
graphene.monolayer(),
pb.rectangle(12),
potential
)


To visualize the potential, there’s the handy Model.onsite_map property which is a StructureMap of the onsite energy of the Hamiltonian matrix.

model.onsite_map.plot_contourf()
pb.pltutils.colorbar(label="U (eV)")


The figure shows a 2D colormap representation of our wavy potential in a square system. The StructureMap.plot_contourf() method we just called is implemented in terms of matplotlib’s contourf function with some slight adjustments for convenience.

To make the potential more flexible, it’s a good idea to enclose it in an outer function, just like this:

def wavy(a, b):
@pb.onsite_energy_modifier
def potential(x, y):
return np.sin(a * x)**2 + np.cos(b * y)**2
return potential

model = pb.Model(
graphene.monolayer(),
wavy(a=0.6, b=0.9)
)
model.onsite_map.plot_contourf()
pb.pltutils.colorbar(label="U (eV)")


Note that we are using a system with hexagonal shape this time (via regular_polygon()). The potential is only plotted inside the area of the actual system.

We can make one more improvement to our wavy function. We’ll add an energy argument:

def wavy2(a, b):
@pb.onsite_energy_modifier
def potential(energy, x, y):
v = np.sin(a * x)**2 + np.cos(b * y)**2
return energy + v
return potential


The energy argument contains the existing onsite energy in the system before the new potential function is applied. By adding to the existing energy, instead of just setting it, we can compose multiple functions. For example, let’s combine the improved wavy2 with a linear potential.

def linear(k):
@pb.onsite_energy_modifier
def potential(energy, x):
return energy + k*x
return potential

model = pb.Model(
graphene.monolayer(),
wavy2(a=0.6, b=0.9),
linear(k=0.2)
)
model.onsite_map.plot_contourf()
pb.pltutils.colorbar(label="U (eV)")


We see a similar wavy pattern as before, but the magnitude increases linearly along the x-axis because of the contribution of the linear potential.

Now that you have a general idea of how to add and compose electric potentials in a model, we should talk about the role of the @onsite_energy_modifier. The full signature of a potential function looks like this:

@pb.onsite_energy_modifier
def potential(energy, x, y, z, sub_id):
return ... # some function of the arguments


This function uses all of the possible arguments of an onsite energy modifier: energy, x, y, z and sub_id. We have already explained the first three. The z argument is, obviously, the z-axis coordinate of the lattice sites. The sub_id argument tells us which sublattice a site belongs to. Its usage will be explained below.

As we have seen before, we don’t actually need to define a function to take all the arguments. They are optional. The @ decorator will recognize a function which takes any of these arguments and it will adapt it for use in a pybinding model. Previously, the linear function accepted only the energy and x arguments, but wavy also included the y argument. The order of arguments is not important, only their names are. Therefore, this is also a valid modifier:

@pb.onsite_energy_modifier
def potential(x, y, energy, sub_id):
return ... # some function


But the argument names must be exact: a typo or an extra unknown argument will result in an error. The decorator checks this at definition time and decides if the given function is a valid modifier or not, so any errors will be caught early.

## 6.3. Opening a band gap¶

The last thing to explain about @onsite_energy_modifier is the use of the sub_id argument. It tells us which sublattice a site belongs to. If you remember from early on in the tutorial, in the process of specifying a lattice, we gave each sublattice a unique name. This name can be used to filter out sites of a specific sublattice. For example, let’s add mass to electrons in graphene:

def mass_term(delta):
"""Break sublattice symmetry with opposite A and B onsite energy"""
@pb.onsite_energy_modifier
def potential(energy, sub_id):
energy[sub_id == 'A'] += delta
energy[sub_id == 'B'] -= delta
return energy
return potential


Note that we don’t need x, y or z arguments because this will be applied everywhere evenly. The mass_term function will add an energy delta to all sites on sublattice A and subtract delta from all B sites. Note that we are indexing the energy array with a condition on the sub_id array of the same length. This is a standard numpy indexing technique which you should be familiar with.

The simplest way to demonstrate our new mass_term is with a graphene nanoribbon. First, let’s just remind ourselves what a pristine zigzag nanoribbon looks like:

model = pb.Model(
graphene.monolayer(),
pb.rectangle(1.2),
pb.translational_symmetry(a1=True, a2=False)
)
model.plot()


And let’s see its band structure:

from math import pi, sqrt

solver = pb.solver.lapack(model)
a = graphene.a_cc * sqrt(3)
bands = solver.calc_bands(-pi/a, pi/a)
bands.plot()


Note that the bands touch at zero energy: there is not band gap. Now, let’s include the mass term and compute the band structure again.

model = pb.Model(
graphene.monolayer(),
pb.rectangle(1.2),
pb.translational_symmetry(a1=True, a2=False),
mass_term(delta=2.5)  # eV
)
solver = pb.solver.lapack(model)
bands = solver.calc_bands(-pi/a, pi/a)
bands.plot()


We set a very high delta value of 2.5 eV for illustration purposes. Indeed, a band gap of 5 eV (delta * 2) is quite clearly visible in the band structure.

## 6.4. PN junction¶

While we’re working with a nanoribbon, let’s add a PN junction along its main axis.

def pn_junction(y0, v1, v2):
@pb.onsite_energy_modifier
def potential(energy, y):
energy[y < y0] += v1
energy[y >= y0] += v2
return energy
return potential


The y0 argument is the position of the junction, while v1 and v2 are the values of the potential (in eV) before and after the junction. Let’s add it to the nanoribbon:

model = pb.Model(
graphene.monolayer(),
pb.rectangle(1.2),
pb.translational_symmetry(a1=True, a2=False),
pn_junction(y0=0, v1=-5, v2=5)
)
pb.pltutils.colorbar(label="U (eV)")


Remember that the Model.onsite_map property is a StructureMap, which has several plotting methods. A contour plot would not look at all good for such a small nanoribbon, but the method StructureMap.plot() is perfect. As before, the ribbon has infinite length along the x-axis and the transparent sites represent the periodic boundaries. The PN junction splits the ribbon in half along its main axis.

We can compute and plot the band structure:

solver = pb.solver.lapack(model)
bands = solver.calc_bands(-pi/a, pi/a)
bands.plot()


Next, let’s create a square potential well. We could define a new modifier function, as before. But lets take a different approach and create the well by composing two PN junctions.

model = pb.Model(
graphene.monolayer(),
pb.rectangle(1.2),
pb.translational_symmetry(a1=True, a2=False),
pn_junction(y0=-0.2, v1=5, v2=0),
pn_junction(y0=0.2, v1=0, v2=5)
)
pb.pltutils.colorbar(label="U (eV)")


It works as expected. This can sometimes be a nice and quick way to extend a model. The square well affects the band structure by breaking electron-hole symmetry:

solver = pb.solver.lapack(model)
bands = solver.calc_bands(-pi/a, pi/a)
bands.plot()


## 6.5. Magnetic field¶

To model a magnetic field, we need to apply the Peierls substitution:

$t_{nm} \rightarrow t_{nm} \text{e}^{i\frac{2\pi}{\Phi_0} \int_n^m \vec{A}_{nm} \cdot d\vec{l}}$

Here $$t_{nm}$$ is the hopping energy between two sites, $$\Phi_0 = h/e$$ is the magnetic quantum, $$h$$ is the Planck constant and $$\vec{A}_{nm}$$ is the magnetic vector potential along the path between sites $$n$$ and $$m$$. We want the magnetic field to be perpendicular to the graphene plane, so we can take the gauge $$\vec{A}(x,y,z) = (By, 0, 0)$$.

This can all be expressed with a @hopping_energy_modifier:

from pybinding.constants import phi0

def constant_magnetic_field(B):
@pb.hopping_energy_modifier
def function(energy, x1, y1, x2, y2):
# the midpoint between two sites
y = 0.5 * (y1 + y2)
# scale from nanometers to meters
y *= 1e-9

# vector potential along the x-axis
A_x = B * y

# integral of (A * dl) from position 1 to position 2
peierls = A_x * (x1 - x2)
# scale from nanometers to meters (because of x1 and x2)
peierls *= 1e-9

# the Peierls substitution
return energy * np.exp(1j * 2*pi/phi0 * peierls)
return function


The energy argument is the existing hopping energy between two sites at coordinates (x1, y1) and (x2, y2). The function computes and returns the Peierls substitution as given by the equation above.

The full signature of a @hopping_energy_modifier is actually:

@pb.hopping_energy_modifier
def function(energy, x1, y1, z1, x2, y2, z2, hop_id):
return ... # some function of the arguments


The hop_id argument tells us which type of hopping it is. Hopping types can be specifically named during the creation of a lattice. This can be used to apply functions only to specific hoppings. However, as with all the modifier arguments, it’s optional, so we only take what we need.

To test out our constant_magnetic_field, we’ll calculate the local density of states (LDOS), where we expect to see peaks corresponding to Landau levels. The computation method used here is explained in detail in the Kernel polynomial method section of the tutorial.

model = pb.Model(
graphene.monolayer(),
pb.rectangle(30),
constant_magnetic_field(B=200)  # Tesla
)
kpm = pb.kpm(model)

ldos = kpm.calc_ldos(energy=np.linspace(-1, 1, 500), broadening=0.015, position=[0, 0])
ldos.plot()
plt.show()


The values of the magnetic field is exaggerated here (200 Tesla), but that is done to keep the computation time low for the tutorial (less than 0.5 seconds for this LDOS calculation).

## 6.7. Example¶

Download source code

"""PN junction and broken sublattice symmetry in a graphene nanoribbon"""
import pybinding as pb
import matplotlib.pyplot as plt
from pybinding.repository import graphene
from math import pi, sqrt

pb.pltutils.use_style()

def mass_term(delta):
"""Break sublattice symmetry with opposite A and B onsite energy"""
@pb.onsite_energy_modifier
def potential(energy, sub_id):
energy[sub_id == 'A'] += delta
energy[sub_id == 'B'] -= delta
return energy

return potential

def pn_juction(y0, v1, v2):
"""PN junction potential

The y0 argument is the position of the junction, while v1 and v2
are the values of the potential (in eV) before and after the junction.
"""
@pb.onsite_energy_modifier
def potential(energy, y):
energy[y < y0] += v1
energy[y >= y0] += v2
return energy

return potential

model = pb.Model(
graphene.monolayer(),
pb.rectangle(1.2),  # width in nanometers
pb.translational_symmetry(a1=True, a2=False),
mass_term(delta=2.5),  # eV
pn_juction(y0=0, v1=-2.5, v2=2.5)  # y0 in [nm] and v1, v2 in [eV]
)
model.plot()
plt.show()

# plot the potential: note that pn_junction cancels out delta on some sites

# compute the bands