Harmonic Oscillator, Current biased Josephson junction

Some Functions  and Definitions:

(* general pure state in polar coordinates *)
State[θ_,ϕ_]:={Cos[θ/2],Sin[θ/2]Exp[I ϕ]}
(*normalized state with real first argument *)
StateNorm[x_List]:=Module[{arg1=Arg[x[[1]]],arg2=Arg[x[[2]]],normx=x/Norm[x],θ,ϕ},
θ=2Chop@ArcCos[normx[[1]]Exp[-I arg1]];
ϕ=Mod[arg2-arg1,2π];
{Cos[θ/2],Sin[θ/2]Exp[I ϕ]}]
(* outputs the polar coordinate from a state vector *)
blochcoord[state_]:=Module[{arg1=Arg[state[[1]]],arg2=Arg[state[[2]]],normx=state/Norm[state],θ,ϕ},
θ=Abs[2Chop@ArcCos[normx[[1]]Exp[-I arg1]]];
ϕ=Abs@Mod[arg2-arg1,2π];
{θ,ϕ}]

(* define a bloch sphere graph *)
With[{r=1},xy=ParametricPlot3D[{r Cos[θ],r Sin[θ],0},{θ,0,2 π},PlotStyle->Dashed,Boxed->False,Axes->False];
yz=ParametricPlot3D[{0,r Cos[θ],r Sin[θ]},{θ,0,2 π},PlotStyle->Dashed,Boxed->False,Axes->False];
xz=ParametricPlot3D[{r Cos[θ],0,r Sin[θ]},{θ,0,2 π},PlotStyle->Dashed,Boxed->False,Axes->False];
bloch=Show[Graphics3D[{Green,Opacity[.25],Sphere[{0,0,0},r],Opacity[1],Black,Thick,Line[{{0,r,0},{0,-r,0}}],Line[{{0,0,r},{0,0,-r}}],Line[{{r,0,0},{-r,0,0}}],Text[Style["z+",FontSize->16,FontWeight->Bold,FontColor->GrayLevel[.3]],{0,0,1.3 r}],Text[Style["z-",FontSize->16,FontWeight->Bold,FontColor->GrayLevel[.3]],{0,0,-1.3 r}],Text[Style["x+",FontSize->16,FontWeight->Bold,FontColor->GrayLevel[.3]],{1.3 r,0,0}],Text[Style["x-",FontSize->16,FontWeight->Bold,FontColor->GrayLevel[.3]],{-1.3 r,0,0}],Text[Style["y+",FontSize->16,FontWeight->Bold,FontColor->GrayLevel[.3]],{0,1.3 r,0}],Text[Style["y-",FontSize->16,FontWeight->Bold,FontColor->GrayLevel[.3]],{0,-1.3 r,0}]}],xy,xz,yz,Boxed->False]];
(* generate vector for bloch sphere representation *)
Options[blochvec]={color->Blue};
blochvec[{θ_,ϕ_},OptionsPattern[]]:=Graphics3D[{OptionValue[color],Thickness[0.02],Line[{{0,0,0},{Sin[θ]Cos[ϕ],Sin[θ]Sin[ϕ],Cos[θ]}}]}]

(* draw vector for bloch sphere taking an arbitrary pure input state *)
blochvecstate[state_]:=blochvec@blochcoord[state]
(*bloch vector for mixed states *)
blochvecmsphere[r_,θ_,ϕ_]:=Graphics3D[{Blue,Thickness[0.02],Line[{{0,0,0},r{Sin[θ]Cos[ϕ],Sin[θ]Sin[ϕ],Cos[θ]}}]}]
blochvecmsphere[{r_,θ_,ϕ_}]:=blochvec[r,θ,ϕ]
(* cartesian coordinates *)
blochvecmcart[x_,y_,z_]:=Graphics3D[{Blue,Thickness[0.02],Line[{{0,0,0},{x,y,z}}]}]
blochvecmcart[{x_,y_,z_}]:=blochvecmcart[x,y,z]
(* show vector with sphere *)
vectoronsphere[vec_]:=Show[bloch,vec,SphericalRegion->True]
(* Bloch-sphere coordinates from arbitrary state *)
MatrixToBloch[x_]:={Tr[σx.x],Tr[σy.x],Tr[σz.x]}

