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 the device through one or more “source terminals” and out of the device through a “drain terminal.”

In this notebook we demonstrate the treatment of terminal currents in a superscreen model. Terminal currents can only be modeled for instances of superscreen.TransportDevice, which is a subclass of superscreen.Device that contains only one film and no bounding box.

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

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

import os

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

import logging

logging.basicConfig(level=logging.INFO)

import numpy as np
import matplotlib.pyplot as plt

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

import superscreen as sc

Film with a constriction

Define the device geometry.

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

film = (
    sc.Polygon("film", layer="base", points=sc.geometry.box(width, height))
    .difference(
        sc.geometry.box(slot_width, slot_height, center=(-(width - slot_width) / 2, 0))
    )
    .difference(
        sc.geometry.box(slot_width, slot_height, center=(+(width - slot_width) / 2, 0))
    )
    .resample(251)
)

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

device = sc.TransportDevice(
    "constriction",
    layer=sc.Layer("base", Lambda=0.1),
    film=film,
    source_terminals=[source_terminal],
    drain_terminal=drain_terminal,
    length_units=length_units,
).translate(dx=dx, dy=dy)

device.make_mesh(min_points=2000, smooth=50)

INFO:superscreen.device.device:Generating mesh...
INFO:superscreen.device.device:Smoothin mesh with 2161 vertices.
INFO:superscreen.device.device:Finished generating mesh with 2161 points and 4069 triangles.
INFO:superscreen.device.device:Calculating weight matrix.
INFO:superscreen.device.device:Calculating Laplace operator.
INFO:superscreen.device.device:Calculating kernel matrix.
INFO:superscreen.device.device:Calculating gradient matrix.

_ = device.plot(mesh=True, mesh_kwargs=dict(color="k", lw=0.5))

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

terminal_currents = {"source": "100 uA"}

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

INFO:superscreen.solve:Calculating base response to applied field.

Plot the stream function.

fig, axes = solution.plot_streams()
for polygon in device.polygons.values():
    polygon.plot(ax=axes[0], lw=2)

Plot the supercurrent distribution.

xs = np.linspace(-width / 2, +width / 2, 201)
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),
)
for polygon in device.polygons.values():
    polygon.plot(ax=axes[0], lw=2)

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

