Quickstart

This notebook is intended to demonstrate the basic usage of superscreen by calculating the magnetic response of several simple single-layer devices.

%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
from superscreen.geometry import circle, box

Device geometry and materials

• Layers represent different physical planes in a device. All superconducting films in a given layer are assumed to have the same thickness and penetration depth.

• Films represent the actual superconducting films, which may fill only part of the Layer that they are in, and may have one or more holes. Each film can have an arbitrary polygonal geometry, and there can be multiple (non-overlapping) films per layer.

• Holes are polygonal regions of vacuum completely surrounded (in 2D) by a superconducting film, which can contain trapped flux.

• Abstract regions are polygonal regions which need not represent any real geometry in the device. Abstract regions will be meshed just like films and holes, and after solving a Device, one can calculate the flux through all abstract regions. Abstract regions can be used to define a “bounding box,” a region of vacuum outside of the convex hull of all the films that will be meshed and simulated.

Superconducting ring with a slit

Here we define a superconducting ring with inner radius 1 \(\mu\)m, outer radius 3 \(\mu\)m, London penetration depth \(\lambda=100\) nm, and thickness \(d=25\) nm, for an effective penetration depth of \(\Lambda=\lambda^2/d=400\) nm. The ring also has a slit of width 250 nm in it.

length_units = "um"
ro = 3  # outer radius
ri = 1  # inner radius
slit_width = 0.25
layer = sc.Layer("base", london_lambda=0.100, thickness=0.025, z0=0)

ring = circle(ro)
hole = circle(ri)
slit = box(slit_width, 1.5 * (ro - ri), center=(0, -(ro + ri) / 2))

film = sc.Polygon.from_difference(
    [ring, slit, hole], name="ring_with_slit", layer="base"
)
# # The above is equivalent to all of the following:
# film = sc.Polygon(
#     name="ring_with_slit", layer="base", points=ring
# ).difference(slit, hole)
# film = sc.Polygon(
#     name="ring_with_slit", layer="base", points=ring
# ).difference(slit).difference(hole)
# film = sc.Polygon(
#     name="ring_with_slit", layer="base", points=ring
# ).difference(sc.Polygon(points=hole).union(slit))

film = film.resample(500)
bounding_box = sc.Polygon("bounding_box", layer="base", points=circle(1.2 * ro))

device = sc.Device(
    film.name,
    layers=[layer],
    films=[film],
    abstract_regions=[bounding_box],
    length_units=length_units,
)

print(device)

Device(
    "ring_with_slit",
    layers=[
        Layer("base", Lambda=0.400, thickness=0.025, london_lambda=0.100, z0=0.000),
    ],
    films=[
        Polygon(name="ring_with_slit", layer="base", points=<ndarray: shape=(500, 2)>, mesh=True),
    ],
    holes=None,
    abstract_regions=[
        Polygon(name="bounding_box", layer="base", points=<ndarray: shape=(101, 2)>, mesh=True),
    ],
    length_units="um",
)

fig, ax = device.draw(exclude="bounding_box", legend=True)

Generate the mesh

device.make_mesh(min_points=4_000, smooth=50)

INFO:superscreen.device.device:Generating mesh...
INFO:superscreen.device.device:Smoothin mesh with 4065 vertices.
INFO:superscreen.device.device:Finished generating mesh with 4065 points and 7971 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.

fig, ax = device.plot(color="k", alpha=0.5, mesh=True, mesh_kwargs=dict(linewidth=0.5))
_ = ax.set_title(
    f"Mesh: {device.points.shape[0]} points, " f"{device.triangles.shape[0]} triangles"
)

Simulate Meissner screening

Here we apply a uniform field field in the \(+z\) direction and calculate the ring’s magnetic response.

applied_field = sc.sources.ConstantField(1)

solutions = sc.solve(
    device=device,
    applied_field=applied_field,
    field_units="uT",
    current_units="uA",
)
assert len(solutions) == 1
solution = solutions[-1]

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

xs = np.linspace(-3.5, 3.5, 401)

cross_section_coords = [
    # [x-coords, y-coords]
    np.stack([xs, 0 * xs], axis=1),  # horizontal cross-section
    np.stack([xs, -2 * np.ones_like(xs)], axis=1),  # horizontal cross-section
    np.stack([0 * xs, xs], axis=1),  # vertical cross-section
]

Visualize the fields

fig, axes = solution.plot_fields(
    cross_section_coords=cross_section_coords, figsize=(4, 5)
)

Visualize the currents

