4. Example #2: Simulating an SIS mixer

  • This example will simulate a simple SIS mixer.

  • The input will consists of 3 signals:

    • a strong local-osciallator (LO) signal at 230 GHz,

    • a weak RF signal at 232 GHz, and

    • an IF output at 2 GHz.

  • An arbitrary embedding circuit is included (therefore, harmonic balance is required).

  • The response function is generated from an experimental DC I-V curve.

  • This simulation will calculate:

    • the pumped I-V curve, and

    • the mixer’s gain.

[1]:
%matplotlib inline

import qmix
import numpy as np
import scipy.constants as sc
import matplotlib.pyplot as plt

from qmix.mathfn.misc import slope
from qmix.qtcurrent import qtcurrent
from qmix.circuit import EmbeddingCircuit
from qmix.harmonic_balance import harmonic_balance

# https://github.com/garrettj403/SciencePlots
plt.style.use(['science', 'notebook'])

4.1. Define the junction properties

  • These are common parameters that describe the electrical properties of SIS junctions. QMix only uses normalized values, so these values will be used to de-normalize whatever voltages, currents and frequencies that QMix calculates.

[2]:
vgap = 2.7e-3              # gap voltage in [V]
rn   = 13.5                # normal resistance in [ohms]
igap = vgap / rn           # gap current in [A]
fgap = sc.e * vgap / sc.h  # gap frequency in [Hz]

4.2. Define the embedding circuit

  • We will load the EmbeddingCircuit class with all of the information about the embedding circuit. This includes the Thevenin voltage, the Thevenin impedance, and the frequency of each signal that is applied to the junction.

  • Note:

  • All of the circuit properties are normalized:

    • voltages are normalized to the gap voltage (vgap),

    • resistances to the normal resistance (rn),

    • currents to the gap current (igap = vgap / rn), and

    • frequencies to the gap frequency (fgap).

  • This example includes the embedding circuit. Each tone/harmonic requires:

    • the Thevenin voltage (normalized to vgap)

    • the Thevenin impedance (normalized to rn)

    • the frequency (normalized to fgap).

  • The normalized frequency is the same thing as the normalized photon voltage.

[4]:
# simulation parameters
num_f = 3  # number of tones
num_p = 1  # number of harmonics
num_b = (10, 5, 10)  # Bessel function summation limits

# LO signal
f_lo         = 230e9 / fgap      # frequency (normalized value)
alpha_lo     = 1.2               # junction drive level (normalized value)
impedance_lo = 0.3 - 0.3*1j      # embedding impedance (normalized value)

# RF signal
f_rf         = 232e9 / fgap      # frequency (normalized value)
alpha_rf     = 0.012             # junction drive level (normalized value)
impedance_rf = 0.3 - 0.3*1j      # embedding impedance (normalized value)

# IF signal
f_if         = 2e9 / fgap        # frequency (normalized value)
impedance_if = 1.                # embedding impedance (normalized value)

# build embedding circuit
cct = EmbeddingCircuit(num_f, num_p, fgap=fgap, vgap=vgap, rn=rn)
cct.comment[1][1] = 'LO'
cct.comment[2][1] = 'USB'
cct.comment[3][1] = 'IF'

# photon voltages (same as normalized frequency)
cct.freq[1] = f_lo
cct.freq[2] = f_rf
cct.freq[3] = f_if

# Thevenin voltages
cct.vt[1, 1] = cct.freq[1] * alpha_lo
cct.vt[2, 1] = cct.freq[2] * alpha_rf
cct.vt[3, 1] = 0.

# Thevenin impedances
cct.zt[1, 1] = impedance_lo
cct.zt[2, 1] = impedance_rf
cct.zt[3, 1] = impedance_if

cct.print_info()
Embedding circuit (Tones:3, Harmonics:1)
   f=1, p=1                 230.0 GHz x 1           LO
        Thev. voltage:          0.4228 * Vgap
                                1.2000 * Vph
        Thev. impedance:        0.30-0.30j * Rn
        Avail. power:           4.02E-08 W
                                -43.956 dBm
   f=2, p=1                 232.0 GHz x 1           USB
        Thev. voltage:          0.0043 * Vgap
                                0.0120 * Vph
        Thev. impedance:        0.30-0.30j * Rn
        Avail. power:           4.09E-12 W
                                -83.881 dBm
   f=3, p=1                 2.0 GHz x 1             IF
        Thev. voltage:          0.0000 * Vgap
                                0.0000 * Vph
        Thev. impedance:        1.00+0.00j * Rn
        Avail. power:           0.00E+00 W
                                -inf dBm

4.3. Generate the response function

  • Here, the response function is generated from an experimental DC I-V curve.

[5]:
# Load DC I-V data
dciv = qmix.exp.RawData0('example-data/dciv.csv')

# Grab response function
resp = dciv.resp

DC I-V data:
        Vgap:             2.72  mV
        fgap:           658.72  GHz

        Rn:              13.41  ohms
        Rsg:            368.05  ohms
        Q:               27.44

        Jc:              13.54  kA/cm^2
        Ileak:            8.71  uA

        Offset:           0.10  mV
                          9.67  uA

        Vint:             0.45  mV

[6]:
# Plot
fig, ax = plt.subplots(figsize=(7,5))
resp.plot(ax=ax);
_images/multi-tone-simulation_8_0.png

