Matrix Integro-differential equation

Question

I'm trying to solve this integro-differential equation to obtain the density matrix elements $\rho$.

$ \frac{d\rho}{dt}=-i[H_{0}(t),\rho(t)]-A^{2}\Big[H_{1},\int_{t1}^{t}e^{-B(t-s)}\Big(e^{-i\int_{s}^{t}H_{0}(t')dt'}\Big)[H_{1},\rho(s)]\Big(e^{i\int_{s}^{t}H_{0}(t')dt'}\Big)ds\Big], $

with $$H_{0}(t)=\begin{pmatrix}v t - J\cos(k) & -{\it i}Jg\sin(k) \\{\it i}Jg\sin(k) & -v t + J\cos(k)\\\end{pmatrix},$$ and $$H_{1}=2X\left(\begin{array}{cc}1 & 0 \\0 & -1\\\end{array}\right)+2 Jg Y\sin(k)\left(\begin{array}{cc}0 & -i \\ i & 0\\\end{array}\right).$$

The parameters $A=1$, $B=1$, $v=5$, $g=1$, $J=1$, $X=1$, $Y=0.5$, $k=\pi/3$, $t_1=h_1/v$, and $t_2=h_2/v$ with $h_1=-10$ and $h_2=10$ where $t_1$ is initial time and $t_2$ is final time. At $t_1$ the initial density matrix is

$$\rho(t_1)=\left(\begin{array}{cc}1 & 0 \\ 0 & 0 \\ \end{array}\right).$$

The matrices $$R=e^{-{\cal i}\int_{s}^{t}H_{0}(t')dt'}$$ and $$R^\dagger=e^{{\cal i}\int_{s}^{t}H_{0}(t')dt'}$$ can be obtained analytically by Mathematica.

Without the second term in the right hand side of the differential equation (integration) the differential equation is easily solved by NDSolveValue. But in the presence of the integration I have to use the interpolation method to calculate the integration and then solve the integro-differential equation.

I used the method which has been used by Alex Trounev @Alex Trounev Solving an integro-differential equation with Mathematica.

Unfortunately my code doesn't work. I was wondering if you would be able to help me. Thank you

J = 1;
h1 = -10;
h2 = 10;
v = 1;
A = 1;
B = 1;
X = 1;
Y = 0.5;
k = \[Pi]/3;
t1 = h1/v;
t2 = h2/v;

del = 10^-6;

dt = (t1 - del)/6 - del;

n = 5;

H0[t_] = 2*{{v*t - J*Cos[k], -I*J*g*Sin[k]}, {I*J*g*Sin[k], -v*t + 
      J*Cos[k]}};

H1 = 2*X*{{1, 0}, {0, -1}} + 2*Y*J*g*Sin[k]*{{0, -I}, {I, 0}};

R = MatrixExp[-I*Integrate[H0[t'], {t', s, t}]];

Rdag = MatrixExp[I*Integrate[H0[t'], {t', s, t}]];

Int[0][t_] := {{0, 0}, {0, 0}};

Do[
 Sol[i] = 
  NDSolveValue[{D[rho[t], t] == -I*(H0[t].rho[t] - rho[t].H0[t]) - 
      A^2*(H1.Int[i - 1][t] - Int[i - 1][t].H1),
    rho[t1] == {{1, 0}, {0, 0}}}, rho, {t, t1, t2},
   Method -> {"PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"TensorProductGrid", 
         "MinPoints" -> 137}}}];
 
 Int[i] = 
  Interpolation[
   Flatten[ParallelTable[{t, 
      NIntegrate[
       e^(-B*(t - s))*(R.(H1.U[i][s] - U[i][s].H1).Rdag), {s, t1, t, 
        t2},
       Method -> "PrincipalValue"]}, {t, t1 + del, t - del, dt}], 
    1]];, {i, 1, n}]

Sol[5][t2][[1, 1]]

Sol[5][t2][[1, 2]]

Sol[5][t2][[2, 1]]

Sol[5][t2][[2, 2]]

Sol[5][t2][[2, 2]]

Answer 1

The problem can easily be solved after reformulation to a system of ODEs. Let us start with

$$ \frac{d\rho}{dt}=-i[H_{0}(t),\rho(t)]-A^{2}\Big[H_{1},\int_{t1}^{t}e^{-B(t-s)}\Big(e^{-i\int_{s}^{t}H_{0}(t')dt'}\Big)[H_{1},\rho(s)]\Big(e^{i\int_{s}^{t}H_{0}(t')dt'}\Big)ds\Big]. $$

We can denote $$ \Pi(t) = \int_{t1}^{t}e^{-B(t-s)}R(s,t)[H_{1},\rho(s)]R^\dagger(s,t)\,ds , $$ where

$$ R(s,t)=e^{-i\int_{s}^{t}H_{0}(t')dt'}. $$ and

$$ \Pi(t_1)=0. $$ It remains to devise an EOM for $\Pi(t)$. By differentiating over $t$ we obtain $$ \frac{d}{dt}\Pi(t)=[H_{1},\rho(t)]-B\,\pi(t)-i\Big[H_0(t), \Pi(t)\Big].\tag{1}\label{1} $$ Other equation to solve is: $$ \frac{d\rho}{dt}=-i[H_{0}(t),\rho(t)]-A^{2}\Big[H_{1},\Pi(t)\Big].\tag{2}\label{2}$$

The implementation should be straightforward, but let me know if you have problems.