fig, axes = solution.plot_currents(
    cross_section_coords=cross_section_coords,
    figsize=(4, 5),
)

Simulate trapped vortices

We represent a trapped vortex as an instance of the superscreen.Vortex class. Here we assume no applied field.

vortices = [
    sc.Vortex(x=1.5, y=1.5, layer="base"),
    sc.Vortex(x=-1.5, y=-1.5, layer="base"),
    sc.Vortex(x=0, y=2.5, layer="base"),
]

solutions = sc.solve(
    device=device, vortices=vortices, field_units="mT", current_units="mA"
)
assert len(solutions) == 1
solution = solutions[-1]

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

fig, axes = solution.plot_fields()
for ax in axes:
    device.plot(ax=ax, color="w", legend=False)

fig, axes = solution.plot_currents()
for ax in axes:
    device.plot(ax=ax, color="w", legend=False)

Superconducting ring without a slit

Let’s see what happens if we add a hole to our device (making it a ring or washer).

device = sc.Device(
    "ring",
    layers=[sc.Layer("base", london_lambda=0.100, thickness=0.025, z0=0)],
    films=[sc.Polygon("ring", layer="base", points=ring)],
    holes=[sc.Polygon("hole", layer="base", points=hole)],
    abstract_regions=[bounding_box],
    length_units=length_units,
)

device.make_mesh(min_points=4_000, smooth=50)

INFO:superscreen.device.device:Generating mesh...
INFO:superscreen.device.device:Smoothin mesh with 4257 vertices.
INFO:superscreen.device.device:Finished generating mesh with 4257 points and 8341 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.

fig, ax = device.draw(exclude="bounding_box", legend=True)

Trapped flux

We can also solve for the field and current distribution from circulating currents associated with flux trapped in the hole.

We assume there is a total current of 1 mA circulating clockwise in the ring (associated with some positive net trapped flux), and that there is otherwise no applied magnetic field. From here we can calculate the current distribution in the ring, the total magnetic field in the plane of the ring, and the flux through the ring.

Note that, although here we are assuming no applied field, we can also solve models with both trapped flux and applied fields.

circulating_currents = {"hole": "1 mA"}

solutions = sc.solve(
    device,
    circulating_currents=circulating_currents,
    field_units="mT",
    current_units="mA",
)
solution = solutions[-1]

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

fig, axes = solution.plot_fields(
    cross_section_coords=cross_section_coords[:1], figsize=(4, 5)
)

Verify the circulating current by integrating the current density \(\vec{J}\)

fig, axes = solution.plot_currents(
    units="mA/um",
    cross_section_coords=cross_section_coords[:1],
    figsize=(4, 5),
)

# horizontal cut
xs = np.linspace(0, 3.5, 101)
xcut = np.stack([xs, 0 * xs], axis=1)
# vertical cut
ys = np.linspace(-3.5, 0, 101)
ycut = np.stack([0 * ys, ys], axis=1)

for i, (cut, label, axis) in enumerate(zip((xcut, ycut), "xy", (1, 0))):
    # Plot the cut coordinates on the current plot
    for ax in axes:
        ax.plot(cut[:, 0], cut[:, 1], "w--")
    # Evaluate the current density at the cut coordinates
    j_interp = solution.interp_current_density(cut, units="mA / um", with_units=False)[
        "base"
    ]
    # Integrate the approriate component of the current density
    # along the cut to get the total circulating current.
    I_tot = np.trapz(j_interp[:, axis], x=cut[:, i])
    I_target = 1  # mA
    print(
        f"{label} cut, total current: {I_tot:.5f} mA "
        f"(error = {100 * abs((I_tot - I_target) / I_target):.3f}%)"
    )

x cut, total current: 1.00095 mA (error = 0.095%)
y cut, total current: 1.00157 mA (error = 0.157%)

Calculate the ring’s fluxoid and self-inductance

The self-inductance of a superconducting loop with a hole \(h\) is given by

\[L = \frac{\Phi_h^f}{I_h},\]

where \(I_h\) is the current circulating around the hole and \(\Phi_h^f\) is the fluxoid for a path enclosing the hole. The fluxoid \(\Phi^f_S\) for a 2D region \(S\) with 1D boundary \(\partial S\) is given by the sum of flux through \(S\) and the line integral of sheet current around \(\partial S\):

\[\Phi^f_S = \underbrace{\int_S\mu_0H_z(\vec{r})\,\mathrm{d}^2r}_\text{flux part} + \underbrace{\oint_{\partial S}\mu_0\Lambda\vec{J}(\vec{r})\cdot\mathrm{d}\vec{r}}_\text{supercurrent part}.\]

