0
$\begingroup$

I have searched almost all similar pages discussing a system of differential equations, but none of them discussed the case of variable coefficients. My problem is a bit more complicated, because instead of functions we have quantum mechanical operators (say, $o_1(t)$ or $o_1^\dagger(t)$, which is the hermitian conjugate of $o_1(t)$. I have shown $\dagger$ and $*$ signs by subscript $s$ in the following codes). One cannot assign an initial value to $o_1$, but it is meaningful to have expressions like $\langle o_1^\dagger(0) o_1(0)\rangle$, which depends on the quantum mechanical state. In other words, $\langle o_1^\dagger(0) o_1(0)\rangle$ is known, but $o_1(0)$ is unknown.

Now consider a system of differential equaitons like enter image description here where

 m11 = -0.25`;
 m12 = 0;
 m13 = 0.75` E^((0.` + 22.8` I) t);
 m14 = 0;
 m15 = (0.097` + 0.193` E^((0.` - 22.8` I) t) - 
 0.29` E^((0.` - 45.6` I) t));
 m16 = (-0.29` + 0.193` E^((0.` + 22.8` I) t) + 
 0.097` E^((0.` + 45.6` I) t));m21 = 0;
 m22 = -0.25`;
 m23 = 0;
 m24 =  0.75` E^((0.` - 22.8` I) t) ; m25 = (-0.29` + 
0.193`  E^((0.` - 22.8` I) t) + 0.0967` E^((0.` - 45.6` I) t));
m26 = (0.097` + 0.193` E^((0.` + 22.8` I) t) - 
0.29` E^((0.` + 45.6` I) t));
m31 = 0.75` E^((0.` - 22.8` I) t);
m32 = 0;
m33 = -0.25`;
m34 = 0;
m35 = (0.097` E^((0.` - 22.8` I) t) + 0.193` E^((0.` - 45.6` I) t) - 
0.29` E^((0.` - 68.4` I) t));
m36 = (0.193` - 0.29` E^((0.` - 22.8` I) t) + 
0.097` E^((0.` + 22.8` I) t));
m41 = 0.75` E^((0.` + 22.8` I) t) ;
m42 = 0;
m43 = -0.25`;
m44 = 0;
m45 = (0.193` + 0.097` E^((0.` - 22.8` I) t) - 
0.29` E^((0.` + 22.8` I) t));
m46 = (0.097` E^((0.` + 22.8` I) t) + 0.193` E^((0.` + 45.6` I) t) - 
0.29` E^((0.` + 68.4` I) t));
m51 = (-0.097` - 0.193` E^((0.` + 22.8` I) t) + 
0.29` E^((0.` + 45.6` I) t));
m52 = (-0.29` + 0.193` E^((0.` + 22.8` I) t) + 
0.097` E^((0.` + 45.6` I) t));
m53 = (-0.097` E^((0.` + 22.8` I) t) - 0.193` E^((0.` + 45.6` I) t) + 
0.29` E^((0.` + 68.4` I) t));
m54 = (0.193` - 0.29` E^((0.` - 22.8` I) t) + 
0.097` E^((0.` + 22.8` I) t));
m55 = -0.037`;
m56 = 0;
m61 = (-0.29` + 0.193` E^((0.` - 22.8` I) t) + 
0.097` E^((0.` - 45.6` I) t));
m62 = (0.097` + 0.193` E^((0.` - 22.8` I) t) - 
0.29` E^((0.` - 45.6` I) t)) ;
m63 = (0.193` + 0.097` E^((0.` - 22.8` I) t) - 
0.29` E^((0.` + 22.8` I) t));
m64 = (0.097` E^((0.` - 22.8` I) t) + 0.193` E^((0.` - 45.6` I) t) - 
0.29` E^((0.` - 68.4` I) t));
m65 = 0;
m66 = -0.037`;

 \[Lambda]1 = (0.` - 
 2907.509` I) + (0.` + 
  5330.436` I) E^((0.` - 22.8` I) t) - (0.` + 
  2422.925 I) E^((0.` + 22.8` I) t);

\[Lambda]1s = (0.` + 
 2907.511` I) - (0.` + 5330.437` I) E^((0.` + 22.8 I) t) + (0.` + 
  2422.925` I) E^((0.` - 22.8 I) t);

\[Lambda]2 = (0.` - 
 1453.754 I) - (0.` + 4845.851` I) E^((0.` - 22.8` I) t) + (0.` + 
  6299.608` I) E^((0.` - 45.6` I) t);

\[Lambda]2s = (0.` + 
 1453.755 I) + (0.` + 1356.17` I) E^((0.` + 22.8` I) t) - (0.` + 
  6299.608` I) E^((0.` + 45.6` I) t);

\[Lambda]3 = (0.` - 
 1502.866` I) - (0.` + 
  1127.149` I) E^((0.` - 22.8` I) t) + (0.` + 
  5260.031` I) E^((0.` + 22.8` I) t) - (0.` + 
  1502.866` I) E^((0.` + 45.6` I) t) - (0.` + 
  1127.149` I) E^((0.` + 68.4` I) t);

\[Lambda]3s = (0.` + 
 1502.866` I) - (0.` + 
  5260.031` I) E^((0.` - 22.8` I) t) + (0.` + 
  1127.149` I) E^((0.` + 22.8` I) t) + (0.` + 
  1502.866` I) E^((0.` - 45.6` I) t) + (0.` + 
  1127.149` I) E^((0.` - 68.4` I) t);

n1 = +0.707` E^((0.` + 34.2` I) t) fas + E^((0.` + 34.2` I) t) fp;
n1s = 0.707` E^((0.` - 34.2` I) t) fa + E^((0.` - 34.2` I) t) fps;
n2 = 0.707` E^((0.` + 11.4` I) t) fas - E^((0.` + 11.4` I) t) fp;
n2s = 0.707` E^((0.` - 11.4` I) t) fa - E^((0.` - 11.4` I) t) fps;
n3 = 0.272` E^((0.` + 22.8` I) t) fm;
n3s = 0.272` E^((0.` - 22.8` I) t) fms;
M = {{m11, m12, m13, m14, m15, m16}, {m21, m22, m23, m24, m25, 
m26}, {m31, m32, m33, m34, m35, m36}, {m41, m42, m43, m44, m45, 
m46}, {m51, m52, m53, m54, m55, m56}, {m61, m62, m63, m64, m65, 
m66}};

\[Lambda] = {{\[Lambda]1}, {\[Lambda]1s}, {\[Lambda]2}, \
{\[Lambda]2s}, {\[Lambda]3}, {\[Lambda]3s}};

n = {{n1}, {n1s}, {n2}, {n2s}, {n3}, {n3s}};

In page 8 of https://arxiv.org/pdf/1609.00075.pdf it is suggested to solve the system by the matrix exponential: enter image description here

In general, it takes very long for Mathematica to calculate the matrix exponential. So, it is suggested to expand $\mathcal T \int_{0}^{t} d\tau M(\tau)$ as follows: enter image description here

where the time interval is divided into $N$ pieces: $h=(t-0)/N$. To ensure the time ordering, $\mathcal T$, the lower limit of the product is $N-1$ and the upper is $0$. The reference has defined enter image description here As you see, the lower limit of integral is $\tau$, not $0$. The whole problem is summarized in finding the above matrix and then calculating expressions like $A=\int_{0}^{t}d\tau d_{11}(t,\tau) \lambda_1(\tau)$. Finally, one can plot $A^*A$ versus time. I have no idea how to divide the exponential integral where the lower limit $\tau$ is another integration variable. As you see, we don't need to know the initial value of $o_1$ and so on. As one can find $d_{ij}(t,\tau)$ and then integrate with respect to $\tau$, everything is done.

Finally, I should emphasize that one can find the exponential integral analytically, but in that case one have to use MatrixExp which for some parameters takes very long.

$\endgroup$
  • 1
    $\begingroup$ So... where is the question? $\endgroup$ – Henrik Schumacher Jan 13 at 6:52
  • $\begingroup$ @HenrikSchumacher I have no idea how to write a code that calculates the above equaion (exponential integral), knowing M. I have tried MatrixExp, but in general it doesn't work well. $\endgroup$ – Saeid Jan 13 at 14:06
  • $\begingroup$ Your code is lacking the first row of $M$. $\endgroup$ – Cesareo Jan 13 at 18:17
  • $\begingroup$ @Cesareo I fixed that. $\endgroup$ – Saeid Jan 13 at 20:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.