4.4. Perform harmonic balance

  • Calculate the voltage drop across the junction for each tone/harmonic using harmonic balance.

[7]:
vj = harmonic_balance(cct, resp, num_b)
Running harmonic balance:
 - 3 tone(s) and 1 harmonic(s)
 - 7 calls to the quasiparticle tunneling current (qtc) function per iteration
 - max. iterations: 10
Estimated time:
 - time per qtc call:     0.09 s /   0.00 min /  0.00 hrs
 - time per iteration:    0.61 s /   0.01 min /  0.00 hrs
 - max sim time:          6.07 s /   0.10 min /  0.00 hrs
Error after 0 iteration(s):
        f:1, p:1,   med. rel. error:     0.248,   max. rel. error:     1.013,     0.0 % complete
        f:2, p:1,   med. rel. error:     0.213,   max. rel. error:     0.768,     0.0 % complete
        f:3, p:1,   med. rel. error: 99999.999,   max. rel. error: 99999.999,     0.0 % complete
Calculating inverse Jacobian |--------------------| 100.0% Complete
Applying correction
Error after 1 iteration(s):
        f:1, p:1,   med. rel. error:     0.003,   max. rel. error:     0.134,    29.4 % complete
        f:2, p:1,   med. rel. error:     0.011,   max. rel. error:     0.235,    14.9 % complete
        f:3, p:1,   med. rel. error:     0.117,   max. rel. error:    18.770,     0.0 % complete
Calculating inverse Jacobian |--------------------| 100.0% Complete
Applying correction
Error after 2 iteration(s):
        f:1, p:1,   med. rel. error:     0.000,   max. rel. error:     0.003,    99.5 % complete
        f:2, p:1,   med. rel. error:     0.000,   max. rel. error:     0.005,    98.0 % complete
        f:3, p:1,   med. rel. error:     0.000,   max. rel. error:     0.304,    88.1 % complete
Calculating inverse Jacobian |--------------------| 100.0% Complete
Applying correction
Error after 3 iteration(s):
        f:1, p:1,   med. rel. error:     0.000,   max. rel. error:     0.000,   100.0 % complete
        f:2, p:1,   med. rel. error:     0.000,   max. rel. error:     0.000,   100.0 % complete
        f:3, p:1,   med. rel. error:     0.000,   max. rel. error:     0.000,   100.0 % complete
Done: Minimum error target was achieved.
 - sim time:               2.08 s /   0.03 min /  0.00 hrs
 - 3 iterations required
 - time per iteration:     0.69 s /   0.01 min /  0.00 hrs

4.5. Calculate tunnelling currents

  • idc is the DC tunneling current

  • ilo is the AC tunneling current at the LO frequency

  • iif is the AC tunneling current at the IF frequency

[9]:
freq_list = [0., f_lo, f_if]
idc, ilo, iif = qtcurrent(vj, cct, resp, freq_list, num_b=num_b)
Calculating tunneling current...
 - 3 tone(s)
 - 1 harmonic(s)
Done.
Time: 0.2079 s

4.5.1. Calculate mixer gain

[10]:
# Available power from upper sideband (USB)
pusb = cct.available_power(f=2)

# Power delivered to load
rload = cct.zt[3,1].real * rn
pload = 0.5 * np.abs(iif * igap)**2 * rload

# Gain
gain = pload / pusb

4.6. Plot Results

[11]:
# Conversion: normalized -> de-normalized
vmv = vgap / sc.milli
iua = igap / sc.micro
[12]:
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14,14))

# Plot I-V curves
ax1.plot(resp.voltage*vmv, resp.current*iua, label='Unpumped')
ax1.plot(cct.vb*vmv, idc.real*iua, 'r--', label='Pumped')
ax1.set(xlim=(0,4.5), xlabel='Bias Voltage (mV)')
ax1.set(ylim=(0,300), ylabel='DC Current (uA)')
ax1.legend(frameon=False)

# Plot AC currents at LO frequency
ax2.plot(cct.vb*vmv, np.abs(ilo), 'k', label='Absolute')
ax2.plot(cct.vb*vmv, ilo.real, 'r--', label='Real')
ax2.plot(cct.vb*vmv, ilo.imag, ls='-.', label='Imaginary')
ax2.set(xlim=(0,4.5), xlabel='Bias Voltage (mV)')
ax2.set(ylabel='AC Current (uA)')
ax2.legend(frameon=False, loc=2, title="@ LO frequency")

# Plot AC currents at IF frequency
ax3.plot(cct.vb*vmv, np.abs(iif), 'k', label='Absolute')
ax3.plot(cct.vb*vmv, iif.real, 'r--', label='Real')
ax3.plot(cct.vb*vmv, iif.imag, ls='-.', label='Imaginary')
ax3.set(xlim=(0,4.5), xlabel='Bias Voltage (mV)')
ax3.set(ylabel='AC Current (uA)')
ax3.legend(frameon=False, loc=1, title="@ IF frequency")

# Plot gain of SIS junction
ax4.plot(cct.vb*vmv, gain*100)
ax4.set(xlim=(0,4.5), xlabel='Bias Voltage (mV)')
ax4.set(ylabel='Gain (%)')
ax4.set_ylim(bottom=0)

fig.savefig('multi-tone-results.png', dpi=500, bbox_inches='tight');
_images/multi-tone-simulation_17_0.png