for coords in sections:
    J = solution.interp_current_density(
        coords,
        units="uA/um",
        with_units=False,
    )["base"]
    _, unit_normals = sc.geometry.path_vectors(coords)
    dr = np.linalg.norm(np.diff(coords, axis=0), axis=1)[0]
    total_current = np.sum(J * dr * unit_normals)
    target_current = solution.terminal_currents["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.455 uA (0.45 % error)
Cross-section: y = 0.30 um, total current = 100.208 uA (0.21 % error)
Cross-section: y = 0.00 um, total current = 99.998 uA (-0.00 % error)

Film with holes

device.holes = {
    "hole1": sc.Polygon(
        "hole1",
        layer="base",
        points=sc.geometry.circle(width / 4, center=(-width / 8, +height / 4)),
    ),
    "hole2": sc.Polygon(
        "hole2",
        layer="base",
        points=sc.geometry.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 = {"source": "100 uA"}

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

INFO:superscreen.device.device:Evaluating 'constriction' mutual inductance matrix column (1/2), source = 'hole1'.
INFO:superscreen.solve:Calculating base response to applied field.
INFO:superscreen.device.device:Evaluating fluxoids for solution 1/1.
INFO:superscreen.device.device:Evaluating 'constriction' mutual inductance matrix column (2/2), source = 'hole2'.
INFO:superscreen.solve:Calculating base response to applied field.
INFO:superscreen.device.device:Evaluating fluxoids for solution 1/1.
INFO:superscreen.solve:Calculating base response to applied field.
INFO:superscreen.solve:Calculating base response to applied field.

fig, axes = solution.plot_streams()
for polygon in device.polygons.values():
    polygon.plot(ax=axes[0], lw=2)

fig, axes = solution.plot_currents(streamplot=True)
for polygon in device.polygons.values():
    polygon.plot(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=sc.geometry.box(60, 24)).union(
    sc.geometry.box(10, 40)
)
w_h, h = 4, 2
base_hole = sc.Polygon(layer="base", points=sc.geometry.box(w_h, 2 * w_h))
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=sc.geometry.box(10, 0.1, center=(0, 20)))
drain = sc.Polygon("drain", points=sc.geometry.box(10, 0.1, center=(0, -20)))

device = sc.TransportDevice(
    "squid_array",
    film=film,
    layer=layer,
    holes=holes,
    source_terminals=[source],
    drain_terminal=drain,
    length_units="um",
)

fig, ax = device.draw()

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

INFO:superscreen.device.device:Generating mesh...
INFO:superscreen.device.device:Smoothin mesh with 3092 vertices.
INFO:superscreen.device.device:Finished generating mesh with 3092 points and 5999 triangles.
INFO:superscreen.device.device:Calculating weight matrix.
INFO:superscreen.device.device:Calculating Laplace operator.
INFO:superscreen.device.device:Calculating kernel matrix.
INFO:superscreen.device.device:Calculating gradient matrix.

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=dict(source=source_current),
        field_units="uT",
    )[-1]

    # Define cross sections
    xs = np.linspace(-64 / 2, +64 / 2, 501)
    ys = np.ones_like(xs)
    sections = [
        np.stack([xs, 15.0 * ys], axis=1),
        np.stack([xs, 10.0 * ys], axis=1),
        np.stack([xs, 7.5 * ys], axis=1),
        np.stack([xs, 0.0 * ys], axis=1),
    ]

    # Plot currents.
    fig, axes = solution.plot_currents(
        streamplot=True,
        min_stream_amp=1e-4,
        cross_section_coords=sections,
        figsize=(5, 4),
        grid_method="linear",
    )
    for polygon in device.polygons.values():
        polygon.plot(ax=axes[0], lw=1, color="w", ls="--")

    # Evaluate the current through each cross section.
    for coords in sections:
        J = solution.interp_current_density(
            coords,
            units="uA/um",
            with_units=False,
        )["base"]
        _, unit_normals = sc.geometry.path_vectors(coords)
        dr = np.linalg.norm(np.diff(coords, axis=0), axis=1)[0]
        total_current = np.sum(J * dr * unit_normals)
        target_current = solution.terminal_currents["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
    dx = dy = 0.5
    xs = np.arange(-30, 30 + dx, dx)
    ys = np.arange(-20, 20 + dy, dy)
    X, Y = np.meshgrid(xs, ys)
    points = np.stack([X.ravel(), Y.ravel()], axis=1)

    fig, axes = solution.plot_field_at_positions(
        points,
        zs=0.75,
        cross_section_coords=sections,
        figsize=(5, 4),
    )
    for polygon in device.polygons.values():
        polygon.plot(ax=axes[0], 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")

INFO:superscreen.solve:Calculating base response to applied field.

Cross-section: y = 15.00 um, total current = 199.964 uA (-0.02 % error)
Cross-section: y = 10.00 um, total current = 199.950 uA (-0.02 % error)
Cross-section: y = 7.50 um, total current = 200.080 uA (0.04 % error)
Cross-section: y = 0.00 um, total current = 199.538 uA (-0.23 % error)

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

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

INFO:superscreen.solve:Calculating base response to applied field.

Cross-section: y = 15.00 um, total current = 200.069 uA (0.03 % error)
Cross-section: y = 10.00 um, total current = 199.752 uA (-0.12 % error)
Cross-section: y = 7.50 um, total current = 199.946 uA (-0.03 % error)
Cross-section: y = 0.00 um, total current = 199.098 uA (-0.45 % 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 = sc.geometry.box(width, height)
bar = sc.Polygon("plus", points=points)
plus = bar.union(bar.rotate(90))
plus.name = "plus"
plus.layer = layer.name
terminal = sc.Polygon(
    points=sc.geometry.box(height, width / 100, 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)
drain, *sources = terminals
device = sc.TransportDevice(
    "plus",
    film=plus,
    layer=layer,
    source_terminals=sources,
    drain_terminal=drain,
    length_units="um",
)

device.make_mesh(min_points=2000, smooth=50)

INFO:superscreen.device.device:Generating mesh...
INFO:superscreen.device.device:Smoothin mesh with 2053 vertices.
INFO:superscreen.device.device:Finished generating mesh with 2053 points and 3848 triangles.
INFO:superscreen.device.device:Calculating weight matrix.
INFO:superscreen.device.device:Calculating Laplace operator.
INFO:superscreen.device.device:Calculating kernel matrix.
INFO:superscreen.device.device:Calculating gradient matrix.

_ = device.plot(mesh=True, mesh_kwargs=dict(color="k", lw=0.5))

def solve_and_plot_model(terminal_currents, applied_field="0 uT"):
    # Define cross-sections for each terminal
    xs = np.linspace(-2, 2, 201)
    ys = -3 * np.ones_like(xs)
    rs = np.stack([xs, ys], axis=1)
    sections = [sc.geometry.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]
    target_currents = list(solution.terminal_currents.values()) + [None]

    # Calculate the total current though each terminal
    for coords, target, name in zip(sections, target_currents, device.terminals):
        J = solution.interp_current_density(
            coords,
            units="uA/um",
            with_units=False,
        )["base"]
        _, unit_normals = sc.geometry.path_vectors(coords)
        dr = np.linalg.norm(np.diff(coords, axis=0), axis=1)[0]
        current = np.sum(J * dr * unit_normals)
        if target is None:
            # This is the drain terminal. It should sink all current.
            target = -sum(solution.terminal_currents.values())
        err = (-current - target) / abs(target) * 100
        print(f"{name}: Total current {-current:.3f} uA ({err:.2e} % error)")

    # Plot currents
    fig, axes = solution.plot_currents(
        streamplot=True,
        min_stream_amp=1e-4,
        cross_section_coords=sections,
        figsize=(5, 7),
        grid_method="linear",
    )
    for polygon in device.polygons.values():
        polygon.plot(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,
}
solution = solve_and_plot_model(terminal_currents)

INFO:superscreen.solve:Calculating base response to applied field.

source1: Total current 1.000 uA (2.38e-04 % error)
source2: Total current 2.000 uA (1.11e-03 % error)
source3: Total current 3.000 uA (2.14e-03 % error)
drain: Total current -6.000 uA (-7.19e-04 % 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",
}
solution = solve_and_plot_model(terminal_currents)

INFO:superscreen.solve:Calculating base response to applied field.

source1: Total current -2.000 uA (-2.01e-03 % error)
source2: Total current 2.000 uA (9.75e-04 % error)
source3: Total current 2.000 uA (1.86e-03 % error)
drain: Total current -2.000 uA (6.91e-04 % error)

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

terminal_currents = {
    "source1": "1 uA",
    "source2": "2 uA",
    "source3": "3 uA",
}
solution = solve_and_plot_model(terminal_currents, applied_field="4 uT")

INFO:superscreen.solve:Calculating base response to applied field.

source1: Total current 1.000 uA (1.92e-02 % error)
source2: Total current 1.999 uA (-6.13e-02 % error)
source3: Total current 2.999 uA (-2.14e-02 % error)
drain: Total current -6.001 uA (-9.63e-03 % error)

sc.version_table()

Software	Version
SuperScreen	0.8.0
Numpy	1.24.1
SciPy	1.10.0
matplotlib	3.6.2
ray	2.2.0
jax	0.4.1
IPython	8.8.0
Python	3.9.15 (main, Oct 26 2022, 11:17:18) [GCC 9.3.0]
OS	posix [linux]
Number of CPUs	Physical: 1, Logical: 2
BLAS Info	OPENBLAS
Wed Jan 11 19:36:19 2023 UTC