id=(({{1, 0}, {0, 1}}));σx=(({{0, 1}, {1, 0}}));σy=(({{0, -I}, {I, 0}}));σz=(({{1, 0}, {0, -1}}));
tps[x_]:=Outer[Times,Conjugate[x],x]
tps[x_,y_]:=Outer[Times,Conjugate[x],y]
(* groundstate *) 
gs={1,0};
(* excited state *)
es={0,1};
(* superpositions state *)
xps={1,1}/(Sqrt[2]);
xms={1,-1}/(Sqrt[2]);
yps={1,I}/(Sqrt[2]);
yms={1,-I}/(Sqrt[2]);

TensorProduct[a_,b_]:=Table[(*a,b are nxn and mxm-martices*)a[[Ceiling[s/Length[b]],Ceiling[t/Length[b]]]]*b[[s-Floor[(s-1)/Length[b]]*Length[b],t-Floor[(t-1)/Length[b]]*Length[b]]],{s,1,Length[a]*Length[b]},{t,1,Length[a]*Length[b]}]

gg={1,0,0,0};ge={0,1,0,0};eg={0,0,1,0};ee={0,0,0,1};

Urot[α_,t_,θ_,ϕ_]:=(({{Cos[α t/2] - I Cos[θ]Sin[α t /2], -I Exp[-I ϕ] Sin[θ] Sin[α t/2]}, {(-I Exp[I ϕ]) Sin[θ] Sin[α t/2], Cos[α t/2] + I Cos[θ]Sin[α t /2]}}))

RX[α_]:=Urot[α,1,π/2,0]
RY[α_]:=Urot[α,1,π/2,π/2]
RZ[α_]:=Urot[α,1,0,0]

ZZ=TensorProduct[σz,σz];
RX[α_]:=Urot[α,1,π/2,0]
RY[α_]:=Urot[α,1,π/2,π/2]
RZ[α_]:=Urot[α,1,0,0]
RX1[α_]:=TensorProduct[Urot[α,1,π/2,0],(({{1, 0}, {0, 1}}))]
RY1[α_]:=TensorProduct[Urot[α,1,π/2,π/2],(({{1, 0}, {0, 1}}))]
RZ1[α_]:=TensorProduct[Urot[α,1,0,0],(({{1, 0}, {0, 1}}))]
RX2[α_]:=TensorProduct[(({{1, 0}, {0, 1}})),Urot[α,1,π/2,0]]
RY2[α_]:=TensorProduct[(({{1, 0}, {0, 1}})),Urot[α,1,π/2,π/2]]
RZ2[α_]:=TensorProduct[(({{1, 0}, {0, 1}})),Urot[α,1,0,0]]

Frequency response of a damped Harmonic Oscillator as the Fourier transform of the energy decay rate

RLC - Circuit :

(* solution for underdamped oscillator *)

IL[t_]:=I0 Exp[((-α t)/2)+I ωr t]/;t≥0;
IL[t_]:=I0/;t<0;

Plot[Re[IL[t]/.{α->0.1,ωr->3,I0->1}],{t,-3,10},PlotRange->All]

Fourier Transform of I (t) :

Clear[ω]

Iω[ω_]=Integrate[I0 Exp[((-α t)/2)+I ωr t]Exp[-I ω t],{t,0,∞},Assumptions->{α>0,ωr∈Reals,Im[ω]==0}]

ISqr[ω_,α_,ωr_,I0_]=Iω[ω]Conjugate[Iω[ω]]//FullSimplify[#,Assumptions->{α>0,ωr∈Reals,Im[ω]==0,Im[I0]==0}]&

Plot[ISqr[o,1,10,1],{o,5,15},PlotRange->All,Frame->True,FrameLabel->{"ω","|I(ω)((|)^2)"},BaseStyle->{FontSize->14},Axes->False]

Full width at half maximum (FWHM) :

