ABSTRACT (original article): We present an exact solution for the linearized junction array modes of the superconducting qubit fluxonium in the absence of array disorder. This solution holds for arrays of any length and ground capacitance and for both differential and grounded devices. Array mode energies are determined by roots of convex combinations of Chebyshev polynomials, and their spatial profiles are modified plane waves. We also provide a simple, approximate solution which estimates array mode properties over a wide range of circuit parameters and an accompanying Mathematica file that implements both the exact and approximate solutions. CITATION (original article): Sorokanich, S., Hays, M., & Warrington, N. C. (2024). Exact and approximate fluxonium array modes. Physical Review B 110 (12), 125404. https://doi.org/10.1103/physrevb.110.125404
In this notebook we provide code implementing and checking results in our paper. This notebook is organized in the following way:
-In “Definitions”, many functions used throughout the rest of the notebook are defined. These are implemented in “Functions”. So that the user can easily see what they do, many of these functions are executed in “Illustrations.
-In “Implementations” several useful tools are coded.
-In “Exact solution”, the exact solution for the inverse capacitance eigensystem is implemented.
-In “Approximate solution”, two approximate solutions for the eigenvalues are implemented.
-In “Eigensystem from sparse matrices”, a fast method for numerically computing the inverse capacitance matrix eigensystem is implemented.
-In “Checking formulas”, many formulas of the text are checked with short code snippets.
Note:
-Throughout this notebook, we call the number of array junctions both “N” and “n”; In writing we use N, because it matches the text of the paper (and the usual terminology), but in code we call it “n” because in Mathematica “N” is a protected variable.
-Throughout this notebook, we call the number of array junctions both “N” and “n”; In writing we use N, because it matches the text of the paper (and the usual terminology), but in code we call it “n” because in Mathematica “N” is a protected variable.
Definitions
Definitions
Implementations
Implementations
Exact solution
Exact solution
In the section below we implement the exact solution of the inverse capacitance matrix, contained in Eqs. 29 & 30 of the paper. The code below implements the following logic:1. Charging energies are supplied (in units of h GHz)2. These are converted to capacitances , , , (in units of fF).3. The quantities ϵ, ϵ’, , (from Eqs. 20 & 27) are formed 4. The eigensystem of -Δ(ϵ,ϵ’) is computed.5. The eigenvectors of -Δ(ϵ,ϵ’) are identical to those of . The eigenvalues of -Δ(ϵ,ϵ’), called l, are related to those of , called λ, via the relation λ = . At the very end of the calculation the ls are converted to λs Computing the eigensystem of -Δ(ϵ,ϵ’) constitutes most of the calculation below. It proceeds in the following steps:1. The N roots of the polynomial (Eq. 29 of the paper ) are computed through an NSolve.2. The N roots of the polynomial are computed through an NSolve.3. 2N normalized vectors are formed from Eq. 30. 4. The locations of the N vectors with the correct symmetry are identified. 5. The lists are trimmed to those eigenvalues/eigenvectors with the correct symmetry.For convenience we plot the eigenvalues of at the end. We also compute the eigenvalues with direct diagonalization and print the difference between the result of the procedure above and direct diagonalization to show the methods agree.FINALLY & IMPORTANTLY: The method below is simply a demonstration. We do not claim it is the most efficient or wise way to implement our exact solution; our aim is to simply provide a framework that users can begin with. The code block below struggles to if N > 20, because Mathematica struggles to find the roots of high-order Chebyshev polynomials. If one desires to use our solution for large N, we suggest two alternatives. Either a more robust root-finding method should be implemented by the user, or one should leverage the sparse matrix representation of the inverse capacitance matrix. For convenience, we implement the second option below, in the section “Eigensystem from sparse matrices”.
a
C
b
C
a
C
g
b
C
b
ϵ
+
ϵ
-
-1
C
-1
C
-1
+l
a
C
a
C
g
P
+
P
-
-1
C
Block{ϵ,ϵp,ϵplus,ϵminus,n},(*numberofarrayjunctions*)n=16;(*examplechargingenergiesfromApp.Fofhttps://arxiv.org/pdf/1208.5747;unitsarehGHz*){Ea,Eb,Eag,Ebg}={1.0,3.7,1750,5};(*converttocapacitances;unitsarefF*){Ca,Cb,Cag,Cbg}={Cap[Ea],Cap[Eb],Cap[Eag],Cap[Ebg]};ϵ=;ϵp=;ϵplus=ϵ-ϵp;ϵminus=ϵ+ϵp;(*findallnrootsofthepolynomial*)rootsPLUS=l/.NSolve[SetPrecision[CharPoly[n,ϵplus,l],16]==0,l];(*findallnrootsofthepolynomial*)rootsMINUS=l/.NSolve[SetPrecision[CharPoly[n,ϵminus,l],16]==0,l];(*evaluateallnvectorsofEq.30&normalizethem*)vecsPLUS=Table,{i,Length@rootsPLUS};(*evaluateallnvectorsofEq.30&normalizethem*)vecsMINUS=Table,{i,Length@rootsMINUS};(*figureouttheindiciesofallevenvectors*)indPLUS={};Do(*testifvectoriseven*)If[Abs[vecsPLUS[[i]][[1]]-vecsPLUS[[i]][[n]]]<,AppendTo[indPLUS,i];,0];,{i,Length@vecsPLUS};(*keeponlytheeveneigenvalues/eigenvectorsandstoretheminvecsPLUS/valsPLUS*)tmpvecsPLUS=Table[vecsPLUS[[indPLUS[[i]]]],{i,Length@indPLUS}];tmpvalsPLUS=Table[rootsPLUS[[indPLUS[[i]]]],{i,Length@indPLUS}];vecsPLUS=tmpvecsPLUS;valsPLUS=tmpvalsPLUS;(*figureouttheindiciesofalloddvectors*)indMINUS={};Do(*testifvectorisodd*)If[Abs[vecsMINUS[[i]][[1]]+vecsMINUS[[i]][[n]]]<,AppendTo[indMINUS,i];,0];,{i,Length@vecsMINUS};(*keeponlytheoddeigenvalues/eigenvectorsandstoretheminvecsMINUS/valsMINUS*)tmpvecsMINUS=Table[vecsMINUS[[indMINUS[[i]]]],{i,Length@indMINUS}];tmpvalsMINUS=Table[rootsMINUS[[indMINUS[[i]]]],{i,Length@indMINUS}];vecsMINUS=tmpvecsMINUS;valsMINUS=tmpvalsMINUS;(*joinalleigenvalues/eigenvectorsintosinglelists*)(*theRiffle&Reverseoperations,andthecarefulorderingof"MINUS"and"PLUS"simplyputseverythinginorderofincreasingeigenvalue;notnecessary,butconvenient*)If[n/2∈Integers,valsALL=Reverse@Riffle[valsMINUS,valsPLUS];vecsALL=Reverse@Riffle[vecsMINUS,vecsPLUS];,valsALL=Reverse@Riffle[valsPLUS,valsMINUS];vecsALL=Reverse@Riffle[vecsPLUS,vecsMINUS];];(*formeigenvaluesfrom-Δ(ϵ,ϵ')eigenvalues*)CapMatEvals=Table,{i,Length@valsALL};(*plottheeigenvaluesoftheinversecapacitancematrix*)PrintListPlotCapMatEvals,AxesLabel->{"index","eval []"},PlotRange->All;(*aquickcrosscheckthatthecorrectvalueswereobtained*)c=C1[n,Ca]+C2[n,Cb,Cbg]+C3[n,Cag]+C4[n,Cag,Cbg];StandardMethod=Sort[Eigenvalues@Inverse[c]];Print["Diagonalization Eigenvalues - Formula Eigenvalues = "];Print[Chop[CapMatEvals-StandardMethod]];;
Cag
Cbg
Cb+Cbg
2Cb+Cbg
Cag
Cbg
Cb
2Cb+Cbg
P
+
P
-
v
+
eigenvector[n,ϵplus,rootsPLUS[[i]]]
1/2
(eigenvector[n,ϵplus,rootsPLUS[[i]]].eigenvector[n,ϵplus,rootsPLUS[[i]]])
v
-
eigenvector[n,ϵminus,rootsMINUS[[i]]]
1/2
(eigenvector[n,ϵminus,rootsMINUS[[i]]].eigenvector[n,ϵminus,rootsMINUS[[i]]])
-12
10
-12
10
-1
C
1
Ca+Cag/valsALL[[i]]
-1
C
-1
fF
Approximate solution
Approximate solution
In this section we implement eigenvalue approximation schemes. The first scheme, summarized in Eq. 38-39 in the text, is coded below as “Approximation1”.
The second scheme, which leverages Eq. 44, is “Approximation2”. We provide an example of the method below as well.
The second scheme, which leverages Eq. 44, is “Approximation2”. We provide an example of the method below as well.
This is the approximation scheme used in the paper. returns a list of approximate eigenvalues, ordered from lowest to highest
This is the alternative approximation scheme NOT used in the paper; it less accurate at small N, but has the correct large N behavior. The following returns a list of approximate eigenvalues, ordered from lowest to highest
Numerical diagonalization from sparse matrices
Numerical diagonalization from sparse matrices
At the end of the code block below we print the computed eigenvalues and the lowest three eigenvectors.
Note:
Using these matrix methods, very large arrays can be diagonalized. In the example below we have set N=10,000, which runs within seconds. In contrast, forming the dense matrix C , then diagonalizing it, is prohibitively expensive.
Using these matrix methods, very large arrays can be diagonalized. In the example below we have set N=10,000, which runs within seconds. In contrast, forming the dense matrix C , then diagonalizing it, is prohibitively expensive.
Again, we do not claim this code is optimized in any way.
Checking formulas
Checking formulas
Verifying Eq. 29 & Eq. 30 (most important check)
Verifying Eq. 29 & Eq. 30 (most important check)
This is the most important check in this notebook. It shows our exact solution holds. Below we show that the eigensystem of -Δ(ϵ,ϵ’) is given by Eqs. 29 & 30. We do this by
1. First computing the eigensystem numerically
2. Using the formulae
2. Using the formulae
N even
N even
N odd
N odd
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Exact and approximate fluxonium array modes
by Stephen Sorokanich, Max Hays & Neill C. Warrington
Wolfram Community, STAFF PICKS, September 19, 2024
https://community.wolfram.com/groups/-/m/t/3277880
by Stephen Sorokanich, Max Hays & Neill C. Warrington
Wolfram Community, STAFF PICKS, September 19, 2024
https://community.wolfram.com/groups/-/m/t/3277880