JosephsonCircuits.jl is a high-performance frequency domain simulator for nonlinear circuits containing Josephson junctions, capacitors, inductors, mutual inductors, and resistors. JosephsonCircuits.jl simulates the frequency domain behavior using a variant [1] of nodal analysis [2] and the harmonic balance method [3-5] with an analytic Jacobian. Noise performance, quantified by quantum efficiency, is efficiently simulated through an adjoint method.
Frequency dependent circuit parameters are supported to model realistic impedance environments or dissipative components. Dissipation can be modeled by capacitors with an imaginary capacitance or frequency dependent resistors.
JosephsonCircuits.jl supports the following:
- Nonlinear simulations in which the user defines a circuit, the drive current, frequency, and number of harmonics and the code calculates the node flux or node voltage at each harmonic.
- Linearized simulations about the nonlinear operating point calculated above. This simulates the small signal response of a periodically time varying linear circuit and is useful for simulating parametric amplification and frequency conversion in the undepleted (strong) pump limit. Calculation of node fluxes (or node voltages) and scattering parameters of the linearized circuit [4-5].
- Linear simulations of linear circuits. Calculation of node fluxes (or node voltages) and scattering parameters.
- Calculation of symbolic capacitance and inverse inductance matrices.
As detailed in [6], we find excellent agreement with Keysight ADS simulations and Fourier analysis of time domain simulation performed by WRSPICE.
Warning: this package is under heavy development and there will be breaking changes. We will keep the examples updated to ease the burden of any breaking changes.
To install the latest release of the package, install Julia using Juliaup, start Julia, and enter the following command:
using Pkg
Pkg.add("JosephsonCircuits")
To install the development version, start Julia and enter the command:
using Pkg
Pkg.add(name="JosephsonCircuits",rev="main")
To run the examples below, you will need to install Plots.jl using the command:
Pkg.add("Plots")
If you get errors when running the examples, please try installing the latest version of Julia and updating to the latest version of JosephsonCircuits.jl by running:
Pkg.update()
Then check that you are running the latest version of the package with:
Pkg.status()
Generate a netlist using circuit components including capacitors C, inductors L, Josephson junctions described by the Josephson inductance Lj, mutual inductors described by the mutual coupling coefficient K, and resistors R. See the SPICE netlist format, docstrings, and examples below for usage. Run the harmonic balance analysis using hbnlsolve to solve a nonlinear system at one operating point, hblinsolve to solve a linear (or linearized) system at one or more frequencies, or hbsolve to run both analyses. Add a question mark ? in front of a function to access the docstring.
A driven nonlinear LC resonator.
using JosephsonCircuits
using Plots
@variables R Cc Lj Cj
circuit = [
("P1","1","0",1),
("R1","1","0",R),
("C1","1","2",Cc),
("Lj1","2","0",Lj),
("C2","2","0",Cj)]
circuitdefs = Dict(
Lj =>1000.0e-12,
Cc => 100.0e-15,
Cj => 1000.0e-15,
R => 50.0)
ws = 2*pi*(4.5:0.001:5.0)*1e9
wp = (2*pi*4.75001*1e9,)
Ip = 0.00565e-6
sources = [(mode=(1,),port=1,current=Ip)]
Npumpharmonics = (16,)
Nmodulationharmonics = (8,)
@time jpa = hbsolve(ws, wp, sources, Nmodulationharmonics,
Npumpharmonics, circuit, circuitdefs)
plot(
jpa.linearized.w/(2*pi*1e9),
10*log10.(abs2.(
jpa.linearized.S(
outputmode=(0,),
outputport=1,
inputmode=(0,),
inputport=1,
freqindex=:
),
)),
xlabel="Frequency (GHz)",
ylabel="Gain (dB)",
) 0.003080 seconds (57.81 k allocations: 6.391 MiB)
Circuit parameters from here.
using JosephsonCircuits
using Plots
@variables Rleft Rright Cg Lj Cj Cc Cr Lr
circuit = Tuple{String,String,String,Num}[]
# port on the input side
push!(circuit,("P$(1)_$(0)","1","0",1))
push!(circuit,("R$(1)_$(0)","1","0",Rleft))
Nj=2048
pmrpitch = 4
#first half cap to ground
push!(circuit,("C$(1)_$(0)","1","0",Cg/2))
#middle caps and jj's
push!(circuit,("Lj$(1)_$(2)","1","2",Lj))
push!(circuit,("C$(1)_$(2)","1","2",Cj))
j=2
for i = 2:Nj-1
if mod(i,pmrpitch) == pmrpitch÷2
# make the jj cell with modified capacitance to ground
push!(circuit,("C$(j)_$(0)","$(j)","$(0)",Cg-Cc))
push!(circuit,("Lj$(j)_$(j+2)","$(j)","$(j+2)",Lj))
push!(circuit,("C$(j)_$(j+2)","$(j)","$(j+2)",Cj))
#make the pmr
push!(circuit,("C$(j)_$(j+1)","$(j)","$(j+1)",Cc))
push!(circuit,("C$(j+1)_$(0)","$(j+1)","$(0)",Cr))
push!(circuit,("L$(j+1)_$(0)","$(j+1)","$(0)",Lr))
# increment the index
j+=1
else
push!(circuit,("C$(j)_$(0)","$(j)","$(0)",Cg))
push!(circuit,("Lj$(j)_$(j+1)","$(j)","$(j+1)",Lj))
push!(circuit,("C$(j)_$(j+1)","$(j)","$(j+1)",Cj))
end
# increment the index
j+=1
end
#last jj
push!(circuit,("C$(j)_$(0)","$(j)","$(0)",Cg/2))
push!(circuit,("R$(j)_$(0)","$(j)","$(0)",Rright))
# port on the output side
push!(circuit,("P$(j)_$(0)","$(j)","$(0)",2))
circuitdefs = Dict(
Lj => IctoLj(3.4e-6),
Cg => 45.0e-15,
Cc => 30.0e-15,
Cr => 2.8153e-12,
Lr => 1.70e-10,
Cj => 55e-15,
Rleft => 50.0,
Rright => 50.0,
)
ws=2*pi*(1.0:0.1:14)*1e9
wp=(2*pi*7.12*1e9,)
Ip=1.85e-6
sources = [(mode=(1,),port=1,current=Ip)]
Npumpharmonics = (20,)
Nmodulationharmonics = (10,)
@time rpm = hbsolve(ws, wp, sources, Nmodulationharmonics,
Npumpharmonics, circuit, circuitdefs)
p1=plot(ws/(2*pi*1e9),
10*log10.(abs2.(rpm.linearized.S(
outputmode=(0,),
outputport=2,
inputmode=(0,),
inputport=1,
freqindex=:),
)),
ylim=(-40,30),label="S21",
xlabel="Signal Frequency (GHz)",
legend=:bottomright,
title="Scattering Parameters",
ylabel="dB")
plot!(ws/(2*pi*1e9),
10*log10.(abs2.(rpm.linearized.S((0,),1,(0,),2,:))),
label="S12",
)
plot!(ws/(2*pi*1e9),
10*log10.(abs2.(rpm.linearized.S((0,),1,(0,),1,:))),
label="S11",
)
plot!(ws/(2*pi*1e9),
10*log10.(abs2.(rpm.linearized.S((0,),2,(0,),2,:))),
label="S22",
)
p2=plot(ws/(2*pi*1e9),
rpm.linearized.QE((0,),2,(0,),1,:)./rpm.linearized.QEideal((0,),2,(0,),1,:),
ylim=(0,1.05),
title="Quantum efficiency",legend=false,
ylabel="QE/QE_ideal",xlabel="Signal Frequency (GHz)");
p3=plot(ws/(2*pi*1e9),
10*log10.(abs2.(rpm.linearized.S(:,2,(0,),1,:)')),
ylim=(-40,30),
xlabel="Signal Frequency (GHz)",
legend=false,
title="All idlers",
ylabel="dB")
p4=plot(ws/(2*pi*1e9),
1 .- rpm.linearized.CM((0,),2,:),
legend=false,title="Commutation \n relation error",
ylabel="Commutation \n relation error",xlabel="Signal Frequency (GHz)");
plot(p1, p2, p3, p4, layout = (2, 2)) 2.959010 seconds (257.75 k allocations: 2.392 GiB, 0.21% gc time)
Circuit parameters from here.
using JosephsonCircuits
using Plots
@variables Rleft Rright Lj Cg Cc Cr Lr Cj
weightwidth = 745
weight = (n,Nnodes,weightwidth) -> exp(-(n - Nnodes/2)^2/(weightwidth)^2)
Nj=2000
pmrpitch = 8
# define the circuit components
circuit = Tuple{String,String,String,Num}[]
# port on the left side
push!(circuit,("P$(1)_$(0)","1","0",1))
push!(circuit,("R$(1)_$(0)","1","0",Rleft))
#first half cap to ground
push!(circuit,("C$(1)_$(0)","1","0",Cg/2*weight(1-0.5,Nj,weightwidth)))
#middle caps and jj's
push!(circuit,("Lj$(1)_$(2)","1","2",Lj*weight(1,Nj,weightwidth)))
push!(circuit,("C$(1)_$(2)","1","2",Cj/weight(1,Nj,weightwidth)))
j=2
for i = 2:Nj-1
if mod(i,pmrpitch) == pmrpitch÷2
# make the jj cell with modified capacitance to ground
push!(circuit,("C$(j)_$(0)","$(j)","$(0)",(Cg-Cc)*weight(i-0.5,Nj,weightwidth)))
push!(circuit,("Lj$(j)_$(j+2)","$(j)","$(j+2)",Lj*weight(i,Nj,weightwidth)))
push!(circuit,("C$(j)_$(j+2)","$(j)","$(j+2)",Cj/weight(i,Nj,weightwidth)))
#make the pmr
push!(circuit,("C$(j)_$(j+1)","$(j)","$(j+1)",Cc*weight(i-0.5,Nj,weightwidth)))
push!(circuit,("C$(j+1)_$(0)","$(j+1)","$(0)",Cr))
push!(circuit,("L$(j+1)_$(0)","$(j+1)","$(0)",Lr))
# increment the index
j+=1
else
push!(circuit,("C$(j)_$(0)","$(j)","$(0)",Cg*weight(i-0.5,Nj,weightwidth)))
push!(circuit,("Lj$(j)_$(j+1)","$(j)","$(j+1)",Lj*weight(i,Nj,weightwidth)))
push!(circuit,("C$(j)_$(j+1)","$(j)","$(j+1)",Cj/weight(i,Nj,weightwidth)))
end
# increment the index
j+=1
end
#last jj
push!(circuit,("C$(j)_$(0)","$(j)","$(0)",Cg/2*weight(Nj-0.5,Nj,weightwidth)))
push!(circuit,("R$(j)_$(0)","$(j)","$(0)",Rright))
push!(circuit,("P$(j)_$(0)","$(j)","$(0)",2))
circuitdefs = Dict(
Rleft => 50.0,
Rright => 50.0,
Lj => IctoLj(1.75e-6),
Cg => 76.6e-15,
Cc => 40.0e-15,
Cr => 1.533e-12,
Lr => 2.47e-10,
Cj => 40e-15,
)
ws=2*pi*(1.0:0.1:14)*1e9
wp=(2*pi*7.9*1e9,)
Ip=1.1e-6
sources = [(mode=(1,),port=1,current=Ip)]
Npumpharmonics = (20,)
Nmodulationharmonics = (10,)
@time floquet = hbsolve(ws, wp, sources, Nmodulationharmonics,
Npumpharmonics, circuit, circuitdefs)
p1=plot(ws/(2*pi*1e9),
10*log10.(abs2.(floquet.linearized.S((0,),2,(0,),1,:))),
ylim=(-40,30),label="S21",
xlabel="Signal Frequency (GHz)",
legend=:bottomright,
title="Scattering Parameters",
ylabel="dB")
plot!(ws/(2*pi*1e9),
10*log10.(abs2.(floquet.linearized.S((0,),1,(0,),2,:))),
label="S12",
)
plot!(ws/(2*pi*1e9),
10*log10.(abs2.(floquet.linearized.S((0,),1,(0,),1,:))),
label="S11",
)
plot!(ws/(2*pi*1e9),
10*log10.(abs2.(floquet.linearized.S((0,),2,(0,),2,:))),
label="S22",
)
p2=plot(ws/(2*pi*1e9),
floquet.linearized.QE((0,),2,(0,),1,:)./floquet.linearized.QEideal((0,),2,(0,),1,:),
ylim=(0.99,1.001),
title="Quantum efficiency",legend=false,
ylabel="QE/QE_ideal",xlabel="Signal Frequency (GHz)");
p3=plot(ws/(2*pi*1e9),
10*log10.(abs2.(floquet.linearized.S(:,2,(0,),1,:)')),
ylim=(-40,30),label="S21",
xlabel="Signal Frequency (GHz)",
legend=false,
title="All idlers",
ylabel="dB")
p4=plot(ws/(2*pi*1e9),
1 .- floquet.linearized.CM((0,),2,:),
legend=false,title="Commutation \n relation error",
ylabel="Commutation \n relation error",xlabel="Signal Frequency (GHz)");
plot(p1, p2, p3,p4,layout = (2, 2)) 2.079267 seconds (456.63 k allocations: 1.997 GiB, 0.48% gc time)
Dissipation due to capacitors with dielectric loss, parameterized by a loss tangent. Run the above code block to define the circuit then run the following:
results = []
tandeltas = [1.0e-6,1.0e-3, 2.0e-3, 3.0e-3]
for tandelta in tandeltas
circuitdefs = Dict(
Rleft => 50,
Rright => 50,
Lj => IctoLj(1.75e-6),
Cg => 76.6e-15/(1+im*tandelta),
Cc => 40.0e-15/(1+im*tandelta),
Cr => 1.533e-12/(1+im*tandelta),
Lr => 2.47e-10,
Cj => 40e-15,
)
wp=(2*pi*7.9*1e9,)
ws=2*pi*(1.0:0.1:14)*1e9
Ip=1.1e-6*(1+125*tandelta)
sources = [(mode=(1,),port=1,current=Ip)]
Npumpharmonics = (20,)
Nmodulationharmonics = (10,)
@time floquet = hbsolve(ws, wp, sources, Nmodulationharmonics,
Npumpharmonics, circuit, circuitdefs)
push!(results,floquet)
end
p1 = plot(title="Gain (S21)")
for i = 1:length(results)
plot!(ws/(2*pi*1e9),
10*log10.(abs2.(results[i].linearized.S((0,),2,(0,),1,:))),
ylim=(-60,30),label="tanδ=$(tandeltas[i])",
legend=:bottomleft,
xlabel="Signal Frequency (GHz)",ylabel="dB")
end
p2 = plot(title="Quantum Efficiency")
for i = 1:length(results)
plot!(ws/(2*pi*1e9),
results[i].linearized.QE((0,),2,(0,),1,:)./results[i].linearized.QEideal((0,),2,(0,),1,:),
ylim=(0.6,1.05),legend=false,
title="Quantum efficiency",
ylabel="QE/QE_ideal",xlabel="Signal Frequency (GHz)")
end
p3 = plot(title="Reverse Gain (S12)")
for i = 1:length(results)
plot!(ws/(2*pi*1e9),
10*log10.(abs2.(results[i].linearized.S((0,),1,(0,),2,:))),
ylim=(-10,1),legend=false,
xlabel="Signal Frequency (GHz)",ylabel="dB")
end
p4 = plot(title="Commutation \n relation error")
for i = 1:length(results)
plot!(ws/(2*pi*1e9),
1 .- results[i].linearized.CM((0,),2,:),
legend=false,
ylabel="Commutation\n relation error",xlabel="Signal Frequency (GHz)")
end
plot(p1, p2, p3,p4,layout = (2, 2)) 3.815835 seconds (470.00 k allocations: 2.303 GiB, 0.22% gc time)
3.800166 seconds (470.59 k allocations: 2.310 GiB, 0.29% gc time)
3.824690 seconds (470.75 k allocations: 2.317 GiB, 0.19% gc time)
3.838721 seconds (470.75 k allocations: 2.317 GiB, 0.18% gc time)
Simulations of the linearized system can be effectively parallelized, so we suggest starting Julia with the number of threads equal to the number of physical cores. See the Julia documentation for the procedure. Check how many threads you are using by calling Threads.nthreads().
- Andrew J. Kerman "Efficient numerical simulation of complex Josephson quantum circuits" arXiv:2010.14929 (2020)
- Jiří Vlach and Kishore Singhal "Computer Methods for Circuit Analysis and Design" 2nd edition, Springer New York, NY (1993)
- Stephen A. Maas "Nonlinear Microwave and RF Circuits" 2nd edition, Artech House (1997)
- José Carlos Pedro, David E. Root, Jianjun Xu, and Luís Cótimos Nunes. "Nonlinear Circuit Simulation and Modeling: Fundamentals for Microwave Design" The Cambridge RF and Microwave Engineering Series, Cambridge University Press (2018)
- David E. Root, Jan Verspecht, Jason Horn, and Mihai Marcu. "X-Parameters: Characterization, Modeling, and Design of Nonlinear RF and Microwave Components" The Cambridge RF and microwave engineering series, Cambridge University Press (2013)
- Kaidong Peng, Rick Poore, Philip Krantz, David E. Root, and Kevin P. O'Brien "X-parameter based design and simulation of Josephson traveling-wave parametric amplifiers for quantum computing applications" IEEE International Conference on Quantum Computing & Engineering (QCE22) (2022)
The motivation for developing this package is to simulate the gain and noise performance of ultra low noise amplifiers for quantum computing applications such as the Josephson traveling-wave parametric amplifier, which have thousands of linear and nonlinear circuit elements.
We prioritize speed (including compile time and time to first use), simplicity, and scalability.
- Design optimization.
- More nonlinear components such as kinetic inductors.
- Time domain simulations.
- Xyce.jl provides a wrapper for Xyce, the open source parallel circuit simulator from Sandia National Laboratories which can perform time domain and harmonic balance method simulations.
- NgSpice.jl and LTspice.jl provide wrappers for NgSpice and LTspice, respectively.
- ModelingToolkit.jl supports time domain circuit simulations from scratch and using their standard library
- ACME.jl simulates electrical circuits in the time domain with an emphasis on audio effect circuits.
- Cedar EDA is a Julia-based commercial cloud service for circuit simulations.
- Keysight ADS, Cadence AWR, Cadence Spectre RF, and Qucs are capable of time and frequency domain analysis of nonlinear circuits. WRSPICE performs time domain simulations of Josephson junction containing circuits and frequency domain simulations of linear circuits.
We gratefully acknowledge funding from the AWS Center for Quantum Computing and the MIT Center for Quantum Engineering (CQE).