Solve[ISqr[ω,a,ωr,I0]==ISqr[ωr,a,ωr,I0]/2,ω]
FWHM=(ω/.Last[%])-(ω/.First[%])//Simplify

Plot[ISqr[o,1,10,1],{o,5,15},PlotRange->All,Frame->True,FrameLabel->{"ω","|I(ω)((|)^2)"},BaseStyle->{FontSize->14},Axes->False,GridLines->{{10-.5,10+0.5},None},Epilog->{Text["a",{10,ISqr[10,1,10,1]/2}]}]

Time-dependent Harmonic Oscillator

Time-independent harmonic oscillator

division=200.;
mass=2;

DiffMatrix[m_,dx_,Vx_]:=(DiagonalMatrix[-((ℏ^2)/(2m))((-2)/(dx^2))+Vx]+DiagonalMatrix[Table[-((ℏ^2)/(2m))(1/(dx^2)),{Length[Vx]-1}],1]+DiagonalMatrix[Table[-((ℏ^2)/(2m))(1/(dx^2)),{Length[Vx]-1}],-1])/.ℏ->1

Vfunc[k_,x_]:=(k/2) (x^2)

Varray[V_,boundary_,dx_]:=Table[Vfunc[V,x],{x,-boundary,boundary,dx}];

harmpotential=Varray[60,5,2 5/division];

Set up the Matrix for the Problem and calculate the Eigenfunctions and Eigenvectors:

MProblem=DiffMatrix[mass,2 5/division,harmpotential];
answers=Eigensystem[MProblem];
{evals,evects}=Sort[Transpose[answers],#1[[1]]<#2[[1]]&]//Transpose;
stateenergies=Take[evals,10]

Plot the Eigenvectors on top of the Eigenenergies :

potplot=ListPlot[harmpotential,Filling->1.1Min[harmpotential],Joined->True,Frame->True,PlotRange->{0,60},GridLines->{None,stateenergies},GridLinesStyle->Directive[Thickness[0.002],Red],Axes->False];
ScaleFactor=10;
stateplot=ListPlot[Table[ScaleFactor evects[[i]]+evals[[i]],{i,1,10}],Joined->True,Filling->Table[n->evals[[n]],{n,1,10}],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]],PlotRange->All,ImageSize->600,Frame->True,Axes->False];
Show[stateplot,potplot,PlotRange->{Automatic,{0,60}}]

Time-dependent harmonic oscillator :

Frequency :

ω=(Sqrt[D[Vfunc[60,x],x,x]/mass]);
tfinal=2π/ω;

Start in the groundstate :

MProblem=DiffMatrix[mass,2 5/division,harmpotential];
uinitial=evects[[1]];

Energy of the ground state :

{uinitial.MProblem.uinitial,N[ω/2]}

(Subscript[Ω, 0])=(ω/2) in the ground state (vacuum energy)

