Terminal currents

Supercurrent in a superscreen.Device can be classified into one of three categories:

1. Screening currents, which are the response of a film to an applied magnetic field.

2. Circulating currents, which are the response of a film to flux trapped in a hole.

3. Terminal currents, which flow into and out of a film through “transport terminals” defined on the film boundary.

In this notebook we demonstrate the treatment of terminal currents in a superscreen model.

Background

We make the following assumptions regarding terminal currents:

1. The current density \(|\vec{J}|\) is uniformly distributed along the terminal.

2. Along the terminal, the current direction \(\vec{J}/|\vec{J}|\) is perpendicular to the terminal direction.

The stream function \(g(\vec{r})=g(x, y)\) and 2D supercurrent density \(\vec{J}(x, y)\) are related according to \(\vec{J}=\vec{\nabla}\times(g\hat{z})=(\partial g/\partial y, -\partial g/\partial x)\). For some assumed supercurrent distribution \(\vec{J},\) the associated stream function \(g\) is given by a line integral,

\[g(\vec{r})=g(\vec{r}_0) + \int_{\vec{r}_0}^\vec{r}(\hat{z}\times\vec{J})\cdot\mathrm{d}\vec{\ell},\]

where \(\hat{z}\times\vec{J}=(-J_y, J_x)\) and \(\vec{r}_0\) is some reference position.

The stream function \(g\) in a film can be decomposed into \(g=g_\mathrm{screening}+g_\mathrm{transport},\) where \(g_\mathrm{screening}\) arises from the response to an applied magnetic field or trapped flux, and \(g_\mathrm{transport}\) corresponds to an applied current bias. Below we describe how to set \(g_\mathrm{transport}(\vec{r})\) for a given current bias configuration.

Stream function on the boundary

The boundary conditions for terminal currents are as follows:

• The boundary of the film is oriented in a counterclockwise direction.

• The stream function \(g_\mathrm{transport}(\vec{r})\) along a source terminal \(s\) injecting current \(I_s\) changes linearly by a total amount \(-I_s,\) and \(g_\mathrm{transport}(\vec{r})\) for \(\vec{r}\) along the terminal is given by the equation above, where \(\vec{r}_0\) is the start of the terminal.

• The stream function \(g_\mathrm{transport}(\vec{r})\) on the boundary between terminals is constant.

• The sum of all terminal currents must be zero to ensure current conservation.

Stream function in the bulk

Once \(g_\mathrm{transport}(\vec{r}_\text{boundary})\) has been defined for all \(\vec{r}_\text{boundary}\) on the film’s boundary, we solve for \(g_\mathrm{transport}(\vec{r}_\text{bulk})\) inside the film as follows.

1. First we find an “effective field” \(H_{z,\,\text{eff}}(\vec{r})\) that would generate such a \(g_\mathrm{transport}(\vec{r}_\text{boundary})\).

2. Once \(H_{z,\,\text{eff}}(\vec{r})\) is found, \(g_\mathrm{transport}(\vec{r}_\text{bulk})\) is given by the response of the film to this applied effective field.

Stream function in holes

For films with holes, we first solve for \(g_\mathrm{transport}(\vec{r})\) as described above assuming that there are no holes. Then, for each hole \(h\) we update \(g_\mathrm{transport}(\vec{r})\) for \(\vec{r}\) inside \(h\) to be equal to the average value over the hole area:

\[g_\mathrm{transport}(\vec{r})\to \frac{1}{\mathrm{area}(h)}\int_h g_\mathrm{transport}(\vec{r}')\,\mathrm{d}^2r'\quad\forall\:\vec{r}\text{ in hole }h\]

Examples

# Automatically install superscreen from GitHub only if running in Google Colab
if "google.colab" in str(get_ipython()):
    %pip install --quiet git+https://github.com/loganbvh/superscreen.git

%config InlineBackend.figure_formats = {"retina", "png"}
%matplotlib inline

import os

os.environ["OPENBLAS_NUM_THREADS"] = "1"

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams["figure.figsize"] = (5, 4)
plt.rcParams["font.size"] = 10

import superscreen as sc
from superscreen.geometry import box, circle, rotate

sc.version_table()

Software	Version
SuperScreen	0.13.0
Numpy	2.4.4
Numba	0.65.1
SciPy	1.17.1
matplotlib	3.10.9
IPython	9.13.0
Python	3.12.12 (main, Apr 27 2026, 17:32:11) [GCC 11.4.0]
OS	posix [linux]
Number of CPUs	Physical: 1, Logical: 2
BLAS Info	Generic
Sat May 02 12:32:50 2026 UTC

Film with a constriction

Define the device geometry.

The source terminal(s) and drain terminal are defined as instances of superscreen.Polygon. All points in the mesh that are inside a terminal polygon and on the boundary of the mesh will have the transport boundary conditions applied.

width = 1
height = width * 2
slot_height = height / 5
slot_width = width / 4
x0, y0 = center = (0, 0)
length_units = "um"

film = sc.Polygon("film", layer="base", points=box(width, height))
slot = sc.Polygon(
    points=box(slot_width, slot_height, center=(-(width - slot_width) / 2, 0))
)
film = film.difference(slot, slot.scale(xfact=-1)).resample(401)

source_terminal = sc.Polygon(
    "source", points=box(width, height / 100, center=(0, height / 2))
)
drain_terminal = source_terminal.translate(dy=-height).set_name("drain")

device = sc.Device(
    "constriction",
    layers=[sc.Layer("base", Lambda=0.1)],
    films=[film],
    terminals={film.name: [source_terminal, drain_terminal]},
    length_units=length_units,
)

device.make_mesh(min_points=3000, smooth=0)

fig, ax = device.plot_mesh(edge_color="k", show_sites=False)
_ = device.plot_polygons(ax=ax, legend=True)

Solve the device with a given source-drain current and no applied field.

terminal_currents = {"film": {"source": "100 uA", "drain": "-100 uA"}}

solution = sc.solve(
    device,
    applied_field=sc.sources.ConstantField(0),
    terminal_currents=terminal_currents,
    current_units="uA",
)[-1]

Plot the stream function.

fig, axes = solution.plot_streams()
_ = device.plot_polygons(ax=axes[0])
_ = axes[0].legend(loc="upper left", bbox_to_anchor=(1.5, 1))

Plot the supercurrent distribution.

xs = np.linspace(-(1.05 * width) / 2, (1.05 * width) / 2, 501)
ys = np.ones_like(xs)

sections = [
    np.stack([xs, 0.6 * ys], axis=1),
    np.stack([xs, 0.3 * ys], axis=1),
    np.stack([xs, 0.0 * ys], axis=1),
]

fig, axes = solution.plot_currents(
    streamplot=False,
    cross_section_coords=sections,
    figsize=(3, 6),
)
_ = device.plot_polygons(ax=axes[0], lw=2)

_ = solution.plot_fields()

Compute the total current flowing through the various cross-sections.

for coords in sections:
    total_current = solution.current_through_path(coords, film="film", units="uA").magnitude
    target_current = solution.terminal_currents["film"]["source"]
    err = (total_current - target_current) / abs(target_current) * 100
    print(
        f"Cross-section: y = {coords[0, 1]:.2f} um, total current = {total_current:.3f} uA"
        f" ({err:.2f}% error)"
    )

Cross-section: y = 0.60 um, total current = 100.007 uA (0.01% error)
Cross-section: y = 0.30 um, total current = 100.392 uA (0.39% error)
Cross-section: y = 0.00 um, total current = 99.892 uA (-0.11% error)

Film with holes

device.holes = {
    "hole1": sc.Polygon(
        "hole1",
        layer="base",
        points=circle(width / 4, center=(-width / 8, +height / 4)),
    ),
    "hole2": sc.Polygon(
        "hole2",
        layer="base",
        points=circle(width / 4, center=(+width / 8, -height / 4)),
    ),
}

_ = device.draw()

Find the zero-fluxoid solution for a given source-drain current.

Note that if the two holes were horizontally aligned in the center of the device (such that the device was mirror-symmetric about \(x=0\)), then by symmetry the zero-fluxoid solution would have no circulating current.

terminal_currents = {"film": {"source": "100 uA", "drain": "-100 uA"}}

model = sc.factorize_model(device=device, terminal_currents=terminal_currents, current_units="uA")

solution = sc.find_fluxoid_solution(
    model,
    applied_field=sc.sources.ConstantField(0),
)

fig, axes = solution.plot_streams(symmetric_color_scale=False, cmap="magma")
_ = device.plot_polygons(ax=axes[0], lw=2)

fig, axes = solution.plot_currents(streamplot=True)
_ = device.plot_polygons(ax=axes[0], lw=2)

Parallel SQUID array model

Here we model a device studied in detail in Theoretical model for parallel SQUID arrays with fluxoid focusing, Phys. Reb. B 103, 054509 (2021) (arXiv:2009.05338).

The device consists of a \(d=0.125\,\mu\mathrm{m}\) thick YBCO film with an estimated London penetration depth of \(\lambda=0.3\,\mu\mathrm{m}\). The film is pattered into a device with approximately the dimensions shown below. The device has source and drain leads at the top and bottom respectively, and 10 rectangular holes along the vertical midpoint.

Among other things, the authors simulated the current distribution for a total source-drain current of \(200\,\mu\mathrm{A}\) with both \(0\,\mu\mathrm{T}\) and \(12\,\mu\mathrm{T}\) out-of-plane applied field. Those results, plotted as current streamlines for the top half of the device, are shown in Figure 5 of the paper linked above (which is reproduced below).

Below we demonstrate how to reproduce these results using SuperScreen.

layer = sc.Layer("base", london_lambda=0.3, thickness=0.125)
film = sc.Polygon("film", layer="base", points=box(68, 24)).union(box(10, 40)).resample(401)
w_h, h = 4, 2
base_hole = sc.Polygon(layer="base", points=box(w_h, 2 * w_h)).resample(51)
holes = []
for i in range(10):
    hole = base_hole.translate(dx=1.5 * h + (i - 5) * w_h * 1.5, dy=0)
    hole.name = f"hole{i}"
    holes.append(hole)

source = sc.Polygon("source", points=box(10, 0.1, center=(0, 20)))
drain = sc.Polygon("drain", points=box(10, 0.1, center=(0, -20)))

device = sc.Device(
    "squid_array",
    films=[film],
    layers=[layer],
    holes=holes,
    terminals={film.name: [source, drain]},
    length_units="um",
)

fig, ax = device.draw()

device.make_mesh(max_edge_length=1.25, smooth=50)

fig, ax = device.plot_mesh()
_ = device.plot_polygons(ax=ax, color="k", lw=0.5)

# Region in which we will evaluate the field from the device.
eval_region = sc.Polygon(points=box(68, 40))
eval_mesh = eval_region.make_mesh(max_edge_length=1.5)

def solve_and_plot_model(source_current="200 uA", applied_field="0 uT"):

    # Solve the model
    terminal_currents = {"film": {"source": source_current, "drain": f"-{source_current}"}}
    applied_field = sc.ureg(applied_field).to("uT").magnitude
    model = sc.factorize_model(
        device=device,
        terminal_currents=terminal_currents,
        current_units="uA",
    )
    solution = sc.solve(
        model=model,
        applied_field=sc.sources.ConstantField(applied_field),
        field_units="uT",
    )[-1]

    # # Uncomment to find the zero-fluxoid solution.
    # solution = sc.find_fluxoid_solution(
    #     model,
    #     applied_field=sc.sources.ConstantField(applied_field),
    #     field_units="uT",
    # )

    # Define cross sections
    xs = np.linspace(-70 / 2, +70 / 2, 501)
    ys = np.ones_like(xs)
    sections = [
        np.array([xs, 15.0 * ys]).T,
        np.array([xs, 10.0 * ys]).T,
        np.array([xs, 7.5 * ys]).T,
        np.array([xs, 0.0 * ys]).T,
    ]

    # Plot currents.
    fig, axes = solution.plot_currents(
        streamplot=True,
        min_stream_amp=1e-4,
        cross_section_coords=sections,
        figsize=(5, 4),
    )
    _ = device.plot_polygons(ax=axes[0], lw=1, color="w", ls="--")

    # Evaluate the current through each cross section.
    for coords in sections:
        total_current = solution.current_through_path(coords, film="film", units="uA").magnitude
        target_current = solution.terminal_currents["film"]["source"]
        err = (total_current - target_current) / abs(target_current) * 100
        print(
            f"Cross-section: y = {coords[0, 1]:.2f} um, total current = {total_current:.3f} uA"
            f" ({err:.2f}% error)"
        )

    # Plot fields
    fig, ax = solution.plot_field_at_positions(
        eval_mesh,
        zs=0.75,
        cross_section_coords=sections,
        figsize=(5, 4),
    )
    _ = device.plot_polygons(ax=ax, lw=1, color="w", ls="--")

    return solution

Applied field = \(0\,\mu\mathrm{T}\)

solution_0uT = solve_and_plot_model(source_current="200 uA", applied_field="0 uT")

Cross-section: y = 15.00 um, total current = 199.409 uA (-0.30% error)
Cross-section: y = 10.00 um, total current = 200.009 uA (0.00% error)
Cross-section: y = 7.50 um, total current = 200.010 uA (0.01% error)
Cross-section: y = 0.00 um, total current = 202.896 uA (1.45% error)

Applied field = \(12\,\mu\mathrm{T}\)

solution_12uT = solve_and_plot_model(source_current="200 uA", applied_field="12 uT")

Cross-section: y = 15.00 um, total current = 199.487 uA (-0.26% error)
Cross-section: y = 10.00 um, total current = 200.162 uA (0.08% error)
Cross-section: y = 7.50 um, total current = 200.473 uA (0.24% error)
Cross-section: y = 0.00 um, total current = 200.932 uA (0.47% error)

Multiple source terminals

As mentioned above, while there can be aribitrarily many source terminals, there can technically only be a single drain terminal. However, it is possible to model systems with multiple current sinks by changing the the sign of terminal currents.

Below we demonstrate this with a “+”-shaped device with four current terminals.

layer = sc.Layer("base", Lambda=1)
width, height = 10, 2
points = box(width, height)
bar = sc.Polygon("plus", points=points)
plus = bar.union(bar.rotate(90)).resample(501)
plus.name = "plus"
plus.layer = layer.name
terminal = sc.Polygon(points=box(height, width / 50, center=(0, -width / 2)))

terminals = []
for i, name in enumerate(["drain", "source1", "source2", "source3"]):
    term = terminal.rotate(i * 90)
    term.name = name
    terminals.append(term)

device = sc.Device(
    "plus",
    films=[plus],
    layers=[layer],
    terminals={"plus": terminals},
    length_units="um",
)

device.make_mesh(max_edge_length=0.3, smooth=50)

_ = device.plot_mesh(edge_color="k")

def solve_and_plot_model(terminal_currents, applied_field="0 uT"):
    # Define cross-sections for each terminal
    xs = np.linspace(-1.5, 1.5, 201)
    ys = -3 * np.ones_like(xs)
    rs = np.stack([xs, ys], axis=1)
    sections = [rotate(rs, i * 90) for i in range(4)]
    # Put the drain at the end of the list
    sections.append(sections.pop(0))

    # Solve the model
    applied_field = sc.ureg(applied_field).to("uT").magnitude
    solution = sc.solve(
        device,
        terminal_currents=terminal_currents,
        applied_field=sc.sources.ConstantField(applied_field),
        current_units="uA",
        field_units="uT",
    )[-1]

    terminal_currents = solution.terminal_currents["plus"]
    target_currents = list(terminal_currents.values()) + [None]

    # Calculate the total current though each terminal
    for coords, target, name in zip(sections, target_currents, terminal_currents):
        current = solution.current_through_path(coords, film="plus", units="uA").magnitude
        if target is None:
            # This is the drain terminal. It should sink all current.
            target = -sum(terminal_currents.values())
        err = (-current - target) / abs(target) * 100
        print(f"{name}: Total current {-current:.3f} uA ({err:.3f}% error)")

    # Plot currents
    fig, axes = solution.plot_currents(
        streamplot=True,
        min_stream_amp=1e-4,
        cross_section_coords=sections,
        figsize=(5, 7),
    )
    _ = device.plot_polygons(ax=axes[0], lw=1, color="w", ls="--")

    return solution

# Each terminal current can be a float, string, or pint.Quantity
terminal_currents = {
    "source1": "1 uA",
    "source2": sc.ureg("2 uA"),
    "source3": 3,
    "drain": "-6 uA",
}
solution = solve_and_plot_model({"plus": terminal_currents})

source1: Total current 1.007 uA (0.685% error)
source2: Total current 2.003 uA (0.161% error)
source3: Total current 3.004 uA (0.120% error)
drain: Total current -6.014 uA (-0.227% error)

In the configuration below, both drain and source1 will sink \(2\,\mu\mathrm{m}\) of current.

terminal_currents = {
    "source1": "-2 uA",
    "source2": "2 uA",
    "source3": "2 uA",
    "drain": "-2 uA",
}
solution = solve_and_plot_model({"plus": terminal_currents})

source1: Total current -2.003 uA (-0.127% error)
source2: Total current 2.003 uA (0.174% error)
source3: Total current 2.002 uA (0.115% error)
drain: Total current -2.004 uA (-0.187% error)

Of course, this all works with a non-zero applied field too.

terminal_currents = {
    "source1": "1 uA",
    "source2": "2 uA",
    "source3": "3 uA",
    "drain": "-6 uA",
}
solution = solve_and_plot_model({"plus": terminal_currents}, applied_field="5 uT")

source1: Total current 0.997 uA (-0.344% error)
source2: Total current 2.002 uA (0.098% error)
source3: Total current 3.005 uA (0.165% error)
drain: Total current -6.016 uA (-0.273% error)