Why does NDSolve fail or slow down when simulating a driven JC model with time-dependent terms, even when the drive is off?

Question

I'm performing a numerical check of an effective Hamiltonian transformation applied to a driven Jaynes–Cummings (JC) model. After applying a rotating frame transformation and rotating-wave approximation (RWA), the system reduces to an effective Rabi model. The procedure is described in detail in this paper:[https://arxiv.org/pdf/1107.5748]

To validate this approximation, I simulate dynamics under both the original driven JC model and the effective Rabi model, and compare an observable (a squeezing-related quantity) over time.

Working example: effective Rabi model

This simulation runs quickly and stably, even with physical GHz frequencies:

\[Omega]m=2\[Pi]*8.16734*10^9;   
\[Omega]q=2\[Pi]*8.22352*10^9;   
G =2.\[Pi]*10.*10^6; 

\[CapitalOmega]1=2\[Pi]*1000.*10^6; 
\[CapitalOmega]2=2\[Pi]*125.0005*10^6;  
\[Omega]1=2\[Pi]*8.16654*10^9; 
\[Omega]2=2\[Pi]*6.16654*10^9; 

\[Omega]=\[Omega]m-\[Omega]1;
\[CapitalOmega]=\[CapitalOmega]2;
g = G/2.; 

n=45;
hRabi= \[Omega] KroneckerProduct[adagTa[n], aIdentity] + \[CapitalOmega]/2     KroneckerProduct[bIdentity[n], sz]+ g KroneckerProduct[aPadag[n], sx];

functionlist=Table[\[CapitalPsi][k][t],{k,1,2n}];
namelist=Table[\[CapitalPsi][k],{k,1,2n}];
eqns=Thread[I*D[functionlist,t]==hRabi . functionlist];
initial=N[Flatten[KroneckerProduct[ReplacePart[Table[0, n], 1 -> 1], {0, 1}]]];
\[CapitalPsi]0=Thread[functionlist==initial]/.t->0;