NDSolve[{-I MProblem.u[t]==u'[t],u[0]==uinitial},u[t],{t,0,2tfinal}]

answertdepground[t_]=u[t]/.First@@NDSolve[{-I MProblem.u[t]==u'[t],u[0]==uinitial},u[t],{t,0,2tfinal}]

ListPlot[Abs[answertdepground[6 tfinal]],PlotRange->All]

Animate[Show[potplot,ListPlot[ evals[[1]]+ScaleFactor Abs[answertdepground[t]],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]]],PlotRange->{-5,30}],{t,0,2tfinal},AnimationRunning->False]

Start in the 1st excited state :

expected frequency (Subscript[Ω, 1])=((3ω)/2)

MProblem=DiffMatrix[mass,2 5/division,harmpotential];
uinitial=evects[[2]];
answertdep1[t_]=u[t]/.First@@NDSolve[{-I MProblem.u[t]==u'[t],u[0]==uinitial},u[t],{t,0,2tfinal}];

Animate[Show[potplot,ListPlot[ evals[[2]]+20 Abs[answertdep1[t]],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]]],PlotRange->{0,30}],{t,0,2tfinal},AnimationRunning->False]

Start in the 2nd excited state :

(Subscript[Ω, 2])=((5ω)/2)

MProblem=DiffMatrix[mass,2 5/division,harmpotential];
uinitial=evects[[3]];
answertdep2[t_]=u[t]/.First@@NDSolve[{-I MProblem.u[t]==u'[t],u[0]==uinitial},u[t],{t,0,2tfinal}];

Animate[Show[potplot,ListPlot[{evals[[3]]+20 Abs[answertdep2[t]],evals[[2]]+20 Abs[answertdep1[t]],evals[[1]]+20 Abs[answertdepground[t]]},PlotStyle->{Black,Red,Green}],PlotRange->{-5,30}],{t,0,2tfinal},AnimationRunning->False]

Start in superposition state :

MProblem=DiffMatrix[mass,2 5/division,harmpotential];
uinitial=1/(Sqrt[2])(evects[[2]]+evects[[1]]);
answertdepsup[t_]=u[t]/.First@@NDSolve[{-I MProblem.u[t]==u'[t],u[0]==uinitial},u[t],{t,0,2tfinal}];

Animate[Show[potplot,ListPlot[ evals[[1]]+ScaleFactor Abs[answertdepsup[t]],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]]],PlotRange->{{50,150},{0,20}}],{t,0,2tfinal},AnimationRate->0.4,AnimationRunning->False]

Particles bounces in between the potential walls.

Start in a localized (Gaussian) state :

PDF[NormalDistribution[0,0.05],x]

ugauss=Table[PDF[NormalDistribution[0,0.05],x],{x,-5,5,2 5/division}];

Normalize to |⟨Ψ|Ψ⟩(|^2)=1:

normsqrt=(Sqrt[ugauss\[Conjugate].ugauss])

uinitial=ugauss/normsqrt;

ListPlot[uinitial,PlotRange->All,Joined->True,Frame->True,FrameLabel->{"x","|Ψ(x)((|)^2)"},BaseStyle->{FontSize->14}]

Energy of the gaussian state:

gaussenergy=uinitial\[Conjugate].MProblem.uinitial

MProblem=DiffMatrix[mass,2 5/division,harmpotential];
answertdeploc[t_]=u[t]/.First@@NDSolve[{-I MProblem.u[t]==u'[t],u[0]==uinitial},u[t],{t,0,tfinal}];

Animate[Show[potplot,ListPlot[gaussenergy+10Abs[answertdeploc[t]],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]]],PlotRange->{{60,150},{40,60}}],{t,0,tfinal},AnimationRunning->False]

Start in a localized state at x≠ 0 :

uinitial=#/Norm[#]&[Table[PDF[NormalDistribution[-.3,0.4],x],{x,-5,5,2 5/division}]];
ListPlot[uinitial,PlotRange->All]

Energy of the displaced gaussian state:

gaussenergy2=uinitial\[Conjugate].MProblem.uinitial

MProblem=DiffMatrix[mass,2 5/division,harmpotential];
answertdeploc2[t_]=u[t]/.First@@NDSolve[{-I MProblem.u[t]==u'[t],u[0]==uinitial},u[t],{t,0,tfinal}];

Animate[Show[potplot,ListPlot[gaussenergy2+10Abs[answertdeploc2[t]],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]]],PlotRange->{{50,150},{4,10}}],{t,0,tfinal},AnimationRunning->False]

Current Biased Josephson junction :

RCSJ (resistively and capacitively shunted junction) model

Tilted Washboard Potential

VCurrentBiased[δ_,Ib_,Ic_]=((-Ic Φ0)/(2π ))(Cos[δ]+(Ib/Ic)δ)

Manipulate[Plot[VCurrentBiased[δ,Ib,1]/.Φ0->1,{δ,0,10π},PlotRange->{-.2-2π Ib,0.2},Frame->True,FrameLabel->{"δ","U(δ)"},BaseStyle->{FontSize->14}],{Ib,0,1}]

Find the height of the potential barrier U0 :

D[VCurrentBiased[δ,Ib,Ic],δ]
Solve[D[VCurrentBiased[δ,Ib,Ic],δ]==0,δ]

has nodal points at δ=(2k+1)π-ArcSin[(Ib/Ic)] and (2k)π+ArcSin[(Ib/Ic)], k=0,1,2,...

Plot[{VCurrentBiased[x,0.1,1]/.Φ0->1,},{x,0,20},GridLines->{{3 π-ArcSin[.1],2π+ArcSin[.1]},None},PlotPoints->30]

V0CB[Ib_,Ic_]=VCurrentBiased[3π-ArcSin[(Ib/Ic)],Ib,Ic]-VCurrentBiased[2π+ArcSin[(Ib/Ic)],Ib,Ic]//FullSimplify

Oscillation Freqency (Subscript[ω, 0 ])given by  (Subscript[ω, 0])=(Sqrt[((U''[(Subscript[δ, min])])/m)]), m=C(((Φ0/2π))^2):

OscFrequency[Cp_,Ib_,Ic_]=(Sqrt[D[VCurrentBiased[x,Ib,Ic],x,x]/(Cp ((Φ0)^2)/(((2π))^2))])/.x->ArcSin[(Ib/Ic)]//FullSimplify[#,Assumptions->{Φ0>0}]&

with the plasma frequency (Subscript[ω, p])=(Sqrt[((2π (Subscript[I, c]))/((Subscript[Φ, 0])C))])(oscillation frequency of the 'phase'-particle at the bottom of the well)

Plot[V0CB[x,1]/.Φ0->1,{x,0,2},Frame->True,FrameLabel->{"Ib/Ic","U0(Ib/Ic)"},BaseStyle->{FontSize->14}]

Set up the tables and matrices for the problem :

division=100;
mass=0.6;
Icconst=10;

VCBpotential[Ib_]=Table[VCurrentBiased[x,Ib,Icconst]/.Φ0->1,{x,0,7.5 2π,2 5/division}];

potplottilt=ListPlot[VCBpotential[0.1],Filling->1.1Min[VCBpotential[0.1]],Joined->True,Frame->True,PlotRange->{-4,4},GridLinesStyle->Directive[Thickness[0.002],Red],Axes->False]

Set up the matrix for the problem and calculate the Eigenfunctions and Eigenvectors:

MProblemCB=DiffMatrix[mass,2 5./division,VCBpotential[0.1]];
answersCB=Eigensystem[MProblemCB];
{evalsCB,evectsCB}=Sort[Transpose[answersCB],#1[[1]]<#2[[1]]&]//Transpose;
stateenergiesCB=Take[evalsCB,13]

ScaleFactor=2;
stateplotCB=ListPlot[Table[ScaleFactor evectsCB[[i]]+evalsCB[[i]],{i,1,13}],Joined->True,Filling->Table[n->evalsCB[[n]],{n,1,13}],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]],PlotRange->All,ImageSize->600,Frame->True,Axes->False];
Show[stateplotCB,potplottilt,PlotRange->{Automatic,{-4,4}}]

Time evolution of a groundstate :

ScaleFactor=10;
gs1=6;
stateplotCB=ListPlot[ScaleFactor evectsCB[[gs1]]+evalsCB[[gs1]],Joined->True,FillingStyle->Directive[Red,Opacity[0.3]],PlotRange->All,ImageSize->300,Frame->True,Axes->False];
Show[stateplotCB,potplottilt,PlotRange->{{0,500},{-4,4}}]

tfinal=Abs[2π/stateenergiesCB[[6]]]

MProblem=DiffMatrix[mass,2 5./division,VCBpotential[0.1]];
uinitial=evectsCB[[gs1]];
answerCB[t_]=u[t]/.First@@NDSolve[{-I MProblem.u[t]==u'[t],u[0]==uinitial},u[t],{t,0,tfinal}];

Animate[Show[potplottilt,ListPlot[ScaleFactor Abs[answerCB[t]]+evalsCB[[gs1]],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]]],PlotRange->{{0,500},{-5,5}}],{t,0,tfinal},AnimationRunning->False]

Time evolution of an excited state :

ScaleFactor=4;
es1=13;
stateplotCB=ListPlot[ScaleFactor evectsCB[[es1]]+evalsCB[[es1]],Joined->True,Filling->evalsCB[[es1]],FillingStyle->Directive[Red,Opacity[0.3]],PlotRange->All,ImageSize->500,Frame->True,Axes->False];
Show[stateplotCB,potplottilt,PlotRange->{{0,500},{-4,4}}]

MProblem=DiffMatrix[mass,2 5./division,VCBpotential[0.1]];
uinitial=evectsCB[[es1]];
answerCBexc[t_]=u[t]/.First@@NDSolve[{-I MProblem.u[t]==u'[t],u[0]==uinitial},u[t],{t,0,tfinal},StartingStepSize->0.3];

Animate[Show[potplottilt,ListPlot[ScaleFactor Abs[answerCBexc[t]]+evalsCB[[es1]],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]]],PlotRange->{{0,500},{-2,2}}],{t,0,tfinal},AnimationRunning->False]

And now we tilt the potential and look at the time evolution of the former excited eigenstate:

potplottilt2=ListPlot[VCBpotential[4],Filling->1.1Min[VCBpotential[4]],Joined->True,Frame->True,PlotRange->{{0,500},{-25,10}},GridLinesStyle->Directive[Thickness[0.002],Red],Axes->False]

MProblemCBtilt=DiffMatrix[mass,2 5./division,VCBpotential[4]];
answersCBtilt=Eigensystem[MProblemCBtilt];
{evalsCBtilt,evectsCBtilt}=Sort[Transpose[answersCBtilt],#1[[1]]<#2[[1]]&]//Transpose;
stateenergiesCBtilt=Take[evalsCBtilt,40]

ScaleFactor=3;
stateplotCBtilt=ListPlot[Table[ScaleFactor evectsCBtilt[[i]]+evalsCBtilt[[i]],{i,35,40}],Joined->True,Filling->Table[n->evalsCBtilt[[n+34]],{n,1,5}],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]],PlotRange->All,ImageSize->600,Frame->True,Axes->False];
Show[stateplotCBtilt,potplottilt2,PlotRange->{All,{-15,0}}]

ScaleFactor=4;
gs1tilt=36;
es1tilt=39;
stateplotCBtilt=ListPlot[{ScaleFactor evectsCBtilt[[es1tilt]]+evalsCBtilt[[es1tilt]],ScaleFactor evectsCBtilt[[gs1tilt]]+evalsCBtilt[[gs1tilt]]},Joined->True,Filling->{1->evalsCBtilt[[es1tilt]],2->evalsCBtilt[[gs1tilt]]},FillingStyle->Directive[Red,Opacity[0.3]],PlotRange->All,ImageSize->500,Frame->True,Axes->False];
Show[stateplotCBtilt,potplottilt2,PlotRange->{{0,500},{-10,-5}}]

MProblem=DiffMatrix[mass,2 5./division,VCBpotential[4]];
uinitial=evectsCB[[es1]];
answerCBexctilt[t_]=u[t]/.First@@NDSolve[{-I MProblem.u[t]==u'[t],u[0]==uinitial},u[t],{t,0,5.}];

Animate[Show[potplottilt2,ListPlot[20 Abs[answerCBexctilt[t]],PlotStyle->Red],PlotRange->{{0,500},{-10,5}}],{t,0,5.},AnimationRunning->False]

What happens to the former groundstate?

MProblem=DiffMatrix[1,2 5./division,VCBpotential[4]];
uinitial=evectsCB[[gs1]];
answerCBgstilt[t_]=u[t]/.First@@NDSolve[{-I MProblem.u[t]==u'[t],u[0]==uinitial},u[t],{t,0,tfinal},StartingStepSize->1];

Animate[Show[potplottilt2,ListPlot[20 Abs[answerCBgstilt[t]],PlotStyle->Red,FillingStyle->Directive[Red,Opacity[0.3]]],PlotRange->{{0,500},{-10,10}}],{t,0,2tfinal},AnimationRunning->False]

© Stefan Filipp, 10/2008