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.

• There is a single drain terminal in the system, along which the stream function \(g_\mathrm{transport}(\vec{r})\) changes linearly by a total amount \(\sum_sI_s,\) ensuring current conservation.

• When the terminals are oriented in a counterclockwise direction, the boundary immediately after the drain terminal is the “reference boundary,” along which \(g_\mathrm{transport}(\vec{r})=0\).

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})\) in a film with (effectively) infinite penetration depth \(\Lambda\to\infty\).

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

The effective field \(H_{z,\,\mathrm{eff}}\) is found by embedding the film in a (hypothetical) loop with a given circulating current \(I_\mathrm{circ}\) (see figure below). For \(\Lambda\) much larger than the size of the film and \(I_\mathrm{circ}=I_s,\) this setup produces a uniform current distribution in the film with total current \(I_s\) as desired.

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.10.5
Numpy	1.26.3
Numba	0.58.1
SciPy	1.11.4
matplotlib	3.8.2
IPython	8.18.1
Python	3.9.18 (main, Oct 18 2023, 22:10:26) [GCC 9.3.0]
OS	posix [linux]
Number of CPUs	Physical: 1, Logical: 2
BLAS Info	Generic
Thu Jan 18 18:40:20 2024 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.translate(dx=width - slot_width)).resample(251)

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=50)

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,
)[-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=True,
    cross_section_coords=sections,
    figsize=(3, 6),
)
_ = device.plot_polygons(ax=axes[0], lw=2)

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 = 99.968 uA (-0.03% error)
Cross-section: y = 0.30 um, total current = 100.019 uA (0.02% error)
Cross-section: y = 0.00 um, total current = 99.955 uA (-0.04% 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"}}

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

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(60, 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(60, 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
    applied_field = sc.ureg(applied_field).to("uT").magnitude
    solution = sc.solve(
        device,
        applied_field=sc.sources.ConstantField(applied_field),
        terminal_currents={"film": {"source": source_current, "drain": f"-{source_current}"}},
        field_units="uT",
    )[-1]

    # Define cross sections
    xs = np.linspace(-64 / 2, +64 / 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.659 uA (-0.17% error)
Cross-section: y = 10.00 um, total current = 199.632 uA (-0.18% error)
Cross-section: y = 7.50 um, total current = 199.841 uA (-0.08% error)
Cross-section: y = 0.00 um, total current = 203.186 uA (1.59% 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.733 uA (-0.13% error)
Cross-section: y = 10.00 um, total current = 199.752 uA (-0.12% error)
Cross-section: y = 7.50 um, total current = 199.521 uA (-0.24% error)
Cross-section: y = 0.00 um, total current = 203.289 uA (1.64% 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.005 uA (0.478% error)
source2: Total current 2.010 uA (0.492% error)
source3: Total current 3.015 uA (0.504% error)
drain: Total current -6.030 uA (-0.503% 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.010 uA (-0.507% error)
source2: Total current 2.010 uA (0.485% error)
source3: Total current 2.010 uA (0.503% error)
drain: Total current -2.010 uA (-0.514% 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 1.005 uA (0.501% error)
source2: Total current 1.999 uA (-0.063% error)
source3: Total current 3.006 uA (0.203% error)
drain: Total current -6.034 uA (-0.562% error)