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|}[H_{1},\rho(s)]ds\Big], $

with $$H_{0}(t)=\begin{pmatrix}v t & -{\it i}Jg \\{\it i}Jg & -v t \\\end{pmatrix},$$ and $$H_{1}=2X\left(\begin{array}{cc}1 & 0 \\0 & -1\\\end{array}\right)+2 Jg Y\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).$$

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{{vt - JCos[k], -IJgSin[k]}, {IJgSin[k], -vt + 
      J*Cos[k]}};
H1 = 2X{{1, 0}, {0, -1}} + 2YJgSin[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.