How to solve time dependent Optical Bloch Equations for a three level system?

Question

I am trying to solve the optical block equations for three level system which are time dependent, which i generated using the "Atomic Density Matrix" package.

When i use the command to solve this 9 linear coupled differential equations with some initial conditions it is showing me no enough memory is available.

So what should be the system requirements to solve them Analytically ?? i have to increase the RAM or hard disk memory or anything else ?

Similarly i want to work out for five level system where there will be 25 such equations.

here are the equations with initial conditions.

variable = {Subscript[ρ, 1, 1][t], Subscript[ρ, 1, 2][t],Subscript[ρ, 1, 3][t], Subscript[ρ, 2, 1][t], Subscript[ρ, 2, 2][t], Subscript[ρ, 2, 3][t],Subscript[ρ, 3, 1][t], Subscript[ρ, 3, 2][t], Subscript[ρ, 3, 3][t]};


init = {Subscript[ρ, 1, 1][0] == 1, Subscript[ρ, 1, 2][0] == 0, Subscript[ρ, 1, 3][0] == 0, Subscript[ρ, 2, 1][0] == 0, Subscript[ρ, 2, 2][0] == 0, Subscript[ρ, 2, 3][0] == 0, Subscript[ρ, 3, 1][0] == 0, Subscript[ρ, 3, 2][0] == 0, Subscript[ρ, 3, 3][0] == 0};



 equations = {Derivative[1][Subscript[ρ, 1, 1]][t] == γ/
 2 - γ Subscript[ρ, 1, 1][t] - 
 I (1/2 ΩRa Subscript[ρ, 1, 3][t] - 
    1/2 ΩRa Subscript[ρ, 3, 1][t]) + 
 1/2 Γ Subscript[ρ, 3, 3][t], 

   Derivative[1][Subscript[ρ, 1, 2]][
     t] == -γ Subscript[ρ, 1, 2][t] - 
     I (Δa Subscript[ρ, 1, 2][
          t] - Δb Subscript[ρ, 1, 2][t] + 
        1/2 ΩRb Subscript[ρ, 1, 3][t] - 
        1/2 ΩRa Subscript[ρ, 3, 2][t]), 
   Derivative[1][Subscript[ρ, 1, 3]][t] == 
    1/2 (-γ Subscript[ρ, 1, 3][
          t] - (γ + Γ) Subscript[ρ, 1, 3][
          t]) - I (1/2 ΩRa Subscript[ρ, 1, 1][t] + 
        1/2 ΩRb Subscript[ρ, 1, 2][
          t] + Δa Subscript[ρ, 1, 3][t] - 
        1/2 ΩRa Subscript[ρ, 3, 3][t]), 
   Derivative[1][Subscript[ρ, 2, 1]][
     t] == -γ Subscript[ρ, 2, 1][t] - 
     I (-Δa Subscript[ρ, 2, 1][
          t] + Δb Subscript[ρ, 2, 1][t] + 
        1/2 ΩRa Subscript[ρ, 2, 3][t] - 
        1/2 ΩRb Subscript[ρ, 3, 1][t]), 
   Derivative[1][Subscript[ρ, 2, 2]][t] == γ/
     2 - γ Subscript[ρ, 2, 2][t] - 
     I (1/2 ΩRb Subscript[ρ, 2, 3][t] - 
        1/2 ΩRb Subscript[ρ, 3, 2][t]) + 
     1/2 Γ Subscript[ρ, 3, 3][t], 
   Derivative[1][Subscript[ρ, 2, 3]][t] == 
    1/2 (-γ Subscript[ρ, 2, 3][
          t] - (γ + Γ) Subscript[ρ, 2, 3][
          t]) - I (1/2 ΩRa Subscript[ρ, 2, 1][t] + 
        1/2 ΩRb Subscript[ρ, 2, 2][
          t] + Δb Subscript[ρ, 2, 3][t] - 
        1/2 ΩRb Subscript[ρ, 3, 3][t]), 
   Derivative[1][Subscript[ρ, 3, 1]][t] == 
    1/2 (-γ Subscript[ρ, 3, 1][
          t] - (γ + Γ) Subscript[ρ, 3, 1][
          t]) - I (-(1/2) ΩRa Subscript[ρ, 1, 1][
          t] - 1/2 ΩRb Subscript[ρ, 2, 1][
          t] - Δa Subscript[ρ, 3, 1][t] + 
        1/2 ΩRa Subscript[ρ, 3, 3][t]), 
   Derivative[1][Subscript[ρ, 3, 2]][t] == 
    1/2 (-γ Subscript[ρ, 3, 2][
          t] - (γ + Γ) Subscript[ρ, 3, 2][
          t]) - I (-(1/2) ΩRa Subscript[ρ, 1, 2][
          t] - 1/2 ΩRb Subscript[ρ, 2, 2][
          t] - Δb Subscript[ρ, 3, 2][t] + 
        1/2 ΩRb Subscript[ρ, 3, 3][t]), 
   Derivative[1][Subscript[ρ, 3, 3]][
     t] == -I (-(1/2) ΩRa Subscript[ρ, 1, 3][t] - 
        1/2 ΩRb Subscript[ρ, 2, 3][t] + 
        1/2 ΩRa Subscript[ρ, 3, 1][t] + 
        1/2 ΩRb Subscript[ρ, 3, 2][
          t]) - (γ + Γ) Subscript[ρ, 3, 3][
       t]};
         TableForm[solution = DSolve[Join[equations, init], variable, t]]

If it has to be solved numerically then how to do that ?

can any one help me out with this ?

Answer 1

Edit

Because of the question's title, I assumed initially that you're interested in the general case of time-dependent driving. This complication is briefly addressed in the first part of my answer.

However, the actual system of differential equations in the question has purely time-independent coefficients and is therefore easier to solve in principle. The only major issue is the appearance of many parameters which Mathematica by default assumes to be general complex numbers satisfying no special assumptions. By making the parameters numerical, you get the solutions easily, as shown in the second part below.

End edit

Here is a simpler example to show that Mathematica can at least get analytical results for the undamped periodically driven two-level system. I compare the result of DSolve with the textbook Rabi formula for the occupation probability of the excited state at the end:

Clear[e0, v0, ω, t, c1, c2];

$Assumptions = {e0 > 0, v0 > 0, ω > 0, t > 0};

h0 = ( {
    {-e0/2, 0},
    {0, e0/2}
   } );

v = ( {
    {0, v0 E^(I ω t)},
    {v0 E^(-I ω t), 0}
   } );

vi = Simplify[MatrixExp[I h0 t].v.MatrixExp[-I h0 t]];

ψ[t_] = {c1[t], c2[t]} /. 
   Simplify@
    First@DSolve[
      Join[Thread[
        vi.{c1[t], c2[t]} == I {c1'[t], c2'[t]}], {c1[0] == 1, 
        c2[0] == 0}], {c1[t], c2[t]}, t];

FullSimplify[
 ComplexExpand[Last[Conjugate[ψ[t]] ψ[t]]] == 
  v0^2/(v0^2 + (ω - e0)^2/4)
    Sin[Sqrt[v0^2 + (ω - e0)^2/4] t]^2]

(* ==> True *)

The last line means that the analytical result derived by DSolve agrees with a know textbook formula.

Here I defined the time-dependent potential v and then transformed it to the interaction picture vi. The coefficients c1, c2 are the interaction-picture amplitudes, and their absolute square yields the occupation probabilities. If I introduce a damping term into this two-level equation, DSolve no longer yields a solution.

This means that for a more general problem than what I'm doing here, it's likely that the best bet is a numerical approach as described in Solving a time-dependent Schroedinger equation.

Addressing the equations in the question

You only need to perform some small additions to your original code to get the numerical result:

ΩRa = 1;
ΩRb = 1.1;
Δa = .1;
Δb = -.1;
Γ = .01;
γ = .001;

NDSolve[Join[equations, init], variable, {t, 0, 1}]

Here I defined all the variable parameters of the problem by assigning them (completely arbitrary) numbers. This is of course where you have to provide additional information to make the correct choice. Then in NDSolve you need to specify a time interval, and you're done.

Edit: Analytical results with numerical coefficients

Because there is no driving term, DSolve actually does produce a result if you just follow the last suggestion and specify all the parameters in terms of numerical values. Here is what that looks like:

solution = DSolve[Join[equations, init], variable, t];
Short[solution[[1,1]],10]

Some of the parameters were specified as machine-precision numbers, so you get corresponding numerical factors in the output as well. One could in principle try with exact numbers as input for DSOlve, but that runs much too slowly so I didn't wait for it to complete.

The fact that DSolve can solve this system also suggests that one could get it to derive a general result with symbolic parameters, but that requires tweaking the Assumptions for those parameters, which is not always successful.

The solution in the present case has the simple form of a matrix exponential, but that exponential is difficult to compute symbolically because of the many parameters in the matrix. However, you can make use of the formal solution to define a straightforward function that will spit out a solution for arbitrary initial conditions as follows:

Clear[ΩRa, ΩRb, Δa, Δb, Γ, γ];

j = D[equations[[All, 2]], {variable}];

initVec = init[[All, 2]]

(* ==> {1, 0, 0, 0, 0, 0, 0, 0, 0} *)

liouville[ΩRa_, ΩRb_, Δa_, Δb_, Γ_, γ_] := MatrixExp[t j];

ΩRa = 1;
ΩRb = 1.1;
Δa = .1;
Δb = -.1;
Γ = .01;
γ = .001;

sol = 
  liouville[ΩRa, ΩRb, Δa, Δb, Γ, γ].initVec;

Plot[Evaluate[Abs@sol], {t, 0, 10}, PlotRange -> {0, 1}]

Here I first defined the matrix of coefficients j for your linear, constant-coefficient differential equation in 9 variables, by taking derivatives with respect to variable. Then I define a function liouville that performs the matrix exponential if it is provided with the necessary parameters. Here you will have to insert the numerical parameters to get an actual result. Then the time evolution is simply obtained by multiplying the vector of initial conditions initVec with this matrix. The plot shows what the result looks like for the same choice of parameters I used above.