How to model quantum wells

This tutorial shows you how to use SigSpace to study one of the most commonly used model systems: a 1D quantum well. You will learn how to create Hamiltonian operators, solve them, and reproduce key properties of quantum wells.

Define the system to solve

The first step is to define the system we want to study. This information is stored in the Hamiltonian. Therefore, we need to define the Hamiltonian in order to solve the problem:

In [1]:
import sigspace as sg
H = sg.Hamiltonian(sdims=1)
H += 1.23
print(H)
-0.5*Δ + 1.23

The first line imports the SigSpace library. Next, we define an empty Hamiltonian. The argument sdims specifies the number of spatial dimensions (in this case, we are working with a 1D problem).

You can treat the Hamiltonian like a mathematical object, adding terms as needed. For example, we could add a constant term, as shown in the third line. Later, we'll see how to add more complex terms to model more interesting problems.

The final line prints a summary of the Hamiltonian, which looks like this:

-0.5*Δ + 1.23 Eh

Every Hamilton operator always includes the kinetic energy term (-0.5*Δ), which is the default and cannot be modified. The second term, 1.23 Eh, represents the constant we added earlier. “Eh” is the symbol for Hartree, the unit of energy in atomic units. Internally, all computations are performed in atomic units to ensure a good numerical stability.

Let's calculate something

Now, let’s define a new Hamiltonian with a custom potential and solve it:

In [2]:
import sigspace as sg
import numpy as np

def Vfun(x):
    return -1/(1+x ** 2)

H = sg.Hamiltonian(sdims=1) + Vfun

solver = sg.Solver(H)
solver.grid(-10, 10, 200)
sol = solver.run()

sg.quickplot(sol)

Let’s walk through the code:

  • Custom-shaped potentials can be defined using callback functions. To achieve this, we define the function Vfun(x), which returns the potential at the location x. In this case, the function generates a negative (attractive) potential with its minimum at zero.
  • Next, we add this callback function to the Hamiltonian. If we print the Hamiltonian, the summary will show that a custom potential has been added. You can add multiple callback functions, and they will be summed up.
  • We then create a solver object using solver = sg.Solver(H). SigSpace provides several solver implementations, and sg.Solver(H) automatically selects the most appropriate one based on the Hamiltonian, returning a new solver instance.
  • The command solver.grid(-10, 10, 200) creates a regular grid from -10 to 10 with 200 points. In some cases (e.g., when specifying atoms), the solver attempts to guess a useful grid. However, when working with custom potentials a grid has to be specified manually.
  • Now, we're ready to solve the problem. The command sol = solver.run() performs the calculation and stores the results in sol. The energy of the wave function is accessible via sol.energy, and the wave function itself is stored in sol.wave.
  • For quick visualization of the results, you can use quickplot. This function automatically selects the most suitable graphical representation based on the problem being solved.

Quantum Wells

Now that we’ve learned how to solve problems with custom potentials, we can apply this knowledge to explore the properties of quantum wells. To streamline the process, we can define a function that handles all the main computations for us:

In [3]:
import sigspace as sg
import numpy as np

def compute_Q_well(wells, offset=0):
    def V_well(x, well):
        return well["height"] * (np.abs(x - well["pos"]) < well["width"] / 2)

    H = sg.Hamiltonian(sdims=1) + offset
    for well in wells:
        H += lambda x, well=well: V_well(x, well)

    solver = sg.Solver(H)
    solver.grid(-10, 10, 200)
    solver.print_info = True
    return solver.run()

The compute_Q_well function is quite similar to the code from the previous example. It takes a list of dictionaries as an argument, where each entry should define the position, height, and width of a quantum well. As before, we define a Hamiltonian and add as many potentials as there are wells in the list.

To pass additional arguments to the potential function, we use a lambda statement: lambda x, well=well: V_well(x, well). It's important to include well in the lambda function’s arguments; otherwise, only the last argument value will be used.

