# Quickstart

This notebook is intended to demonstrate the basic usage of `superscreen`

by calculating the magnetic response of several simple single-layer devices.

```
[1]:
```

```
%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"] = (8, 6)
plt.rcParams["font.size"] = 14
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.

```
[2]:
```

```
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,
)
```

```
[3]:
```

```
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",
)
```

```
[4]:
```

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

### Generate the mesh

```
[5]:
```

```
device.make_mesh(min_points=4_000, optimesh_steps=10)
```

```
INFO:superscreen.device.device:Generating mesh...
INFO:superscreen.device.device:Optimizing mesh with 4364 vertices.
INFO:superscreen.device.device:Finished generating mesh with 4364 points and 8559 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.
```

```
[6]:
```

```
fig, ax = device.plot(color="k", alpha=0.5, mesh=True)
_ = 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.

```
[7]:
```

```
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.
```

```
[8]:
```

```
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

```
[9]:
```

```
fig, axes = solution.plot_fields(
cross_section_coords=cross_section_coords, figsize=(6, 8)
)
```

#### Visualize the currents

```
[10]:
```

```
fig, axes = solution.plot_currents(
cross_section_coords=cross_section_coords,
figsize=(6, 8),
)
```

### Simulate trapped vortices

We represent a trapped vortex as an instance of the `superscreen.Vortex`

class. Here we assume no applied field.

```
[11]:
```

```
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.
```

```
[12]:
```

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

```
[13]:
```

```
fig, axes = solution.plot_currents(figsize=(6, 5))
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).

```
[14]:
```

```
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,
)
```

```
[15]:
```

```
device.make_mesh(min_points=4_000, optimesh_steps=10)
```

```
INFO:superscreen.device.device:Generating mesh...
INFO:superscreen.device.device:Optimizing mesh with 4153 vertices.
INFO:superscreen.device.device:Finished generating mesh with 4153 points and 8125 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.
```

```
[16]:
```

```
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.

```
[17]:
```

```
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.
```

```
[18]:
```

```
fig, axes = solution.plot_fields(
cross_section_coords=cross_section_coords[:1], figsize=(6, 8)
)
```

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

```
[19]:
```

```
fig, axes = solution.plot_currents(
units="mA/um",
cross_section_coords=cross_section_coords[:1],
figsize=(6, 8),
)
# 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: 0.99889 mA (error = 0.111%)
y cut, total current: 0.99979 mA (error = 0.021%)
```

### 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\):

```
[20]:
```

```
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.69264496, 'magnetic_flux_quantum')>, supercurrent_part=<Quantity(1.12679239, 'magnetic_flux_quantum')>)
The total fluxoid is: 2.819 Φ_0
```

```
[21]:
```

```
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.69264496, 'magnetic_flux_quantum')>, supercurrent_part=<Quantity(1.12679239, 'magnetic_flux_quantum')>)
The total fluxoid is: 2.819 Φ_0
```

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

```
[22]:
```

```
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 2819.437 Φ_0/A = 5.830 pH.
```

Calculcate self-inductance using `Device.mutual_inductance_matrix()`

:

```
[23]:
```

```
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 2819.437 Φ_0/A = 5.830 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.

```
[24]:
```

```
# 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.805 mA.
Total fluxoid: -0.000000 Phi_0.
```

```
[25]:
```

```
fig, axes = solution.plot_fields(
cross_section_coords=cross_section_coords[:1], figsize=(6, 8)
)
```

```
[26]:
```

```
fig, axes = solution.plot_currents(
cross_section_coords=cross_section_coords[:1], figsize=(6, 8)
)
```

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.

```
[27]:
```

```
# 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.355 mA.
Total fluxoid: 1.000000 Phi_0.
```

```
[28]:
```

```
fig, axes = solution.plot_fields(
cross_section_coords=cross_section_coords[:1], figsize=(6, 8)
)
```

## 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.

```
[29]:
```

```
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,
)
```

```
[30]:
```

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

```
[31]:
```

```
device.make_mesh(min_points=4_000, optimesh_steps=None)
```

```
INFO:superscreen.device.device:Generating mesh...
INFO:superscreen.device.device:Finished generating mesh with 4084 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.
```

```
[32]:
```

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

### Full mutual inductance matrix

```
[33]:
```

```
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.411673813197216 -0.9926728089995726] |
---|---|

Units | picohenry |

As promised, the mutual inductance matrix is approximately symmetric:

```
[34]:
```

```
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.089%
```

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

```
[35]:
```

```
# 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.
```

```
[36]:
```

```
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.1518894168825073, 'hole1': -1.5054947232162428} mA.
Total fluxoid: [7.216449660063518e-15, -1.1102230246251565e-16] Phi_0.
```

```
[37]:
```

```
fig, axes = solution.plot_fields(figsize=(8, 3))
fig, axes = solution.plot_currents(figsize=(8, 3))
```

```
[38]:
```

```
sc.version_table()
```

```
[38]:
```

Software | Version |
---|---|

SuperScreen | 0.7.0 |

Numpy | 1.23.2 |

SciPy | 1.9.1 |

matplotlib | 3.5.3 |

ray | 2.0.0 |

jax | 0.3.17 |

IPython | 8.5.0 |

Python | 3.8.6 (default, Oct 19 2020, 15:10:29) [GCC 7.5.0] |

OS | posix [linux] |

Number of CPUs | Physical: 1, Logical: 2 |

BLAS Info | OPENBLAS |

Thu Sep 08 22:07:27 2022 UTC |

```
[ ]:
```

```
```