The method Solution.polygon_fluxoid() calculates the fluxoid for an arbitrary closed region in a Device, and the method Solution.hole_fluxoid() calculates the fluxoid for a given hole in a Device (automatically generating an appropriate integration region \(S\) if needed).

For a device with \(N\) holes, there are \(N^2\) mutual inductances between them. This mutual inductance matrix is given by:

\[M_{ij}=\frac{\Phi_i^f}{I_j},\]

where \(\Phi_i^f\) is the fluxoid for a region enclosing hole \(i\) and \(I_j\) is the current circulating around hole \(j\). Note that the mutual inductance matrix is symmetric: \(M_{ij}=M_{ji}\). The diagonals of the mutual inductance matrix are the hole self-inductances: \(M_{ii} = L_i\), the self-inductance of hole \(i\).

The above process of modeling a circulating current for each hole, calculating the fluxoid, and extracting the inductance is automated in the method Device.mutual_inductance_matrix(). In this case, for our device with a single hole, the mutual inductance matrix will be a \(1\times1\) matrix containing only the self-inductance \(L\).

Calculate the fluxoid \(\Phi^f\):

region = circle(2)
fluxoid = solution.polygon_fluxoid(region)["base"]
print(fluxoid)
print(f"The total fluxoid is: {sum(fluxoid):~.3fP}")

Fluxoid(flux_part=<Quantity(1.68418618, 'magnetic_flux_quantum')>, supercurrent_part=<Quantity(1.11921012, 'magnetic_flux_quantum')>)
The total fluxoid is: 2.803 Φ_0

hole_name = list(device.holes)[0]
fluxoid = solution.hole_fluxoid(hole_name)
print(fluxoid)
print(f"The total fluxoid is: {sum(fluxoid):~.3fP}")

Fluxoid(flux_part=<Quantity(1.68418618, 'magnetic_flux_quantum')>, supercurrent_part=<Quantity(1.11921012, 'magnetic_flux_quantum')>)
The total fluxoid is: 2.803 Φ_0

Calculcate self-inductance \(L=\Phi^f / I_\mathrm{circ}\):

I_circ = device.ureg(circulating_currents["hole"])
L = (sum(fluxoid) / I_circ).to("Phi_0 / A")
print(f"The ring's self-inductance is {L:.3f~P} = {L.to('pH'):~.3fP}.")

The ring's self-inductance is 2803.396 Φ_0/A = 5.797 pH.

Calculcate self-inductance using Device.mutual_inductance_matrix():

M = device.mutual_inductance_matrix(units="Phi_0 / A")
print(f"Mutual inductance matrix shape:", M.shape)
L = M[0, 0]
print(f"The ring's self-inductance is {L:.3f~P} = {L.to('pH'):~.3fP}.")

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

Mutual inductance matrix shape: (1, 1)
The ring's self-inductance is 2803.396 Φ_0/A = 5.797 pH.

Solve for a specific fluxoid state: \(\Phi^f=n\Phi_0\)

Current and field distributions for a given fluxoid state \(\Phi^f=n\Phi_0\), where \(\Phi_0\) is the superconducting flux quantum, can be modeled by adjusting the circulating current \(I_\mathrm{circ}\) to realize the desired fluxoid value. This calculation is performed by the function superscreen.find_fluxoid_solution().

Here we solve for the current distribution in the ring for the \(n=0\) fluxoid state (i.e. Meissner state), which can be achieved by cooling the ring through its superconducting transition with no applied field. If a small field is then applied, it is screened by the ring such that the fluxoid remains zero.

# n = 0 fluxoid state, apply a field of 1 mT
solution = sc.find_fluxoid_solution(
    device,
    fluxoids=dict(hole=0),
    applied_field=sc.sources.ConstantField(1),
    field_units="mT",
    current_units="mA",
)
I_circ = solution.circulating_currents["hole"]
fluxoid = sum(solution.hole_fluxoid("hole")).to("Phi_0").magnitude
print(f"Total circulating current: {I_circ:.3f} mA.")
print(f"Total fluxoid: {fluxoid:.6f} Phi_0.")

INFO:superscreen.device.device:Evaluating 'ring' mutual inductance matrix column (1/1), source = 'hole'.
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.

Total circulating current: -1.809 mA.
Total fluxoid: 0.000000 Phi_0.

fig, axes = solution.plot_fields(
    cross_section_coords=cross_section_coords[:1], figsize=(4, 5)
)

fig, axes = solution.plot_currents(
    cross_section_coords=cross_section_coords[:1], figsize=(4, 5)
)

Similarly, we can solve for a non-zero fluxoid state, in this case with no applied field. A non-zero fluxoid state can be realized by cooling a ring through its superconducting transition with an applied field.

# n = 1 fluxoid state, apply a field of 0 mT
solution = sc.find_fluxoid_solution(
    device,
    fluxoids=dict(hole=1),
    applied_field=sc.sources.ConstantField(0),
    field_units="mT",
    current_units="mA",
)
I_circ = solution.circulating_currents["hole"]
fluxoid = sum(solution.hole_fluxoid("hole")).to("Phi_0").magnitude
print(f"Total circulating current: {I_circ:.3f} mA.")
print(f"Total fluxoid: {fluxoid:.6f} Phi_0.")

INFO:superscreen.device.device:Evaluating 'ring' mutual inductance matrix column (1/1), source = 'hole'.
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.

Total circulating current: 0.357 mA.
Total fluxoid: 1.000000 Phi_0.

fig, axes = solution.plot_fields(
    cross_section_coords=cross_section_coords[:1], figsize=(4, 5)
)

Film with multiple holes

Here we simulate a device with fewer symmetries than the ring, namely a rectangular film with two off-center rectangular holes.

length_units = "um"

layers = [
    sc.Layer("base", Lambda=0.1, z0=0),
]

films = [
    sc.Polygon("film", layer="base", points=box(8, 4)),
]

holes = [
    sc.Polygon("hole0", layer="base", points=box(5, 1, center=(0.5, -0.25))).resample(
        101
    ),
    sc.Polygon("hole1", layer="base", points=box(1, 2.5, center=(-3, 0.25))).resample(
        51
    ),
]

abstract_regions = [
    sc.Polygon("bounding_box", layer="base", points=box(9, 5)),
]

device = sc.Device(
    "rect",
    layers=layers,
    films=films,
    holes=holes,
    abstract_regions=abstract_regions,
    length_units=length_units,
)

fig, ax = device.draw(exclude="bounding_box")

device.make_mesh(min_points=4_000)

INFO:superscreen.device.device:Generating mesh...
INFO:superscreen.device.device:Finished generating mesh with 4360 points and 8506 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.

fig, ax = device.plot(mesh=True, mesh_kwargs=dict(linewidth=0.5))
_ = ax.set_title(
    f"Mesh: {device.points.shape[0]} points, " f"{device.triangles.shape[0]} triangles"
)

Full mutual inductance matrix

M = device.mutual_inductance_matrix(units="pH")
print(f"Mutual inductance matrix shape:", M.shape)
display(M)

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

Mutual inductance matrix shape: (2, 2)

Magnitude	[[6.429888172860313 -0.9957703908130666] [-0.9946826119944995 4.875344901946636]]
Units	picohenry

As promised, the mutual inductance matrix is approximately symmetric:

asymmetry = float(np.abs((M[0, 1] - M[1, 0]) / min(M[0, 1], M[1, 0])))
print(f"Mutual inductance matrix fractional asymmetry: {100 * asymmetry:.3f}%")

Mutual inductance matrix fractional asymmetry: 0.109%

Model both holes in the \(n=0\) fluxoid state

# n = 0 fluxoid state, apply a field of 1 mT
solution = sc.find_fluxoid_solution(
    device,
    fluxoids=dict(hole0=0, hole1=0),
    applied_field=sc.sources.ConstantField(1),
    field_units="mT",
    current_units="mA",
)

INFO:superscreen.device.device:Evaluating 'rect' mutual inductance matrix column (1/2), source = 'hole0'.
INFO:superscreen.solve:Calculating base response to applied field.
INFO:superscreen.device.device:Evaluating fluxoids for solution 1/1.
INFO:superscreen.device.device:Evaluating 'rect' mutual inductance matrix column (2/2), source = 'hole1'.
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.

I_circ = solution.circulating_currents
fluxoids = [
    sum(solution.hole_fluxoid(hole)).to("Phi_0").magnitude
    for hole in ("hole0", "hole1")
]
print(f"Total circulating current: {I_circ} mA.")
print(f"Total fluxoid: {fluxoids} Phi_0.")

Total circulating current: {'hole0': -2.1459112092423704, 'hole1': -1.4983020233796058} mA.
Total fluxoid: [-5.662137425588298e-15, 2.55351295663786e-15] Phi_0.

fig, axes = solution.plot_fields(figsize=(6, 2))
fig, axes = solution.plot_currents(figsize=(6, 2))

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:35:06 2023 UTC