CPB in phase basis

In the phase - basis solving the Schrödinger Equation for a  Cooper pair box superconducting electronic circuit leads to a second - order differential equation, which can be solved in terms of Mathieu functions:

DSolve[-Ec ϕ''[δ]-Ej Cos[δ]ϕ[δ]-a ϕ[δ]==0,ϕ[δ],δ]

The Mathieu functions are only periodic in δ for certain values of the value a, the Mathieu characteristic values a(r,-2(Ej/Ec)). This characteristic value determines the energy of the system, where r is an integer or rational number.

En[r_,Ec_,Ej_]:=(Ec/4) MathieuCharacteristicA[r,-2(Ej/Ec)]

Transforming back to the state vector Ψ we obtain the solution :

ψ[δ_,Ec_,Ej_,ng_,r_]=C[1] Exp[I ng δ]MathieuC[((4 En[r,Ej,Ec])/Ec),-((2 Ej)/Ec),(δ/2)]+C[2]  Exp[I ng δ]MathieuS[((4 En[r,Ej,Ec])/Ec),-((2 Ej)/Ec),(δ/2)]

From Floquet theory, Mathieu Functions are periodic only if r is integer or rational, which one can see by writing them in the form 
(Subscript[M, C])(a,q,(δ/2))=exp[i r (δ/2)] mc((δ/2))
(Subscript[M, S])(a,q,(δ/2))=exp[i r (δ/2)] ms((δ/2))

with r(a,q) a real parameter and mc(z) and ms(z) 2π-periodic on δ.
This gives 

Ψ(δ,r,..)= (Subscript[C, 1])(e^(i ((Subscript[n, g])+(r/2))δ))mc((δ/2))+(Subscript[C, 2])(e^(i ((Subscript[n, g])+(r/2))δ))ms((δ/2))

Since the function we search shall satisfy the periodic boundary conditions ψ[δ] = ψ[δ + 2 π] this implies that ng+r/2 shall be integer, i.e. for each ng only a discrete set of r=(Subscript[r, k]) with integer k is possible. This yields a set of associated eigenfunctions (Subscript[Ψ, k]).

energies=Table[Tooltip[En[2(j-ng),1,0.1],ToString[j]],{j,-3,3}];

Plot[energies,{ng,-1,2},BaseStyle->{FontSize->14},Frame->True,FrameLabel->{"ng","(Subscript[E, n])"}]

Correct ordering :

tri[ng_]:=((2ng-Floor[2ng])*(1-2Mod[Floor[2ng],2])+Mod[Floor[2ng],2])/2
MathieuCharacteristicR[ng_,ib_]:=(ib+1)-Mod[ib+1,2]-2(((-1))^((ib+1)))tri[ng]

Plot with different EJ/EC ratios :

Manipulate[
Plot[Evaluate[Table[En[MathieuCharacteristicR[ng,k],1,EJCratio],{k,0,3}]],{ng,0,3},BaseStyle->{FontSize->14},Frame->True,FrameLabel->{"ng","(Subscript[E, n])"},ImageSize->550,Axes->False],
{{EJCratio,0.2},0.,50.}]

Look at the % variation of the 01 transition with gate charge :

TransitionP[ng_,EJCratio_]:=-100*(En[MathieuCharacteristicR[ng,1],1,EJCratio]-En[MathieuCharacteristicR[ng,0],1,EJCratio])/(En[MathieuCharacteristicR[0.001,1],1,EJCratio]-En[MathieuCharacteristicR[0.001,0],1,EJCratio])+100

Manipulate[
Plot[TransitionP[ng,EJCratio],{ng,0,3.},BaseStyle->{FontSize->14},Frame->True,PlotRange->{-5.,20.},FrameLabel->{"ng","% variation of transition 01"},ImageSize->550,Axes->False],
{{EJCratio,2.},1.,10.}]

CPB Two level approximation

Clear[EC,EJ,ng]
HCPB=(({{EC(1-2 ng), EJ}, {EJ, EC(2 ng-1)}}))

Find the eigenstates and eigenvalues of the 2-level Hamiltonian :

Eigs=Eigensystem[HCPB]/.EC->1

Manipulate[
Plot[{Eigs[[1,1]]/.EJ->0,Eigs[[1,2]]/.EJ->0,Eigs[[1,1]]/.EJ->EJ1,Eigs[[1,2]]/.EJ->EJ1},{ng,0,1},Frame->True,PlotStyle->{Directive[Black,Dashed],Directive[Black,Dashed],Blue,Red},Axes->{True,False},AxesStyle->LightGray],{{EJ1,0.1},0,1}]

© Stefan Filipp, Peter Leek 11/2009