tfine= 5.*10^-6;
sol=NDSolveValue[Join[eqns,\[CapitalPsi]0],namelist,{t,0,tfine},AccuracyGoal->8,PrecisionGoal->8,MaxSteps->Infinity];
\[Rho]t[t_]:=rdm[(KroneckerProduct[#, #\[Conjugate]] &@Normalize[Chop[Table[sol[[i]][t],{i,1,2n}]]]),{n,2,1}];

x\[Theta]=1/Sqrt[2.] (annihilation[n]E^(-I \[Theta])+adag[n]E^(I \[Theta]));
var[t_]:=NMinimize[Chop@Re[Tr[\[Rho]t[t] . (x\[Theta] . x\[Theta])]-Tr[\[Rho]t[t] . (x\[Theta])]^2],{\[Theta],0,2Pi}][[1]];
\[Xi][t_]:=-10. Log[10,2. var[t]];

data1=Table[{t,var[t]},{t,0.,tfine,tfine/100}];
data2=Table[{t,\[Xi][t]},{t,0.,tfine,tfine/100}];

ListLinePlot/@{data1,data2}

the output result is

Problematic version: full time-dependent driven JC model

The only change is replacing hRabi with a time-dependent Hamiltonian with drives. Even when drive amplitudes are set to 0, this code either:

fails with NDSolveValue::ndnum, or

runs very slowly, even at t = 0

Here is the part code:

hjc[t_?NumericQ] := Module[{h1, hd},
   h1 = \[Omega]m KroneckerProduct[adagTa[n], aIdentity] + \[Omega]q/2 KroneckerProduct[bIdentity[n], sz] + g (KroneckerProduct[annihilation[n], sp] + KroneckerProduct[adag[n], sm]);
   hd = Evaluate[\[CapitalOmega]1 KroneckerProduct[bIdentity[n], sp Exp[-I \[Omega]1 t] + sm Exp[I \[Omega]1 t]]+ \[CapitalOmega]2 KroneckerProduct[bIdentity[n], sp Exp[-I \[Omega]2 t] + sm Exp[I \[Omega]2 t]]];
   Chop[h1 + hd]
   ]

functionlist=Table[\[CapitalPsi][k][t],{k,1,2n}];
namelist=Table[\[CapitalPsi][k],{k,1,2n}];
eqns=Thread[I*D[functionlist,t]==hjc[t] . functionlist];
initial=N[Flatten[KroneckerProduct[ReplacePart[Table[0, n], 1 -> 1], {0, 1}]]];
\[CapitalPsi]0=Thread[functionlist==initial]/.t->0;

My questions:

Why is this failing or extremely slow, even when drive terms are zero?

Are there standard tricks (e.g. frequency rescaling, frame transformations, solver options) for high-frequency, time-dependent quantum systems in NDSolve?

Is NDSolveValue fundamentally unsuitable for such stiff time-dependent systems at GHz?

Any suggestions or experience with this kind of problem would be greatly appreciated. Thank you!

Some basic functions defined by myself:

annihilation::usage = "annihilation operator"; 
adag::usage = "creation operator"; 
aPadag::usage = "creation operator plus annihilation operator    a+adag"; 
adagTa::usage = "creation operator times annihilation (number) operator adag.a"; 
sx::usage = "pauli operator sx"; 
sy::usage = "pauli operator sy"; 
sz::usage = "pauli operator sz"; 
sp::usage = "up ladder operator sp"; 
sm::usage = "down ladder operator sm"; 
aIdentity::usage = "two level identity matrix"; 
bIdentity::usage = "bosic field identity matrix"; 
rdm::usage = "\[Rho]A=rdm[\[Rho]AB,{dA,dB,1}]; \[Rho]B=rdm[\[Rho]AB,{1,dA,dB}];eg,\[Rho]=KroneckerProduct[aTa[n],sp],\[Rho]A=rdm[\[Rho],{n,2,1}];"; 

rdm[(\[Rho]ABC_)?MatrixQ, {dA_Integer /; dA >= 1, dB_Integer /; dB >= 1, dC_Integer /; dC >= 1}] /; Dimensions[\[Rho]ABC] == {dA*dB*dC, dA*dB*dC} := Flatten[TensorContract[ArrayReshape[\[Rho]ABC, {dA, dB, dC, dA, dB, dC}], {2, 5}], {{1, 2}, {3, 4}}]; 

annihilation[n_] := SparseArray[(Flatten[#1, 1] & )[Table[{{i, i + 1} -> Sqrt[i], {i + 1, i} -> 0}, {i, n - 1}]]]; 
adag[n_] := SparseArray[(Flatten[#1, 1] & )[Table[{{i, i + 1} -> 0, {i + 1, i} -> Sqrt[i]}, {i, n - 1}]]]; 
adagTa[n_] := SparseArray[Table[{i + 1, i + 1} -> i, {i, n - 1}]]; 
aPadag[n_] := SparseArray[(Flatten[#1, 1] & )[Table[{{i, i + 1} -> Sqrt[i], {i + 1, i} -> Sqrt[i]}, {i, n - 1}]]]; 
{sx, sy, sz} = Table[SparseArray[PauliMatrix[i]], {i, 3}]; 
sp = (sx + I*sy)/2; 
sm = (sx - I*sy)/2; 
aIdentity = IdentityMatrix[2, SparseArray]; 
bIdentity[n_] := IdentityMatrix[n, SparseArray];

When environmental effects are neglected, I can obtain the system’s evolution by directly solving the Schrödinger equation. However, in more general situations where dissipation is considered, I need to solve the Lindblad master equation to extract the system’s dynamics.

So far, I’ve successfully implemented the Lindblad evolution for the Rabi model. Since the Hamiltonian is time-independent, the computation runs very efficiently. However, when I switch to a time-dependent Jaynes-Cummings (JC) model, the simulation becomes significantly slower, even with a small Hilbert space truncation (e.g., n = 10). This slowdown is especially evident when solving the master equation using NDSolve.

Master equation code to solve the Rabi model as follow:

rdm[(\[Rho]ABC_)?MatrixQ, {dA_Integer /; dA >= 1,      dB_Integer /; dB >= 1, dC_Integer /; dC >= 1}] /;    Dimensions[\[Rho]ABC] == {dA*dB*dC, dA*dB*dC} :=   Flatten[TensorContract[    ArrayReshape[\[Rho]ABC, {dA, dB, dC, dA, dB, dC}], {2, 5}], {{1,      2}, {3, 4}}];

annihilation[n_] :=   SparseArray[(Flatten[#1, 1] &)[
Table[{{i, i + 1} -> Sqrt[i], {i + 1, i} -> 0}, {i, n - 1}]]];
adag[n_] :=   SparseArray[(Flatten[#1, 1] &)[
Table[{{i, i + 1} -> 0, {i + 1, i} -> Sqrt[i]}, {i, n - 1}]]];
adagTa[n_] := SparseArray[Table[{i + 1, i + 1} -> i, {i, n - 1}]];
aPadag[n_] :=   SparseArray[(Flatten[#1, 1] &)[
Table[{{i, i + 1} -> Sqrt[i], {i + 1, i} -> Sqrt[i]}, {i, 
  n - 1}]]];
{sx, sy, sz} = Table[SparseArray[PauliMatrix[i]], {i, 3}];
sp = (sx + I*sy)/2;
sm = (sx - I*sy)/2;
aIdentity = IdentityMatrix[2, SparseArray];
bIdentity[n_] := IdentityMatrix[n, SparseArray];
\[Omega]m = 2 \[Pi]*8.01*10^9;
\[Omega]q = 2 \[Pi]*8.*10^9;
G = 2. \[Pi]*10.*10^6;

\[CapitalOmega]1 = 2 \[Pi]*.7*10^9;
\[CapitalOmega]2 = 2 \[Pi]*10^7;
\[Omega]1 = 2 \[Pi]*8.*10^9;
\[Omega]2 = 2 \[Pi]*6.6*10^9;

\[Omega] = \[Omega]m - \[Omega]1;
\[CapitalOmega] = \[CapitalOmega]2;
g = \[Omega];
n = 45;
hRabi = KroneckerProduct[adagTa[n], 
aIdentity] + \[CapitalOmega]/\[Omega]/2 KroneckerProduct[
 bIdentity[n], sz] - g /\[Omega] KroneckerProduct[aPadag[n], sx];

lindblad[o1_,o2_, rho_] :=  Module[{op, oo},        op = SparseArray[KroneckerProduct[o1, o2]];         2. op . rho . op\[ConjugateTranspose] - rho . op\[ConjugateTranspose] . op - op\[ConjugateTranspose] . op . rho     ];
\[Rho] = Array[Subscript[rho, #1, #2][t] &, {2 n, 2 n}];
L\[Sigma]m = N@lindblad[bIdentity[n],sm, \[Rho]];
L\[Sigma]z = N@lindblad[bIdentity[n],sz, \[Rho]];
La = N@lindblad[annihilation[n],aIdentity, \[Rho]];

\[Gamma] = 2. \[Pi]*0.*10^6;
\[Kappa] = 2. \[Pi]*0.*10^6;

tfinal = 4. Pi;

L\[Rho]=N[-I (hRabi . \[Rho] - \[Rho] . hRabi) + (\[Gamma] /\[Omega])/2. L\[Sigma]m + (\[Gamma] /\[Omega])/2. L\[Sigma]z + (\[Kappa] /\[Omega])/2. La];
initial =  KroneckerProduct[#\[Conjugate], #] &@N[Flatten[KroneckerProduct[ReplacePart[Table[0, n], 1 -> 1], {0, 1}]]];

sol=Flatten@NDSolve[L\[Rho]==\!\(\*SubscriptBox[\(\[PartialD]\), \(t\)]\[Rho]\) && (\[Rho]/.t->0)==initial,Flatten@\[Rho],{t,0,tfinal}];
\[Rho]t=Flatten[\[Rho]/.sol,1];
\[Rho]a[t_]:=Evaluate[rdm[Partition[\[Rho]t,2n],{n,2,1}]];
x\[Theta] =   1/Sqrt[2.] (annihilation[n] E^(-I \[Theta]) +      adag[n] E^(I \[Theta]));
var[t_] :=   NMinimize[
Chop@Re[Tr[\[Rho]a[t] . (x\[Theta] . x\[Theta])] - 
   Tr[\[Rho]a[t] . (x\[Theta])]^2], {\[Theta], 0, 2 Pi}][[1]];
\[Xi][t_] := -10. Log[10, 2. var[t]];

data11 = Table[{t, Chop@var[t]}, {t, 0., tfinal, tfinal/100}];
data12 = Table[{t, Chop@\[Xi][t]}, {t, 0., tfinal, tfinal/100}];

ListPlot[{data11, data12}, PlotStyle -> {Red, Blue},  PlotLegends -> Automatic]

When considering dissipation is zero, the result is the same as solving the Schrödinger equation. The output figure is

But when rabi Hamiltonian is replaced by jc Hamiltonian :

hjc=    \[Omega]m/\[Omega] KroneckerProduct[adagTa[n], aIdentity] + 
\[Omega]q/\[Omega]/2 KroneckerProduct[bIdentity[n], sz]
-    2 g/\[Omega] (KroneckerProduct[annihilation[n], sp] + 
 KroneckerProduct[adag[n], sm]) 
-   \[CapitalOmega]1/\[Omega] KroneckerProduct[
 bIdentity[n], sp Exp[-I \[Omega]1/\[Omega] t] + 
 sm Exp[I \[Omega]1/\[Omega] t]]
-    \[CapitalOmega]2/\[Omega] KroneckerProduct[
 bIdentity[n], sp Exp[-I \[Omega]2/\[Omega] t] + 
 sm Exp[I \[Omega]2/\[Omega] t]]//Re//Chop;

The code run is very slowly, and the output result is different

Apologies for the oversight — I had forgotten to include some of my custom-defined functions in the original post. I've updated the post and added them now. Sorry again for the inconvenience, and thank you for your patience. — gang liu
– gang liu, Commented Jul 10 at 1:17
Thank you for your update (+1). Now first code is working, but second code looks like incomplete. — Alex Trounev
– Alex Trounev, Commented Jul 10 at 13:29
Thank you for helping test my code. The only difference between the first and second code lies in the Hamiltonian. I made a mistake in the second code: in the eqns, the Hamiltonian should be hjc[t] instead of hjc1. — gang liu
– gang liu, Commented Jul 11 at 3:32
I have tested second code. It computes same results as the first code, but 5 times faster. — Alex Trounev
– Alex Trounev, Commented Jul 12 at 1:42

Alex Trounev · Accepted Answer · 2025-07-13 07:29:49Z

We used same parameters and scaling as in the paper Quantum Simulation of the Ultrastrong Coupling Dynamics in Circuit QED. First code with Hamiltonian from the paper

annihilation::usage = "annihilation operator";
adag::usage = "creation operator";
aPadag::usage = 
  "creation operator plus annihilation operator    a+adag";
adagTa::usage = 
  "creation operator times annihilation (number) operator adag.a";
sx::usage = "pauli operator sx";
sy::usage = "pauli operator sy";
sz::usage = "pauli operator sz";
sp::usage = "up ladder operator sp";
sm::usage = "down ladder operator sm";
aIdentity::usage = "two level identity matrix";
bIdentity::usage = "bosic field identity matrix";
rdm::usage = 
  "\[Rho]A=rdm[\[Rho]AB,{dA,dB,1}]; \
\[Rho]B=rdm[\[Rho]AB,{1,dA,dB}];eg,\[Rho]=KroneckerProduct[aTa[n],sp],\
\[Rho]A=rdm[\[Rho],{n,2,1}];";

rdm[(\[Rho]ABC_)?MatrixQ, {dA_Integer /; dA >= 1, 
     dB_Integer /; dB >= 1, dC_Integer /; dC >= 1}] /; 
   Dimensions[\[Rho]ABC] == {dA*dB*dC, dA*dB*dC} := 
  Flatten[TensorContract[
    ArrayReshape[\[Rho]ABC, {dA, dB, dC, dA, dB, dC}], {2, 5}], {{1, 
     2}, {3, 4}}];

annihilation[n_] := 
  SparseArray[(Flatten[#1, 1] &)[
    Table[{{i, i + 1} -> Sqrt[i], {i + 1, i} -> 0}, {i, n - 1}]]];
adag[n_] := 
  SparseArray[(Flatten[#1, 1] &)[
    Table[{{i, i + 1} -> 0, {i + 1, i} -> Sqrt[i]}, {i, n - 1}]]];
adagTa[n_] := SparseArray[Table[{i + 1, i + 1} -> i, {i, n - 1}]];
aPadag[n_] := 
  SparseArray[(Flatten[#1, 1] &)[
    Table[{{i, i + 1} -> Sqrt[i], {i + 1, i} -> Sqrt[i]}, {i, 
      n - 1}]]];
{sx, sy, sz} = Table[SparseArray[PauliMatrix[i]], {i, 3}];
sp = (sx + I*sy)/2;
sm = (sx - I*sy)/2;
aIdentity = IdentityMatrix[2, SparseArray];
bIdentity[n_] := IdentityMatrix[n, SparseArray];
\[Omega]m = 2 \[Pi]*8.01*10^9;
\[Omega]q = 2 \[Pi]*8.*10^9;
G = 2. \[Pi]*20.*10^6;

\[CapitalOmega]1 = 2 \[Pi]*.7*10^9;
\[CapitalOmega]2 = 2 \[Pi]*10^7;
\[Omega]1 = 2 \[Pi]*8.*10^9;
\[Omega]2 = 2 \[Pi]*6.6*10^9;

\[Omega] = \[Omega]m - \[Omega]1;
\[CapitalOmega] = \[CapitalOmega]2;
g = \[Omega];
n = 45;
hRabi = KroneckerProduct[adagTa[n], 
    aIdentity] + \[CapitalOmega]/\[Omega]/2 KroneckerProduct[
     bIdentity[n], sz] - g /\[Omega] KroneckerProduct[aPadag[n], sx];

functionlist = Table[\[CapitalPsi][k][t], {k, 1, 2 n}];
namelist = Table[\[CapitalPsi][k], {k, 1, 2 n}];
eqns = Thread[I*D[functionlist, t] == hRabi.functionlist];
initial = 
  N[Flatten[
    KroneckerProduct[ReplacePart[Table[0, n], 1 -> 1], {0, 1}]]];
\[CapitalPsi]0 = Thread[functionlist == initial] /. t -> 0;

tfine = 4 Pi;
sol = NDSolveValue[Join[eqns, \[CapitalPsi]0], 
    namelist, {t, 0, tfine}]; // AbsoluteTiming
\[Rho]t[t_] := 
  rdm[(KroneckerProduct[#, #\[Conjugate]] &@
     Normalize[Chop[Table[sol[[i]][t], {i, 1, 2 n}]]]), {n, 2, 1}];
x\[Theta] = 
  1/Sqrt[2.] (annihilation[n] E^(-I \[Theta]) + 
     adag[n] E^(I \[Theta]));
var[t_] := 
  NMinimize[
    Chop@Re[Tr[\[Rho]t[t].(x\[Theta].x\[Theta])] - 
       Tr[\[Rho]t[t].(x\[Theta])]^2], {\[Theta], 0, 2 Pi}][[1]];
\[Xi][t_] := -10. Log[10, 2. var[t]];

data1 = Table[{t, var[t]}, {t, 0., tfine, tfine/100}];
data2 = Table[{t, \[Xi][t]}, {t, 0., tfine, tfine/100}];

Second code with Hamiltonian from the paper

hjc1 = \[Omega]m/\[Omega] KroneckerProduct[adagTa[n], 
     aIdentity] + \[Omega]q/\[Omega]/2 KroneckerProduct[bIdentity[n], 
     sz] - 2 g /\[Omega] (KroneckerProduct[annihilation[n], sp] + 
      KroneckerProduct[adag[n], 
       sm]) - \[CapitalOmega]1/\[Omega] KroneckerProduct[bIdentity[n],
      sp Exp[-I \[Omega]1/\[Omega] t] + 
      sm Exp[I \[Omega]1/\[Omega] t]] - \[CapitalOmega]2/\[Omega] \
KroneckerProduct[bIdentity[n], 
     sp Exp[-I \[Omega]2/\[Omega] t] + sm Exp[I \[Omega]2/\[Omega] t]];


functionlist = Table[\[CapitalPsi][k][t], {k, 1, 2 n}];
namelist = Table[\[CapitalPsi][k], {k, 1, 2 n}];
eqns1 = Thread[I*D[functionlist, t] == hjc1.functionlist];
initial = 
  N[Flatten[
    KroneckerProduct[ReplacePart[Table[0, n], 1 -> 1], {0, 1}]]];
\[CapitalPsi]0 = Thread[functionlist == initial] /. t -> 0;


sol1 = NDSolveValue[Join[eqns1, \[CapitalPsi]0], 
    namelist, {t, 0, tfine}]; // AbsoluteTiming

Off[General::munfl]

\[Rho]t1[t_] := 
  rdm[(KroneckerProduct[#, #\[Conjugate]] &@
     Normalize[Chop[Table[sol1[[i]][t], {i, 1, 2 n}]]]), {n, 2, 1}];

x\[Theta] = 
  1/Sqrt[2.] (annihilation[n] E^(-I \[Theta]) + 
     adag[n] E^(I \[Theta]));
var[t_] := 
  NMinimize[
    Chop@Re[Tr[\[Rho]t1[t].(x\[Theta].x\[Theta])] - 
       Tr[\[Rho]t1[t].(x\[Theta])]^2], {\[Theta], 0, 2 Pi}][[1]];
\[Xi][t_] := -10. Log[10, 2. var[t]];

data11 = Table[{t, Chop@var[t]}, {t, 0., tfine, tfine/100}];
data12 = Table[{t, Chop@\[Xi][t]}, {t, 0., tfine, tfine/100}];

Visualization

Show[ListLinePlot[{data1, data2}, PlotLegends -> Automatic, 
  PlotStyle -> {Blue, Red}, Frame -> True, 
  FrameLabel -> {"\[Omega]t", ""}], 
 ListPlot[{data11, data12}, PlotStyle -> {Red, Blue}, 
  PlotLegends -> Automatic]]

We see from these data that there is no complete match, but nevertheless the agreement is quite good.

Thank you so much for helping me with this question! Not only that, I also noticed that you kindly helped me solve another question I posted earlier. I truly appreciate the time you've taken to share your expertise and insight — it has been incredibly helpful to me. — gang liu
– gang liu, Commented Jul 14 at 3:59
When environmental effects are neglected, I can obtain the system’s evolution by directly solving the Schrödinger equation. However, in more general situations where dissipation is considered, I need to solve the Lindblad master equation to extract the system’s dynamics. — gang liu
– gang liu, Commented Jul 14 at 4:00
So far, I’ve successfully implemented the Lindblad evolution for the Rabi model. Since the Hamiltonian is time-independent, the computation runs very efficiently. However, when I switch to a time-dependent Jaynes-Cummings (JC) model, the simulation becomes significantly slower, even with a small Hilbert space truncation (e.g., n = 10). This slowdown is especially evident when solving the master equation using NDSolve. I’d like to ask: could this slowdown be due to the time dependence in the Hamiltonian, which might make the internal method used by NDSolve less efficient? — gang liu
– gang liu, Commented Jul 14 at 4:00
Is there any way to speed up the simulation for the time-dependent JC model? Any advice or suggestions would be greatly appreciated! Thank you again for your generous help! I’ve added the code I used for solving the master equation to the main post. — gang liu
– gang liu, Commented Jul 14 at 4:02
@gangliu You can check number of points used in simulation for tfine=1 as follows points = Length[#] & /@ {sol[[1]][[3]][[1]], sol1[[1]][[3]][[1]]}. On my old computer with Mathematica 11.3 it gives {39, 93095}. It is why computational times in this case are {0.0196606, 9.5627048}. We can't speed up due to large number of points used in sol1 which in turn depends on $\omega_1/\omega=800, \omega_2/\omega= 660$. — Alex Trounev
– Alex Trounev, Commented Jul 14 at 6:22

Stack Exchange Network

Why does NDSolve fail or slow down when simulating a driven JC model with time-dependent terms, even when the drive is off?

1 Answer 1

Your Answer

Hot Network Questions

Why does NDSolve fail or slow down when simulating a driven JC model with time-dependent terms, even when the drive is off?

1 Answer 1

Your Answer

Sign up or log in

Post as a guest

Related

Hot Network Questions