Quickstart
This notebook is intended to demonstrate the basic usage of superscreen by calculating the magnetic response of several simple single-layer and multi-layer devices.
[1]:
# 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
[2]:
%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 circle, box
[3]:
sc.version_table()
[3]:
| Software | Version |
|---|---|
| SuperScreen | 0.10.2 |
| Numpy | 1.25.2 |
| Numba | 0.58.0 |
| SciPy | 1.11.3 |
| matplotlib | 3.8.0 |
| IPython | 8.16.0 |
| Python | 3.9.17 (main, Jul 27 2023, 09:05:49) [GCC 9.3.0] |
| OS | posix [linux] |
| Number of CPUs | Physical: 1, Logical: 2 |
| BLAS Info | OPENBLAS |
| Fri Sep 29 17:16:04 2023 UTC | |
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, each of which exists in a
Layerand 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 within a film which need not represent any real geometry in the device. Abstract regions can be used to locally increase the finite element mesh density inside a film or hole.
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.
[4]:
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))
device = sc.Device(
film.name,
layers=[layer],
films=[film],
length_units=length_units,
)
[5]:
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=(263, 2)>),
],
holes=None,
terminals=None,
abstract_regions=None,
length_units="um",
)
[6]:
fig, ax = device.draw(legend=True)
Generate the mesh
[7]:
device.make_mesh(max_edge_length=0.15, smooth=100)
[8]:
device.mesh_stats()
[8]:
| length_units | 'um' | |
| 'ring_with_slit' | num_sites | 5639 |
| num_elements | 11014 | |
| min_edge_length | 1.130e-03 | |
| max_edge_length | 1.425e-01 | |
| min_vertex_area | 1.374e-06 | |
| max_vertex_area | 1.441e-02 |
[9]:
fig, ax = device.plot_mesh(show_sites=False)
_ = device.plot_polygons(ax=ax, color="k")
Simulate Meissner screening
Here we apply a uniform field field in the \(+z\) direction and calculate the ring’s magnetic response.
[10]:
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]
[11]:
xs = np.linspace(-3.5, 3.5, 401)
cross_section_coords = [
# [x-coords, y-coords]
np.array([xs, 0 * xs]).T, # horizontal cross-section
np.array([xs, -2 * np.ones_like(xs)]).T, # horizontal cross-section
np.array([0 * xs, xs]).T, # vertical cross-section
]
Visualize the fields
[12]:
fig, axes = solution.plot_fields(
cross_section_coords=cross_section_coords, figsize=(4, 5)
)
for ax in axes:
_ = device.plot_polygons(ax=ax, color="w", ls="--", lw=1)
Visualize the currents
[13]:
fig, axes = solution.plot_currents(
cross_section_coords=cross_section_coords,
figsize=(4, 5),
)
for ax in axes:
_ = device.plot_polygons(ax=ax, color="w", ls="--", lw=1)
Simulate trapped vortices
We represent a trapped vortex as an instance of the superscreen.Vortex class. Here we assume no applied field.
[14]:
vortices = [
sc.Vortex(x=1.5, y=1.5, film="ring_with_slit"),
sc.Vortex(x=-1.5, y=-1.5, film="ring_with_slit"),
sc.Vortex(x=0, y=2.5, film="ring_with_slit"),
]
solutions = sc.solve(
device=device, vortices=vortices, field_units="mT", current_units="mA"
)
assert len(solutions) == 1
solution = solutions[-1]
[15]:
fig, axes = solution.plot_fields(shading="gouraud")
for ax in axes:
_ = device.plot_polygons(ax=ax, color="w", ls="--", lw=1)
[16]:
fig, axes = solution.plot_currents(shading="gouraud")
for ax in axes:
_ = device.plot_polygons(ax=ax, color="w", ls="--", lw=1)
Superconducting ring without a slit
Let’s see what happens if we add a hole to our device (making it a ring or washer).
[17]:
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)],
length_units=length_units,
)
[18]:
device.make_mesh(max_edge_length=0.15, smooth=100)
[19]:
fig, ax = device.draw(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.
[20]:
circulating_currents = {"hole": "1 mA"}
solutions = sc.solve(
device,
circulating_currents=circulating_currents,
field_units="mT",
current_units="mA",
)
solution = solutions[-1]
[21]:
fig, axes = solution.plot_fields(
cross_section_coords=cross_section_coords[:1], figsize=(4, 5)
)
for ax in axes:
_ = device.plot_polygons(ax=ax, color="w", ls="--", lw=1)
Verify the circulating current by integrating the current density \(\vec{J}\)
[22]:
fig, axes = solution.plot_currents(
units="mA/um",
cross_section_coords=cross_section_coords[:1],
figsize=(4, 5),
)
for ax in axes:
_ = device.plot_polygons(ax=ax, color="w", ls="--", lw=1)
# horizontal cut
xs = np.linspace(0, 3.15, 151)
xcut = np.array([xs, 0 * xs]).T
# vertical cut
ys = np.linspace(-3.15, 0, 101)
ycut = np.array([0 * ys, ys]).T
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, film="ring", units="mA / um", with_units=False)
# 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: 0.99683 mA (error = 0.317%)
y cut, total current: 0.99223 mA (error = 0.777%)
Calculate the ring’s fluxoid and self-inductance
The self-inductance of a superconducting loop with a hole \(h\) is given by
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\):
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:
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\):
[23]:
region = circle(2)
fluxoid = solution.polygon_fluxoid(region, film="ring")
print(fluxoid)
print(f"The total fluxoid is: {sum(fluxoid):~.3fP}")
Fluxoid(flux_part=<Quantity(1.67636329, 'magnetic_flux_quantum')>, supercurrent_part=<Quantity(1.12062425, 'magnetic_flux_quantum')>)
The total fluxoid is: 2.797 Φ_0
[24]:
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.67636329, 'magnetic_flux_quantum')>, supercurrent_part=<Quantity(1.12062425, 'magnetic_flux_quantum')>)
The total fluxoid is: 2.797 Φ_0
Calculcate self-inductance \(L=\Phi^f / I_\mathrm{circ}\):
[25]:
I_circ = solution.circulating_currents["hole"] * device.ureg(solution.current_units)
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 2796.988 Φ_0/A = 5.784 pH.
Calculcate self-inductance using Device.mutual_inductance_matrix():
[26]:
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}.")
Mutual inductance matrix shape: (1, 1)
The ring's self-inductance is 2796.988 Φ_0/A = 5.784 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.
[27]:
# 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.")
Total circulating current: -1.817 mA.
Total fluxoid: -0.000001 Phi_0.
[28]:
fig, axes = solution.plot_fields(
cross_section_coords=cross_section_coords[:1], figsize=(4, 5)
)
for ax in axes:
_ = device.plot_polygons(ax=ax, color="w", ls="--", lw=1)
[29]:
fig, axes = solution.plot_currents(
cross_section_coords=cross_section_coords[:1], figsize=(4, 5)
)
for ax in axes:
_ = device.plot_polygons(ax=ax, color="w", ls="--", lw=1)
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.
[30]:
# 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.")
Total circulating current: 0.358 mA.
Total fluxoid: 1.000000 Phi_0.
[31]:
fig, axes = solution.plot_fields(
cross_section_coords=cross_section_coords[:1], figsize=(4, 5)
)
for ax in axes:
_ = device.plot_polygons(ax=ax, color="w", ls="--", lw=1)
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.
[32]:
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))).buffer(0),
sc.Polygon("hole1", layer="base", points=box(1, 2.5, center=(-3, 0.25))).buffer(0),
]
device = sc.Device(
"rect",
layers=layers,
films=films,
holes=holes,
length_units=length_units,
)
[33]:
fig, ax = device.draw()
[34]:
device.make_mesh(max_edge_length=0.15)
[35]:
device.mesh_stats()
[35]:
| length_units | 'um' | |
| 'film' | num_sites | 6330 |
| num_elements | 12387 | |
| min_edge_length | 5.000e-02 | |
| max_edge_length | 1.494e-01 | |
| min_vertex_area | 1.499e-03 | |
| max_vertex_area | 1.261e-02 |
Full mutual inductance matrix
[36]:
M = device.mutual_inductance_matrix(units="pH")
print(f"Mutual inductance matrix shape:", M.shape)
display(M)
Holes: 100%|██████████| 2/2 [00:01<00:00, 1.77it/s]
Mutual inductance matrix shape: (2, 2)
| Magnitude | [[6.553566078634796 -0.9880575716896793] |
|---|---|
| Units | picohenry |
As promised, the mutual inductance matrix is approximately symmetric:
[37]:
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.413%
Model both holes in the \(n=0\) fluxoid state
[38]:
# 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",
)
Holes: 100%|██████████| 2/2 [00:01<00:00, 1.86it/s]
[39]:
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.103772447656946, 'hole1': -1.4472227512991362} mA.
Total fluxoid: [1.2230384086819868e-05, 2.5776829958412684e-05] Phi_0.
[40]:
fig, axes = solution.plot_fields(figsize=(6, 2))
fig, axes = solution.plot_currents(figsize=(6, 2))
for ax in axes:
_ = device.plot_polygons(ax=ax, color="w", ls="--", lw=1)
Devices with multiple films
In this example we model an on-chip flux bias line (in layer W0), which is partially covered by a shield in another layer (W1). We use SuperScreen to calculate the mutual inductance between the bias line and a pickup loop which also lies in W0. Both layers are assumed to have a London penetration depth of 100 nm and thickness of 50 nm, and the vertical spacing between the layers is 250 nm.
[41]:
# Wiring layers
W0 = sc.Layer("W0", london_lambda=0.1, thickness=0.05, z0=0.0)
W1 = sc.Layer("W1", london_lambda=0.1, thickness=0.05, z0=0.25)
# Pickup loop geometry
pickup_loop = sc.Polygon("pickup_loop", layer="W0", points=box(5, 10))
pickup_loop_center = sc.Polygon("pickup_loop_center", layer="W0", points=box(3, 8))
# Bias line and shield
bias_line = sc.Polygon(
"bias_line", layer="W0", points=box(10, 5)
).difference(box(10, 4, center=(0.5, 0))).resample(501).translate(dx=8)
bias_line_shield = sc.Polygon(
"bias_line_shield", layer="W1", points=box(9, 6)
).translate(dx=10)
# Transport terminals
source = sc.Polygon("source", points=box(0.1, 1)).translate(dx=13, dy=2.25)
drain = source.scale(yfact=-1).set_name("drain")
device = sc.Device(
"bias_loop",
layers=[W0, W1],
films=[bias_line, bias_line_shield, pickup_loop],
holes=[pickup_loop_center],
terminals={"bias_line": [source, drain]},
length_units="um",
)
[42]:
fig, ax = device.draw(figsize=(6, 4))
_ = device.plot_polygons(ax=ax, legend=True)
[43]:
device.make_mesh(max_edge_length=0.25)
[44]:
device.mesh_stats()
[44]:
| length_units | 'um' | |
| 'bias_line' | num_sites | 1169 |
| num_elements | 1836 | |
| min_edge_length | 7.087e-02 | |
| max_edge_length | 2.195e-01 | |
| min_vertex_area | 1.298e-03 | |
| max_vertex_area | 2.475e-02 | |
| 'bias_line_shield' | num_sites | 3544 |
| num_elements | 6888 | |
| min_edge_length | 9.834e-02 | |
| max_edge_length | 2.490e-01 | |
| min_vertex_area | 6.199e-03 | |
| max_vertex_area | 3.847e-02 | |
| 'pickup_loop' | num_sites | 3484 |
| num_elements | 6764 | |
| min_edge_length | 8.666e-02 | |
| max_edge_length | 2.484e-01 | |
| min_vertex_area | 6.703e-03 | |
| max_vertex_area | 3.663e-02 |
[45]:
fig, ax = device.plot_mesh()
_ = device.plot_polygons(ax=ax, color="k")
Solve the model with 1 mA flowing through the bias line.
[46]:
I_bias = "1 mA"
solutions = sc.solve(
device,
terminal_currents={"bias_line": {"source": I_bias, "drain": f"-{I_bias}"}},
iterations=10,
)
Solver iterations: 100%|██████████| 10/10 [00:09<00:00, 1.09it/s]
Plot the sheet current density in each film.
[47]:
for film in device.films:
fig, axes = solutions[-1].plot_currents(films=[film])
_ = device.plot_polygons(ax=axes[0], color="w", ls="--", lw=1)
Evaluate the field at a plane 500 nm above the W0 layer and 250 nm above the W1 layer. The effect of the shield is clearly visible.
[48]:
# Generate a mesh of points at which to evaluate the field from the films.
eval_region = sc.Polygon(points=box(20, 15, center=(5, 0)))
eval_mesh = eval_region.make_mesh(min_points=2000)
fig, ax = solutions[-1].plot_field_at_positions(
eval_mesh,
zs=0.5,
figsize=(6, 3),
symmetric_color_scale=True,
cmap="coolwarm",
)
for film in device.films.values():
film.plot(ax=ax, color="k", ls="--", lw=1)
for hole in device.holes.values():
hole.plot(ax=ax, color="k", ls="--", lw=1)
Evaluate and plot the mutual inductance between the bias line and the loop as a function of the solver iteration. The solution converges very quickly to a value around -0.7 pH.
[49]:
mutual_inductances = []
for solution in solutions:
fluxoid = sum(solution.hole_fluxoid("pickup_loop_center"))
mutual_inductance = fluxoid / sc.ureg(I_bias)
mutual_inductances.append(mutual_inductance.to("pH").magnitude)
print(f"Mutual inductance = {mutual_inductances[-1]:.4f} pH")
fig, (ax, bx) = plt.subplots(1, 2, sharex=True, figsize=(6, 3), constrained_layout=True)
ax.plot(range(len(mutual_inductances)), mutual_inductances, "o--")
ax.grid(True)
ax.set_title("Mutual inductance value")
ax.set_xlabel("Solver iteration")
ax.set_ylabel("Mutual inductance [pH]")
diff = np.diff(mutual_inductances)
bx.plot(range(1, len(mutual_inductances)), np.abs(diff / mutual_inductances[1:]), ".-")
bx.set_yscale("log")
bx.grid(True)
bx.set_title("Mutual inductance convergence")
bx.set_xlabel("Solver iteration")
_ = bx.set_ylabel("Fractional change\nin mutual inductance")
Mutual inductance = -0.7177 pH
[ ]: