# Time dependent coupled first order differential density matrix equations

I want to solve the following time dependent coupled differential equations numerically. Please guide me.

{
{Derivative[1][Subscript[\[Rho], 1, 1]][t] == \[Gamma]/
2 - \[Gamma] Subscript[\[Rho], 1, 1][t] -
I (1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 3][t] -
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 1][t]) +
1/2 \[CapitalGamma] Subscript[\[Rho], 3, 3][t]},
{Derivative[1][Subscript[\[Rho], 1, 2]][
t] == -\[Gamma] Subscript[\[Rho], 1, 2][t] -
I (\[CapitalDelta]a Subscript[\[Rho], 1, 2][
t] - \[CapitalDelta]b Subscript[\[Rho], 1, 2][t] +
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 1, 3][t] -
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 2][t])},
{Derivative[1][Subscript[\[Rho], 1, 3]][t] ==
1/2 (-\[Gamma] Subscript[\[Rho], 1, 3][
t] - (\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 1, 3][
t]) - I (1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 1][t] +
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 1, 2][
t] + \[CapitalDelta]a Subscript[\[Rho], 1, 3][t] -
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 3][t])},
{Derivative[1][Subscript[\[Rho], 2, 1]][
t] == -\[Gamma] Subscript[\[Rho], 2, 1][t] -
I (-\[CapitalDelta]a Subscript[\[Rho], 2, 1][
t] + \[CapitalDelta]b Subscript[\[Rho], 2, 1][t] +
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 2, 3][t] -
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 1][t])},
{Derivative[1][Subscript[\[Rho], 2, 2]][t] == \[Gamma]/
2 - \[Gamma] Subscript[\[Rho], 2, 2][t] -
I (1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 3][t] -
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 2][t]) +
1/2 \[CapitalGamma] Subscript[\[Rho], 3, 3][t]},
{Derivative[1][Subscript[\[Rho], 2, 3]][t] ==
1/2 (-\[Gamma] Subscript[\[Rho], 2, 3][
t] - (\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 2, 3][
t]) - I (1/2 \[CapitalOmega]Ra Subscript[\[Rho], 2, 1][t] +
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 2][
t] + \[CapitalDelta]b Subscript[\[Rho], 2, 3][t] -
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 3][t])},
{Derivative[1][Subscript[\[Rho], 3, 1]][t] ==
1/2 (-\[Gamma] Subscript[\[Rho], 3, 1][
t] - (\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 3, 1][
t]) - I (-(1/2) \[CapitalOmega]Ra Subscript[\[Rho], 1, 1][
t] - 1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 1][
t] - \[CapitalDelta]a Subscript[\[Rho], 3, 1][t] +
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 3][t])},
{Derivative[1][Subscript[\[Rho], 3, 2]][t] ==
1/2 (-\[Gamma] Subscript[\[Rho], 3, 2][
t] - (\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 3, 2][
t]) - I (-(1/2) \[CapitalOmega]Ra Subscript[\[Rho], 1, 2][
t] - 1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 2][
t] - \[CapitalDelta]b Subscript[\[Rho], 3, 2][t] +
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 3][t])},
{Derivative[1][Subscript[\[Rho], 3, 3]][
t] == -I (-(1/2) \[CapitalOmega]Ra Subscript[\[Rho], 1, 3][t] -
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 3][t] +
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 1][t] +
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 2][
t]) - (\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 3, 3][t]}
}


while the boundary conditions are

{
{Subscript[\[Rho], 1, 1][0] == 1},
{Subscript[\[Rho], 1, 2][0] == 0},
{Subscript[\[Rho], 1, 3][0] == 0},
{Subscript[\[Rho], 2, 1][0] == 0},
{Subscript[\[Rho], 2, 2][0] == 0},
{Subscript[\[Rho], 2, 3][0] == 0},
{Subscript[\[Rho], 3, 1][0] == 0},
{Subscript[\[Rho], 3, 2][0] == 0},
{Subscript[\[Rho], 3, 3][0] == 0}
}


using the numerical values of

{
\[Gamma] = 0.1;

\[CapitalGamma] = 1;

\[CapitalOmega]Ra = 1;

\[CapitalOmega]Rb = 8;

\[CapitalDelta]a = 1;
\[CapitalDelta]b = 0.5;
}


These equations are the elements of the Von Neuman equation for three level atomic system. I want to solve the density matrix elements time dependently to see the time dependent behavior of electromagnetically induced transparency.

• Your ode is complex and the boundary conditions are real valued. Are you looking for real solution? Aug 26, 2022 at 8:33
• @Ulrich Neumann Thanks for your response. I am looking for complex solution. Later on i will use the real and imaginary parts to plot the susceptibility of the system. Regards Aug 26, 2022 at 9:15

Edit: Added code to generate 3D plots of selected solution as function of the parameters below

First initialize the constants then flatten the array of equations and initial conditions:

\[Gamma] = 0.1;
\[CapitalGamma] = 1;
\[CapitalOmega]Ra = 1;
\[CapitalOmega]Rb = 8;
\[CapitalDelta]a = 1;
\[CapitalDelta]b = 0.5;

theEqns = {{Derivative[1][Subscript[\[Rho], 1, 1]][
t] == \[Gamma] (-Subscript[\[Rho], 1, 1][t]) +
1/2 \[CapitalGamma] Subscript[\[Rho], 3, 3][t] -
I (1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 3][t] -
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 1][t]) + \[Gamma]/
2}, {Derivative[1][Subscript[\[Rho], 1, 2]][
t] == -(\[Gamma] Subscript[\[Rho], 1, 2][t]) -
I (\[CapitalDelta]a Subscript[\[Rho], 1, 2][
t] - \[CapitalDelta]b Subscript[\[Rho], 1, 2][t] -
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 2][t] +
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 1, 3][
t])}, {Derivative[1][Subscript[\[Rho], 1, 3]][t] ==
1/2 (-((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 1, 3][
t]) - \[Gamma] Subscript[\[Rho], 1, 3][t]) -
I (\[CapitalDelta]a Subscript[\[Rho], 1, 3][t] +
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 1][t] -
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 3][t] +
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 1, 2][
t])}, {Derivative[1][Subscript[\[Rho], 2, 1]][
t] == -(\[Gamma] Subscript[\[Rho], 2, 1][t]) -
I (-(\[CapitalDelta]a Subscript[\[Rho], 2, 1][
t]) + \[CapitalDelta]b Subscript[\[Rho], 2, 1][t] +
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 2, 3][t] -
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 1][
t])}, {Derivative[1][Subscript[\[Rho], 2, 2]][
t] == \[Gamma] (-Subscript[\[Rho], 2, 2][t]) +
1/2 \[CapitalGamma] Subscript[\[Rho], 3, 3][t] -
I (1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 3][t] -
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 2][t]) + \[Gamma]/
2}, {Derivative[1][Subscript[\[Rho], 2, 3]][t] ==
1/2 (-((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 2, 3][
t]) - \[Gamma] Subscript[\[Rho], 2, 3][t]) -
I (\[CapitalDelta]b Subscript[\[Rho], 2, 3][t] +
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 2, 1][t] +
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 2][t] -
1/
2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 3][
t])}, {Derivative[1][Subscript[\[Rho], 3, 1]][t] ==
1/2 (-((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 3, 1][
t]) - \[Gamma] Subscript[\[Rho], 3, 1][t]) -
I (-(\[CapitalDelta]a Subscript[\[Rho], 3, 1][t]) -
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 1][t] +
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 3][t] -
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 1][
t])}, {Derivative[1][Subscript[\[Rho], 3, 2]][t] ==
1/2 (-((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 3, 2][
t]) - \[Gamma] Subscript[\[Rho], 3, 2][t]) -
I (-(\[CapitalDelta]b Subscript[\[Rho], 3, 2][t]) -
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 2][t] -
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 2][t] +
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 3][
t])}, {Derivative[1][Subscript[\[Rho], 3, 3]][
t] == -((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 3, 3][
t]) - I (-(1/
2) (\[CapitalOmega]Ra Subscript[\[Rho], 1, 3][t]) +
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 1][t] -
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 3][t] +
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 2][t])}} // Flatten

theInit = {{Subscript[\[Rho], 1, 1][0] ==
1}, {Subscript[\[Rho], 1, 2][0] ==
0}, {Subscript[\[Rho], 1, 3][0] ==
0}, {Subscript[\[Rho], 2, 1][0] ==
0}, {Subscript[\[Rho], 2, 2][0] ==
0}, {Subscript[\[Rho], 2, 3][0] ==
0}, {Subscript[\[Rho], 3, 1][0] ==
0}, {Subscript[\[Rho], 3, 2][0] ==
0}, {Subscript[\[Rho], 3, 3][0] == 0}} // Flatten


Then use NDSolveValue to solve the system. Variable solArray contains the nine numerical solutions. First element in solArray, $$\texttt{solArray[[1]][t]}$$ is the first function $$\rho_{1,1}$$ and so forth. I use {t,0,1} as an example plotting the real components of the solutions in one plot then a single plot : For example, if you want the value of say $$\rho_{2,1}(1/2)$$ then use solArray[4][1/2].

solArray =
NDSolveValue[
Join[theEqns, theInit], {Subscript[\[Rho], 1, 1], Subscript[\[Rho],
1, 2], Subscript[\[Rho], 1, 3], Subscript[\[Rho], 2, 1],
Subscript[\[Rho], 2, 2], Subscript[\[Rho], 2, 3],
Subscript[\[Rho], 3, 1], Subscript[\[Rho], 3, 2],
Subscript[\[Rho], 3, 3]}, {t, 0, 1}];
Length@solArray
Plot[Re[#[t]] & /@ solArray, {t, 0, 1}, PlotRange -> All]
Plot[Re[solArray[[1]][t]], {t, 0, 1}]


Also, if the array of functions is a bit awkward to deal with, you can individually specify each function as the following and then use them separately (have to use different names on the left though):

{p11, p12, p13, p21, p22, p23, p31, p32, p33} =
NDSolveValue[
Join[theEqns, theInit], {Subscript[\[Rho], 1, 1], Subscript[\[Rho],
1, 2], Subscript[\[Rho], 1, 3], Subscript[\[Rho], 2, 1],
Subscript[\[Rho], 2, 2], Subscript[\[Rho], 2, 3],
Subscript[\[Rho], 3, 1], Subscript[\[Rho], 3, 2],
Subscript[\[Rho], 3, 3]}, {t, 0, 1}];
Length@solArray
Plot[Re@p11[t], {t, 0, 1}]


Edit: Added code to write DE as a function of parameters

First write the system as a function of parameters. Now the 9 equations are functions of the parameters of $$\texttt{theEqnsF}$$:

ClearAll["Global*"]
theEqns2[\[Gamma]_, \[CapitalGamma]_, \[CapitalOmega]Ra_, \
\[CapitalOmega]Rb_, \[CapitalDelta]a_, \[CapitalDelta]b_] := \
{{Derivative[1][Subscript[\[Rho], 1, 1]][
t] == \[Gamma] (-Subscript[\[Rho], 1, 1][t]) +
1/2 \[CapitalGamma] Subscript[\[Rho], 3, 3][t] -
I (1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 3][t] -
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 1][t]) + \[Gamma]/
2}, {Derivative[1][Subscript[\[Rho], 1, 2]][
t] == -(\[Gamma] Subscript[\[Rho], 1, 2][t]) -
I (\[CapitalDelta]a Subscript[\[Rho], 1, 2][
t] - \[CapitalDelta]b Subscript[\[Rho], 1, 2][t] -
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 2][t] +
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 1, 3][
t])}, {Derivative[1][Subscript[\[Rho], 1, 3]][t] ==
1/2 (-((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 1, 3][
t]) - \[Gamma] Subscript[\[Rho], 1, 3][t]) -
I (\[CapitalDelta]a Subscript[\[Rho], 1, 3][t] +
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 1][t] -
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 3][t] +
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 1, 2][
t])}, {Derivative[1][Subscript[\[Rho], 2, 1]][
t] == -(\[Gamma] Subscript[\[Rho], 2, 1][t]) -
I (-(\[CapitalDelta]a Subscript[\[Rho], 2, 1][
t]) + \[CapitalDelta]b Subscript[\[Rho], 2, 1][t] +
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 2, 3][t] -
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 1][
t])}, {Derivative[1][Subscript[\[Rho], 2, 2]][
t] == \[Gamma] (-Subscript[\[Rho], 2, 2][t]) +
1/2 \[CapitalGamma] Subscript[\[Rho], 3, 3][t] -
I (1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 3][t] -
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 2][t]) + \[Gamma]/
2}, {Derivative[1][Subscript[\[Rho], 2, 3]][t] ==
1/2 (-((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 2, 3][
t]) - \[Gamma] Subscript[\[Rho], 2, 3][t]) -
I (\[CapitalDelta]b Subscript[\[Rho], 2, 3][t] +
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 2, 1][t] +
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 2][t] -
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 3][
t])}, {Derivative[1][Subscript[\[Rho], 3, 1]][t] ==
1/2 (-((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 3, 1][
t]) - \[Gamma] Subscript[\[Rho], 3, 1][t]) -
I (-(\[CapitalDelta]a Subscript[\[Rho], 3, 1][t]) -
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 1][t] +
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 3][t] -
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 1][
t])}, {Derivative[1][Subscript[\[Rho], 3, 2]][t] ==
1/2 (-((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 3, 2][
t]) - \[Gamma] Subscript[\[Rho], 3, 2][t]) -
I (-(\[CapitalDelta]b Subscript[\[Rho], 3, 2][t]) -
1/2 \[CapitalOmega]Ra Subscript[\[Rho], 1, 2][t] -
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 2][t] +
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 3][
t])}, {Derivative[1][Subscript[\[Rho], 3, 3]][
t] == -((\[Gamma] + \[CapitalGamma]) Subscript[\[Rho], 3, 3][
t]) - I (-(1/2) (\[CapitalOmega]Ra Subscript[\[Rho], 1, 3][
t]) + 1/2 \[CapitalOmega]Ra Subscript[\[Rho], 3, 1][t] -
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 2, 3][t] +
1/2 \[CapitalOmega]Rb Subscript[\[Rho], 3, 2][t])}} // Flatten

theInit = {{Subscript[\[Rho], 1, 1][0] ==
1}, {Subscript[\[Rho], 1, 2][0] ==
0}, {Subscript[\[Rho], 1, 3][0] ==
0}, {Subscript[\[Rho], 2, 1][0] ==
0}, {Subscript[\[Rho], 2, 2][0] ==
0}, {Subscript[\[Rho], 2, 3][0] ==
0}, {Subscript[\[Rho], 3, 1][0] ==
0}, {Subscript[\[Rho], 3, 2][0] ==
0}, {Subscript[\[Rho], 3, 3][0] == 0}} // Flatten


Can now supply the system function $$\texttt{theEqnsF}$$ to NDSolveValue where now must supply values for the parameters:

   (* theEqns[\[Gamma],\[CapitalGamma],\[CapitalOmega]Ra,\[CapitalOmega]\
Rb,\[CapitalDelta]a,\[CapitalDelta]b] *)
{p11, p12, p13, p21, p22, p23, p31, p32, p33} =
NDSolveValue[
Join[theEqnsF[0.1, 1, 1, 8, 1, 0.5],
theInit], {Subscript[\[Rho], 1, 1], Subscript[\[Rho], 1, 2],
Subscript[\[Rho], 1, 3], Subscript[\[Rho], 2, 1],
Subscript[\[Rho], 2, 2], Subscript[\[Rho], 2, 3],
Subscript[\[Rho], 3, 1], Subscript[\[Rho], 3, 2],
Subscript[\[Rho], 3, 3]}, {t, 0, 1}];
Plot[Re@p33[t], {t, 0, 1}]


Now generate a series of plots for $$\texttt{p33}$$ when $$\Delta b$$ varies from 0.1 to 5:

plotTable = Table[
(* theEqns[\[Gamma],\[CapitalGamma],\[CapitalOmega]Ra,\
\[CapitalOmega]Rb,\[CapitalDelta]a,\[CapitalDelta]b] *)
{p11, p12, p13, p21, p22, p23, p31, p32, p33} =
NDSolveValue[
Join[theEqns[0.1, 1, 1, 8, 1, currentB],
theInit], {Subscript[\[Rho], 1, 1], Subscript[\[Rho], 1, 2],
Subscript[\[Rho], 1, 3], Subscript[\[Rho], 2, 1],
Subscript[\[Rho], 2, 2], Subscript[\[Rho], 2, 3],
Subscript[\[Rho], 3, 1], Subscript[\[Rho], 3, 2],
Subscript[\[Rho], 3, 3]}, {t, 0, 1}];
ParametricPlot3D[{t, currentB, Re@p33[t]}, {t, 0, 1}],
{currentB, 0.1, 5, 0.1}
];
selectedFunction = "p33";
Show[plotTable, BoxRatios -> {1, 1, 1}, PlotRange -> All,
AxesLabel -> {Style["t", 14], , Style["\[CapitalDelta]b", 14],
Style[selectedFunction, 14]}]


Finally generate a table of values for $$\texttt{p33}$$ in the form needed for the function ListPlot3D and generate a smooth function of $$\texttt{p33}$$ over the indicated range:

plotTable = Table[
(* theEqns[\[Gamma],\[CapitalGamma],\[CapitalOmega]Ra,\
\[CapitalOmega]Rb,\[CapitalDelta]a,\[CapitalDelta]b] *)
{p11, p12, p13, p21, p22, p23, p31, p32, p33} =
NDSolveValue[
Join[theEqns[0.1, 1, 1, 8, 1, currentB],
theInit], {Subscript[\[Rho], 1, 1], Subscript[\[Rho], 1, 2],
Subscript[\[Rho], 1, 3], Subscript[\[Rho], 2, 1],
Subscript[\[Rho], 2, 2], Subscript[\[Rho], 2, 3],
Subscript[\[Rho], 3, 1], Subscript[\[Rho], 3, 2],
Subscript[\[Rho], 3, 3]}, {t, 0, 1}];
if = p33;
Table[{t, currentB, if[t]} // N, {t, 0, 1, 1/50}],
{currentB, 0.1, 5, 0.1}
];

ListPlot3D[Flatten[plotTable, 1], BoxRatios -> {1, 1, 1},
AxesLabel -> {Style["t", 14], , Style["\[CapitalDelta]b", 14],
Style[selectedFunction, 14]}]
`

• Thank you for your reply and help. can we use different value of [CapitalDelta]b to plot the solution against time and [CapitalDelta]b Aug 26, 2022 at 17:46
• @Assad Hafiz: added some code to generate smooth funtion of a selected solution as a function of the parameters. You will need to review the code carefully and adjust it to your needs.
– josh
Aug 26, 2022 at 21:21
• Again thank you for your reply and help. Can we use finite difference method or runge kutta method to solve these equations. Aug 27, 2022 at 5:09