3
$\begingroup$

We know that a unitary transformation should not change any physical information. In the program I’m working on, my model is like the driven Jaynes-Cummings model. I want to compute the dynamics governed by the two Hamiltonians before and after the unitary transformation. When the parameters involved are relatively small, the fidelity between the two dynamics is exactly 1, which is what I expected.

The relative Hamiltonian is: $$ H_1 =\ \omega_M m^\dagger m + \frac{\omega_Q}{2}\sigma_z + G (\sigma_+ m + \sigma_- m^\dagger) + \frac{\eta_1}2 (\sigma_+ e^{-i \omega_1 t} + \sigma_- e^{i \omega_1 t}) + \frac{\eta_2}2 (\sigma_+ e^{-i \omega_2 t} + \sigma_- e^{i \omega_2 t}) + \frac{\eta_3}2 (m^\dagger e^{-i \omega_3 t} + m e^{i \omega_3 t}) $$ where $\eta_i$ is the amplitude, and $\omega_i$ is the frequency of the $i$-th drive. Moving into the rotating frame of the drive frequency $\omega_1$, we perform a unitary transformation $U_1 = \exp[i\omega_1(m^\dagger m + \sigma_z/2)]$, which leads to the following Hamiltonian: $$ H_2 =\ U_1 H_1 U_1^\dagger - i U_1 \dot{U}_1^\dagger =\ \Delta_m m^\dagger m + \frac{\Delta_q}{2}\sigma_z + G (\sigma_+ m + \sigma_- m^\dagger) + \frac{\eta_1}{2} (\sigma_+ + \sigma_-) + \frac{\eta_2}{2} (\sigma_+ e^{i(\omega_1 - \omega_2) t} + \sigma_- e^{-i(\omega_1 - \omega_2) t}) + \frac{\eta_3}{2} (m^\dagger e^{i(\omega_1 - \omega_3) t} + m e^{-i(\omega_1 - \omega_3) t})+ \frac{\eta_1}{2} (\sigma_+ e^{2i\omega_1 t} + \sigma_- e^{-2i\omega_1 t}) + \frac{\eta_2}{2} (\sigma_+ e^{i(\omega_1 + \omega_2) t} + \sigma_- e^{-i(\omega_1 + \omega_2) t}) + \frac{\eta_3}{2} (m^\dagger e^{i(\omega_1 + \omega_3) t} + m e^{-i(\omega_1 + \omega_3) t}) $$ where $\Delta_m = \omega_M - \omega_1$ and $\Delta_q = \omega_Q - \omega_1$.

The formula $F=Tr[\rho_1.\rho_2]$ is used to evaluate the fidelity between two density matrices. The code is as follows:

n =10;\[Omega]m=0.8;\[Omega]q=0.9;G=0.01;\[Eta]1=0.2;\[Eta]2=0.3;\[Eta]3=0.4;\[CapitalOmega]=0.5; \[Omega]1=0.1;\[Omega]2=0.15;\[Omega]3=0.18;\[Omega]p=0.21; aIdentity = IdentityMatrix[2]; bIdentity= IdentityMatrix[n]; annihilation= Normal[SparseArray[(Flatten[#1, 1] & )[Table[{{i, i + 1} -> Sqrt[i], {i + 1, i} -> 0}, {i, n - 1}]]]]; adag = Normal[SparseArray[(Flatten[#1, 1] & )[Table[{{i, i + 1} -> 0, {i + 1, i} -> Sqrt[i]}, {i, n - 1}]]]]; adagTa = Normal[SparseArray[Table[{i + 1, i + 1} -> i, {i, n - 1}]]]; {sx, sy, sz} = Table[PauliMatrix[i], {i, 3}]; sp = (sx + I*sy)/2.; sm = (sx - I*sy)/2.; hamiltonian[t_]:=Module[{h1, h2,hi ,hd1 ,hd2 ,hd3 ,hd}, h1 = \[Omega]m KroneckerProduct[adagTa, aIdentity] + \[Omega]q/2 KroneckerProduct[bIdentity, sz] + G (KroneckerProduct[annihilation, sp] + KroneckerProduct[adag, sm]); hd1=\[Eta]1/2. KroneckerProduct[bIdentity, sp E^(-I \[Omega]1 t)+sm E^(I \[Omega]1 t)]; hd2=\[Eta]2/2. KroneckerProduct[bIdentity, sp E^(-I \[Omega]2 t)+sm E^(I \[Omega]2 t)]; hd3=\[Eta]3/2. KroneckerProduct[adag E^(-I \[Omega]3 t) + annihilation E^(I \[Omega]3 t), aIdentity]; hd = \[CapitalOmega] KroneckerProduct[bIdentity, sp E^(-I \[Omega]p t) + sm E^(I \[Omega]p t)]; Chop[h1 + hd1 + hd2 + hd3 + had] ] functionlist=Table[\[CapitalPsi][k][t],{k,1,2n}]; namelist=Table[\[CapitalPsi][k],{k,1,2n}]; eqns=Thread[I*D[functionlist,t]==hamiltonian[t].functionlist]; initial=N[Flatten[KroneckerProduct[ReplacePart[Table[0, n], 1 -> 1], {0, 1}]]]; \[CapitalPsi]0=Thread[functionlist==initial]/.t->0; Clear[initial]; tf= 45./(2. \[Pi]*16.*10^6); sol=NDSolveValue[Join[eqns,\[CapitalPsi]0],namelist,{t,0,tf},Method->"StiffnessSwitching",AccuracyGoal->15,PrecisionGoal->15,MaxSteps->Infinity]; wavfn3=Table[Table[sol[[i]][t],{i,1,2n}]//Chop,{t,0,tf,0.005tf}]; hamiltonian[t_]:=Module[{h1,hd1 ,hd2 ,hd3 ,h5 ,h6 ,hd}, h1 = (\[Omega]m-\[Omega]1) KroneckerProduct[adagTa, aIdentity]+ (\[Omega]q-\[Omega]1)/2 KroneckerProduct[bIdentity, sz]+ G (KroneckerProduct[annihilation, sp] + KroneckerProduct[adag, sm]); hd1=\[Eta]1/2. KroneckerProduct[bIdentity, sp E^(I (\[Omega]1-\[Omega]1) t) + sm E^(-I (\[Omega]1-\[Omega]1) t)]; hd2=\[Eta]2/2. KroneckerProduct[bIdentity, sp E^(I (\[Omega]1-\[Omega]2) t) + sm E^(-I (\[Omega]1-\[Omega]2) t)]; hd3=\[Eta]3/2. KroneckerProduct[adag E^(I (\[Omega]1-\[Omega]3) t)+annihilation E^(-I (\[Omega]1-\[Omega]3) t), aIdentity]; hd = \[CapitalOmega] KroneckerProduct[bIdentity, sp E^(I (\[Omega]1-\[Omega]p) t) + sm E^(-I (\[Omega]1-\[Omega]p) t)]; Chop[h1 + hd1 +hd2 + hd3 + had] ] functionlist=Table[\[CapitalPsi][k][t],{k,1,2n}]; namelist=Table[\[CapitalPsi][k],{k,1,2n}]; eqns=Thread[I*D[functionlist,t]==hamiltonian[t].functionlist]; initial=N[Flatten[KroneckerProduct[ReplacePart[Table[0, n], 1 -> 1], {0, 1}]]]; \[CapitalPsi]0=Thread[functionlist==initial]/.t->0; Clear[initial]; tf= 45./(2. \[Pi]*16.*10^6); sol=NDSolveValue[Join[eqns,\[CapitalPsi]0],namelist,{t,0,tf},Method->"StiffnessSwitching",AccuracyGoal->15,PrecisionGoal->15,MaxSteps->Infinity]; wavfn4=Table[Table[sol[[i]][t],{i,1,2n}]//Chop,{t,0,tf,0.005tf}]; fun[vector_]:=Module[{vec,rho}, vec=Normalize[vector]; rho = ConjugateTranspose[{vec}].{vec}; rho ] data34=Re@Table[{0.01 i,N[Tr[fun[wavfn3[[i]]].fun[wavfn4[[i]]]],5]},{i,1,Length[wavfn4],1}]; plt34=ListLinePlot[data34,Frame->True,PlotStyle->{Darker@Red,Thickness[0.0035]}, ImageSize->350,PlotRange->{All,{0,1.05}}, FrameTicksStyle->Directive[Black,12],FrameStyle->Directive[Thickness[0.0035],Black,FontSize->20], FrameLabel->(Style[#,FontFamily->"Times",FontColor->Black,FontSize->15]&/@{"time","Fidelity"}) ]

`

and the result is the fidelity

However, when the parameters become large (considering actual system parameters, such as (G) Hz) such as

G= 2\[Pi]*10.*10^6; \[Omega]m=2\[Pi]*8.16764*10^9; \[Omega]q=2\[Pi]*8.22352*10^9; \[Eta]2=2\[Pi]*0.05*10^9; \[Eta]3=2\[Pi]*0.05*10^9; \[Omega]1=2\[Pi]*8.16*10^9; \[Omega]2=2\[Pi]*7.66*10^9; \[Eta]1=\[Omega]1-\[Omega]2; \[Omega]3=\[Omega]1; \[Omega]p=Sqrt[2]/2. \[Eta]2+\[Omega]1; \[CapitalOmega]=Sqrt[2]/4. G; \[Nu]=Sqrt[2]/2. \[Eta]2; gx=G/2. Sqrt[2.]/2.; geff = (4.gx^2)/\[Nu]; \[CapitalDelta]=8. gx^2/(3\[Nu]); \[CapitalDelta]m=\[Omega]m-\[Omega]1; \[CapitalDelta]p=\[Omega]p-\[Omega]1; \[Delta]z=\[Eta]2/2.; \[Delta]x=(\[Eta]3 G)/\[CapitalDelta]m;

I noticed that the computed fidelity is quite poor. Like this

enter image description here

I understand that this might be due to computational errors arising from the numerical solution of the equations. In my code, I am using NDSolveValue to solve the Schrödinger equation and have tried various function options to improve accuracy, but none of these changes have resolved the issue and some even led to errors in the program. I was wondering if anyone could offer insights into where the problem might lie? Any discussion would be greatly appreciated. Thank you so much for your time!"

Improve this question
$\endgroup$
4
  • $\begingroup$ It's not clear what you're trying to figure out? It looks like you used the same names for different solutions. $\endgroup$ Commented Jan 11 at 14:12
  • $\begingroup$ For these two different Hamiltonians (related by a unitary transformation), I want to evaluate the fidelity. However, for actual system parameters, the fidelity is not equal to 1. Mathematically, the fidelity should strictly be 1, which indeed holds for smaller parameters. (I use the function names wavfn3 and wavfn4 to differentiate the computed results.). I don't know why this is happening. I suspect it might be due to insufficient computational precision causing errors in the results, but I haven't been able to resolve the issue. $\endgroup$ Commented Jan 11 at 17:08
  • $\begingroup$ There are several typos in your code. It could be better to add some equations in Latex form to your post. $\endgroup$ Commented Jan 11 at 23:46
  • $\begingroup$ Thanks. The Hamiltonian formula in the LaTeX form is added to my question. Thank you very much again. $\endgroup$ Commented Jan 12 at 0:54

1 Answer 1

Reset to default
1
$\begingroup$

Using your code as basis, I removed the "had" and merged the two hamiltonian definitions, and the output seem okay for large parameters which are different from yours:

n = 30; G = 2 \[Pi]*10.*10^6; \[Omega]m = 2 \[Pi]*8.16764*10^7; \[Omega]q = 2 \[Pi]*8.22352*10^7; \[Eta]2 = 2 \[Pi]*0.05*10^7; \[Eta]3 = 2 \[Pi]*0.05*10^7; \[Omega]1 = 2 \[Pi]*8.16*10^7; \[Omega]2 = 2 \[Pi]*7.66*10^7; \[Eta]1 = \[Omega]1 - \[Omega]2; \[Omega]3 = \[Omega]1; \[Omega]p = Sqrt[2]/2. \[Eta]2 + \[Omega]1; \[CapitalOmega] = Sqrt[2]/4. G; \[Nu] = Sqrt[2]/2. \[Eta]2; gx = G/2. Sqrt[2.]/2.; geff = (4. gx^2)/\[Nu]; \[CapitalDelta] = 8. gx^2/(3 \[Nu]); \[CapitalDelta]m = \[Omega]m - \[Omega]1; \[CapitalDelta]p = \[Omega]p - \[Omega]1; \[Delta]z = \[Eta]2/2.; \[Delta]x = (\[Eta]3 G)/\[CapitalDelta]m; aIdentity = IdentityMatrix[2]; bIdentity = IdentityMatrix[n]; annihilation = Normal[SparseArray[ Flatten[Table[{{i, i + 1} -> Sqrt[i], {i + 1, i} -> 0}, {i, n - 1}], 1]]]; adag = Normal[ SparseArray[ Flatten[Table[{{i, i + 1} -> 0, {i + 1, i} -> Sqrt[i]}, {i, n - 1}], 1]]]; adagTa = Normal[SparseArray[Table[{i + 1, i + 1} -> i, {i, n - 1}]]]; {sx, sy, sz} = Table[PauliMatrix[i], {i, 3}]; sp = (sx + I sy)/2.; sm = (sx - I sy)/2.; hamiltonian[t_] := Module[{h1, hd1, hd2, hd3, hd}, h1 = \[Omega]m KroneckerProduct[adagTa, aIdentity] + \[Omega]q/ 2 KroneckerProduct[bIdentity, sz] + G (KroneckerProduct[annihilation, sp] + KroneckerProduct[adag, sm]); hd1 = \[Eta]1/2 KroneckerProduct[bIdentity, sp Exp[-I \[Omega]1 t] + sm Exp[I \[Omega]1 t]]; hd2 = \[Eta]2/2 KroneckerProduct[bIdentity, sp Exp[-I \[Omega]2 t] + sm Exp[I \[Omega]2 t]]; hd3 = \[Eta]3/2 KroneckerProduct[ adag Exp[-I \[Omega]3 t] + annihilation Exp[I \[Omega]3 t], aIdentity]; hd = \[CapitalOmega] KroneckerProduct[bIdentity, sp Exp[-I \[Omega]p t] + sm Exp[I \[Omega]p t]]; Chop[h1 + hd1 + hd2 + hd3 + hd]]; functionlist = Table[\[CapitalPsi][k][t], {k, 1, 2 n}]; namelist = Table[\[CapitalPsi][k], {k, 1, 2 n}]; eqns = Thread[I*D[functionlist, t] == hamiltonian[t] . functionlist]; initial = N[Flatten[ KroneckerProduct[ReplacePart[Table[0, n], 1 -> 1], {0, 1}]]]; \[CapitalPsi]0 = Thread[functionlist == initial] /. t -> 0; Clear[initial]; tf = 45./(2. \[Pi]*16.*10^6); sol = NDSolveValue[Join[eqns, \[CapitalPsi]0], namelist, {t, 0, tf}, Method -> "StiffnessSwitching", AccuracyGoal -> 15, PrecisionGoal -> 15, MaxSteps -> Infinity]; wavfn3 = Table[Table[sol[[i]][t], {i, 1, 2 n}] // Chop, {t, 0, tf, 0.005 tf}]; fun[vector_] := Module[{vec, rho}, vec = Normalize[vector]; rho = ConjugateTranspose[{vec}] . {vec}; rho]; data34 = Re@Table[{0.01 i, N[Tr[fun[wavfn3[[i]]] . fun[wavfn3[[i]]]], 5]}, {i, 1, Length[wavfn3], 1}]; plt34 = ListLinePlot[data34, Frame -> True, PlotStyle -> {Darker@Red, Thickness[0.0035]}, ImageSize -> 350, PlotRange -> {All, {0, 1.05}}, FrameTicksStyle -> Directive[Black, 12], FrameStyle -> Directive[Thickness[0.0035], Black, FontSize -> 20], FrameLabel -> (Style[#, FontFamily -> "Times", FontColor -> Black, FontSize -> 15] & /@ {"time", "Fidelity"})]

Fidelity stays constant

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ Thank you for helping me with my question about the fidelity check. In your code as if only contain a Hamiltonian and evaluate the fidelity themselves. In my question, I need to evaluate the dynamic fidelity under the different Hamiltonians, and the two Hamiltonians are connected by a unitary transformation. $\endgroup$ Commented Jan 12 at 6:19
  • $\begingroup$ Ok I understand your original question more now, your second Hamiltonian is a duplicate of the first, perhaps you can type in the second Hamiltonian as you intended as given in the second equation. $\endgroup$ Commented Jan 12 at 6:40
  • $\begingroup$ i am trying again $\endgroup$ Commented Jan 12 at 9:52

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.