The second argument, offset, is a constant that can be added to the potential.

By setting solver.print_info = True, we instruct the solver to print a brief summary when the computation is completed.

Now, let's test it with a simple problem:

In [4]:
wells = [{"height": -1, "pos": 0, "width": 1}]

sol = compute_Q_well(wells, offset=1)

x = sol.wave.x
sg.quickplot((x, sol.wave()), (x, sol.H.potential(x)))
Iter: 2488, energy 0.687, residual 0.009989, duration 0.002943s

And it works as expected: we’ve defined a positive potential (shown above) that reaches zero at the center. The wave function peaks inside this well, where the potential is at its minimum.

Additionally, the solver provides a concise summary of the computation, including:

  • the number of iterations needed to solve the problem
  • the resulting energy
  • the remaining error (the default threshold is 1e-2, but this can be adjusted in the solver settings)
  • the duration of the computation

When passing tuples of x- and y-coordinates to quickplot, it will plot them together in one figure. This is especially useful for comparing wave functions when parameters change.

Note that the offset is not essential for showing the general shape of the wave function. If we were to remove it:

In [5]:
sol2 = compute_Q_well(wells, offset = 0)

print(sol.energy)
print(sol2.energy+1)
Iter: 2488, energy -0.313, residual 0.009989, duration 0.002931s
0.6870869029958183
0.6870869029958208

we would still obtain the same result by simply adding the offset to the resulting energy. The wave functions are identical, as constant values in the Hamiltonian just add to the energy.

A common observation in the calculation above is that the wave function is not limited strictly to the quantum well. It extends in the forbidden (repulsive) region. This indicates that, although with reduced probability, electrons can still be found in this region.

The extent to which the wave function penetrates into the repulsive area can be adjusted by altering the depth of the quantum well:

In [6]:
wells = [{"height": -1, "pos": 0, "width": 1}]
sol1 = compute_Q_well(wells)

wells = [{"height": -2, "pos": 0, "width": 1}]
sol2 = compute_Q_well(wells)

wells = [{"height": -10, "pos": 0, "width": 1}]
sol3 = compute_Q_well(wells)

x = sol1.wave.x

sg.quickplot((x, sol1.wave()), (x, sol2.wave()), (x, sol3.wave()), (x, sol1.H.potential(x)))
Iter: 2488, energy -0.313, residual 0.009989, duration 0.002912s
Iter: 1728, energy -0.912, residual 0.009948, duration 0.002095s
Iter: 690, energy -7.738, residual 0.009698, duration 0.000877s

The deeper the quantum well or the higher the energy barrier, the smaller the wave function's extension into the repulsive region (|x| > 0.5). A similar effect occurs when the width of the quantum well is adjusted.

This characteristic of the wave function is crucial for a phenomenon known as quantum tunneling. It allows particles to pass through energy barriers, even when they lack the energy required to overcome them classically.

Now, let’s consider a system with two quantum wells:

In [7]:
wells = [{"height": -1, "pos": -3, "width": 1},
         {"height": -2, "pos":  3, "width": 1}]


sol = compute_Q_well(wells)
sg.quickplot((x, sol.wave()), (x, sol.H.potential(x)))
Iter: 3726, energy -0.912, residual 0.009981, duration 0.004283s

In this simulation, the two quantum wells are relatively far apart. Since one well is deeper, the electron in the stationary state is most likely to be found there.

The wave function extends over the quantum well, but in the second one the probability is zero. However, if we bring the wells closer together, the wave function begins to change:

In [8]:
wells = [{"height": -1, "pos": -0.6, "width": 1},
         {"height": -2, "pos":  0.6, "width": 1}]


sol = compute_Q_well(wells)
sg.quickplot((x, sol.wave()), (x, sol.H.potential(x)))
Iter: 1844, energy -1.040, residual 0.009977, duration 0.002172s

Now, the wave function not only overlaps with the neighboring quantum well but also undergoes a change in shape.