# Matrix Integro-differential equation

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*{{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]]

• The method you try to implement with your code is not converge in this case since function rho is highly oscillating, and parameter A is not small. Can you give a link to the paper with this model explanation? Commented Jun 27, 2022 at 16:04
• Do you try to simulate dichotomous Markov noise effect on a quantum system? What this system is? Commented Jun 28, 2022 at 3:08

## 1 Answer

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.

• Thank you very much for your comments and note. I will leave message if I have any problem. Commented Jun 27, 2022 at 15:13