Support for loops threaded by external flux? (SQUIDs, SNAILs, etc.)

[ttshaw1]

Sorry if this already exists in the program, but I don't see any evidence in the examples that we can specify external flux through a loop of JJs. I can imagine that this would be hard to integrate into the SPICE netlist format. If this is any more feasible, what we do in AWR Microwave Office is to represent a loop with a nonlinear inductor element, with inductance L(I)=L0+IL1+I^2L2+... Then based on the elements in the loop and the external flux we can calculate the parameters that will make the nonlinear inductor equivalent to the loop.

[kpobrien]

Hi @ttshaw1, thanks for posting this discussion topic. The examples show functionality which has been validated through careful comparison with other software, ideally in a publication. More advanced features like flux tuning, 3WM, and multi-tone harmonic balance exist, but have not undergone that level of testing.

Given that disclaimer, below is an unvalidated example of a flux pumped JPA in which flux is coupled into the SQUID loop through mutual inductors with parameters based on Kaiyong He et al 2023 Supercond. Sci. Technol. 36 045010.

The next steps in my opinion would be:
(1) Crosscheck the upper simulation of gain vs frequency using WRspice or another simulation software.
(2) Understand why the solver isn't returning the correct solution for DC currents resulting in a more than half a flux quantum, in the lower simulation.

using JosephsonCircuits
using Plots

# Circuit and parameters from Kaiyong He et al 2023 Supercond. Sci. Technol. 36 045010

# define an inductively JPA with a squid
circuit = Array{Tuple{String,String,String,Union{Complex{Float64},Symbol,Int64}},1}(undef,0)

# first port
push!(circuit,("P1_0","1","0",1))
push!(circuit,("R1_0","1","0",:Rs))

# capacitor to ground
push!(circuit,("C1_0","1","0",:Cp))

# first JJ and inductor
push!(circuit,("L1_2","1","2",:Lg1))
push!(circuit,("R2_0","2","0",:RL))
push!(circuit,("Lj2_0","2","0",:Lj1)) 
push!(circuit,("C2_0","2","0",:Cj))

# second JJ and inductor
push!(circuit,("L1_3","1","3",:Lg2))
push!(circuit,("R3_0","3","0",:RL))
push!(circuit,("Lj3_0","3","0",:Lj2)) 
push!(circuit,("C3_0","3","0",:Cj))

# second port
push!(circuit,("P2","4","0",2))
push!(circuit,("R3","4","0",:Rp))
push!(circuit,("L4_0","4","0",:Lg2))
push!(circuit,("K1","L1_3","L4_0",:K))
                
circuitdefs = Dict{Symbol,Complex{Float64}}(
    :Lj1 => IctoLj(1.32e-6),
    :Lj2 => IctoLj(1.32e-6),
    :Lg1 => 33e-12,
    :Lg2 => 33e-12,
    :Cj => 0.002e-12,
    :Cp => 2.0e-12,
    :Rs => 50.0,
    :Rp => 50.0,
    :RL => 2000.0,
    :K => 0.05,
)

Idc = 59.1e-5
Ip = 49.0e-6
wp = 2*pi*10.6201e9
ws = 2*pi*(5.0:0.001:5.6)*1e9
Npumpharmonics = (10,)
Nmodulationharmonics = (10,)

sources = [
      (mode=(0,),port=2,current=Idc),
      (mode=(1,),port=2,current=Ip),
  ]

sol = hbsolve(ws,(wp,),sources,Nmodulationharmonics, Npumpharmonics,
    circuit, circuitdefs;dc=true,threewavemixing=true,fourwavemixing=true)

p1 = plot(sol.linearized.w/(2*pi*1e9),
    10*log10.(abs2.(sol.linearized.S(
            outputmode=(0,),
            outputport=1,
            inputmode=(0,),
            inputport=1,
            freqindex=:),
    )),
    ylim=(-5,30),
    label="S11",
    xlabel="Signal Frequency (GHz)",
    legend=:bottomright,
    title="Scattering Parameters",
    ylabel="dB"
)

p2=plot(sol.linearized.w/(2*pi*1e9),
    sol.linearized.QE((0,),1,(0,),1,:)./sol.linearized.QEideal((0,),1,(0,),1,:),    
    title="Quantum efficiency",legend=false,
    ylabel="QE/QE_ideal",xlabel="Signal Frequency (GHz)");

p3=plot(sol.linearized.w/(2*pi*1e9),
    10*log10.(abs2.(sol.linearized.S(:,1,(0,),1,:)')),
    ylim=(-40,30),label="S21",
    xlabel="Signal Frequency (GHz)",
    legend=false,
    title="All idlers",
    ylabel="dB")


p4=plot(sol.linearized.w/(2*pi*1e9),
    1 .- sol.linearized.CM((0,),1,:),
    legend=false,title="Commutation \n relation error",
    ylabel="Commutation \n relation error",xlabel="Signal Frequency (GHz)");

plot(p1, p2, p3,p4,layout = (2, 2))

ws = 2*pi*(0.01:0.01:10.0)*1e9
currentvals = (0:0.25:150)*1e-5
outvals = zeros(Complex{Float64},length(ws),length(currentvals))
Ip=0.0

Npumpharmonics = (1,)
Nmodulationharmonics = (1,)

for (k,Idc) in enumerate(currentvals)
    sources = [
          (mode=(0,),port=2,current=Idc),
          (mode=(1,),port=2,current=Ip),
      ]
    sol = hbsolve(ws,(wp,),sources,Nmodulationharmonics, Npumpharmonics,
        circuit, circuitdefs;dc=true,threewavemixing=true,fourwavemixing=true)
    outvals[:,k]=sol.linearized.S((0,),1,(0,),1,:)
end

plot(
    currentvals/(1e-3),
    ws/(2*pi*1e9),
    10*log10.(abs2.(outvals)),
    seriestype=:heatmap,
    xlabel="bias current (mA)",
    ylabel="frequency (GHz)",
    title="Gain (dB) (pump off)",
)