Appendix 1
The perturbed Hamiltonian for the system is given by
$$\begin{aligned} H(t)= \frac{(p+eA)^2+q^2}{2}-e\varPhi = \frac{p^2+q^2}{2}+ \frac{e}{2}\big ((p,A)+(A,p)\big )-e\varPhi + \frac{{e}^2A^2}{2} \end{aligned}$$
(34)
Since
$$\begin{aligned} (p,A)+(A,p)&=-2iA\cdot \nabla -i(\nabla \cdot A)\nonumber \\ (\nabla \cdot A)&=\frac{1}{2}\nabla \cdot (B(t)\times {r}) =-\frac{1}{2}B(t)\cdot (\nabla \times {r})=0, \end{aligned}$$
(35)
we get
$$\begin{aligned} H(t)&=H_0+eA\cdot p+\frac{{e}^2A^2}{2}-e\varPhi =H_0+e\left( \frac{1}{2}B(t)\times r\cdot p+E(t)\cdot r\right) \nonumber \\&\quad +\,\frac{{e}^2}{8}\left( B(t)\times r\right) ^2 \end{aligned}$$
(36)
Since
$$\begin{aligned} \left( B(t)\times r\right) ^2=B^2(t)q^2-(B(t),q)^2 \end{aligned}$$
(37)
and with \( L=r\times p \triangleq q\times p\), we obtain
$$\begin{aligned} H(t)=H_0+e\bigg (\frac{1}{2}B(t)\cdot L+E(t)\cdot q\bigg )+\frac{{e}^2}{8} \left( B^2(t)q^2-(B(t),q)^2\right) \end{aligned}$$
(38)
Appendix 2
Referring to Eq. (14)
$$\begin{aligned} W(t)&=I- ie\int _0^t\left( \frac{1}{2}B\big (t_1\big )\cdot L\big (t_1\big )+E\big (t_1\big )\cdot q\big (t_1\big )\right) {\text {d}}t_1\nonumber \\&\quad - i\frac{{e}^2}{8}\int _0^t\left( B\big (t_1\big )\times q\big (t_1\big )\right) ^2{\text {d}}t_1\nonumber \\&\quad -\,\nonumber {e}^2 \int _{0<t_2<t_1<t}\left( \frac{1}{2}B\big (t_1\big )\cdot L\big (t_1\big )+E\big (t_1\big )\cdot q\big (t_1\big ) \right) \nonumber \\&\quad \left( \frac{1}{2}B\big (t_2\big )\cdot L\big (t_2\big )+E\big (t_2\big )\cdot q\big (t_2\big )\right) {\text {d}}t_1{\text {d}}t_2+\,O\left( {e}^3\right) \end{aligned}$$
(39)
Denoting
$$\begin{aligned} \langle n_1n_2n_3|L_1(t)|m_1m_2m_3\rangle \triangleq \langle n|L_1(t)|m\rangle =\mathrm{e}^{itE(n,m)}\langle n|L_1|m\rangle \end{aligned}$$
(40)
where
$$\begin{aligned} E(n,m)=E(n)-E(m)=sum_{k=1}^3{(n_k-m_k)} \end{aligned}$$
To get general element for the angular momentum matrix, we have
$$\begin{aligned} \langle n|L_1|m\rangle =i\langle n|a_2a_3^\dag -a_2^\dag a_3|m\rangle \end{aligned}$$
and
$$\begin{aligned} a_3|m_1m_2m_3\rangle= & {} \sqrt{m}_3|m_1,m_2,m_3-1\rangle \\ \langle n_1n_2n_3|a_2^\dag= & {} \sqrt{n}_2\langle n_1,n_2-1,n_3|\\ \langle n|a_2^\dag a_3|m\rangle= & {} \sqrt{n_2m_3}\delta \big [n_1-m_1\big ]\delta \big [n_2-1-m_2\big ]\delta \big [n_3-m_3+1\big ]\\ \langle n|a_2a_3^\dag |m\rangle= & {} \langle n|a_3^\dag a_2|m\rangle =\sqrt{n_3m_2}\delta \big [n_1-m_1\big ]\delta \big [n_2-m_2+1\big ]\delta \big [n_3-1-m_3\big ] \end{aligned}$$
For the concise formulation of the equations, we define an operator \({\mathcal {Z}}^{-1}_k\), which acts on function \(f:{\mathbb {Z}}^3_+\rightarrow {\mathbb {R}}\) by the rule
$$\begin{aligned} {\mathcal {Z}}^{-1}_1 f\big (n_1n_2n_3\big )= & {} f\big (n_1-1,n_2,n_3\big )\\ {\mathcal {Z}}^{-1}_2 f\big (n_1n_2n_3\big )= & {} f\big (n_1,n_2-1,n_3\big )\\ {\mathcal {Z}}^{-1}_3 f\big (n_1n_2n_3\big )= & {} f\big (n_1,n_2,n_3-1\big ) \end{aligned}$$
Thus, we have
$$\begin{aligned} \langle n|L_1|m\rangle =-i\left( \sqrt{n_2m_3}{\mathcal {Z}}^{-1}_2{\mathcal {Z}}_3-\sqrt{n_3m_2}{\mathcal {Z}}_2{\mathcal {Z}}^{-1}_3\right) \delta [n-m] \end{aligned}$$
Likewise,
$$\begin{aligned} \langle n|L_2|m\rangle= & {} -i\left( \sqrt{n_3m_1}{\mathcal {Z}}^{-1}_3{\mathcal {Z}}_1-\sqrt{n_1m_3}{\mathcal {Z}}_3{\mathcal {Z}}^{-1}_1\right) \delta [n-m]\\ \langle n|L_3|m\rangle= & {} -i\left( \sqrt{n_1m_2}{\mathcal {Z}}^{-1}_1{\mathcal {Z}}_2-\sqrt{n_2m_1}{\mathcal {Z}}_1{\mathcal {Z}}^{-1}_2\right) \delta [n-m] \end{aligned}$$
The general expression can be written as
$$\begin{aligned} \langle n|L_k|m \rangle = -i\sum _{r,s} \epsilon ({\text {krs}}) \sqrt{n_rm_s}{\mathcal {Z}}^{-1}_r{\mathcal {Z}}_s \delta [n-m] \end{aligned}$$
We also need the matrix elements of \(B(t)\cdot L(t)\)
$$\begin{aligned} \langle n|B(t)\cdot L(t)|m\rangle&=-\iota \text {e}^{\iota tE(n-m)}\left\{ B_1(t)\left( \sqrt{n_2m_3}{\mathcal {Z}}^{-1}_2{\mathcal {Z}}_3-\sqrt{n_3m_2}{\mathcal {Z}}_2{\mathcal {Z}}^{-1}_3\right) \right. \nonumber \\&\quad +\, B_2(t)\left( \sqrt{n_3m_1}{\mathcal {Z}}^{-1}_3{\mathcal {Z}}_1-\sqrt{n_1m_3}{\mathcal {Z}}_3{\mathcal {Z}}^{-1}_1\right) \nonumber \\&\quad \left. +\, B_3(t)\left( \sqrt{n_1m_2}{\mathcal {Z}}^{-1}_1{\mathcal {Z}}_2-\sqrt{n_2m_1}{\mathcal {Z}}_1{\mathcal {Z}}^{-1}_2\right) \right\} \delta [n-m] \end{aligned}$$
(41)
To get the general matrix elements for position matrix, we have
$$\begin{aligned} \langle n|q_k(t)|m\rangle&= \mathrm{e}^{itE(n,m)}\langle n|q_k|m\rangle \nonumber \\ \langle n|q_1|m\rangle&=\langle n|\frac{\left( a_1+a_1^\dag \right) }{\sqrt{2}}|m\rangle =\left( \sqrt{\frac{m_1}{2}} {\mathcal {Z}}_1+ \sqrt{\frac{n_1}{2}}{\mathcal {Z}}^{-1}_1\right) \delta [n-m] \end{aligned}$$
(42)
and likewise for \(q_2\) and \(q_3\). The general expression comes out to be,
$$\begin{aligned} \langle n|q_k|m\rangle =\left( \sqrt{\frac{m_k}{2}}{\mathcal {Z}}_k+ \sqrt{\frac{n_k}{2}}{\mathcal {Z}}^{-1}_k\right) \delta [n-m],\quad k=1,2,3 \end{aligned}$$
We also need matrix elements of
$$\begin{aligned} \langle n|q_k(t)q_l(t)|m\rangle =\mathrm{e}^{itE(n,m)}\langle n|q_kq_l|m\rangle \end{aligned}$$
For \(k \ne l\),
$$\begin{aligned}&\langle n|q_kq_l|m\rangle =\langle n|\frac{\left( a_k+a_k^\dag \right) }{\sqrt{2}} \frac{\left( a_l+a_l^\dag \right) }{\sqrt{2}}|m \rangle =\frac{1}{2}\left\{ \langle n|a_ka_l +a_ka_l^\dag +a_k^\dag a_l + a_k^\dag a_l^\dag |m \rangle \right\} \nonumber \\&\quad =\frac{1}{2}\left\{ \sqrt{m_km_l}{\mathcal {Z}}_k{\mathcal {Z}}_l + \sqrt{m_kn_l}{\mathcal {Z}}_k{\mathcal {Z}}^{-1}_l +\sqrt{m_ln_k}{\mathcal {Z}}^{-1}_k{\mathcal {Z}}_l +\sqrt{n_kn_l}{\mathcal {Z}}^{-1}_k{\mathcal {Z}}^{-1}_l\right\} \delta [n-m] \end{aligned}$$
For \(k=l\), the required matrix elements is given by
$$\begin{aligned} \langle n|q_k^2|m\rangle&=\frac{1}{2}\langle n|\left( a_k+a_k^\dag \right) ^2|m\rangle =\frac{1}{2}\left\{ \langle n|a_k^2 +a_k^{\dag ^2} +a_ka_k^\dag +a_k^\dag a_k|m \rangle \right\} \\&=\frac{1}{2}\left\{ \sqrt{m_k(m_k-1)}{\mathcal {Z}}_k^2+ \sqrt{n_k(n_k-1)}{\mathcal {Z}}^{-2}_k+(2m_k+1)\right\} \delta [n-m] \end{aligned}$$
We redefine the Dyson series as follows,
$$\begin{aligned} W(t)=I+eW_1(t)+{e}^2W_2(t)+O\left( {e}^3\right) \end{aligned}$$
where
$$\begin{aligned} W_1(t)=-i\int _0^t\left( \frac{1}{2}B\big (t_1\big )\cdot L\big (t_1\big )+E\big (t_1\big )\cdot q\big (t_1\big )\right) {\text {d}}t_1 \end{aligned}$$
and
$$\begin{aligned} \langle n|W_1(t)|m \rangle&=-i\int _0^t\frac{1}{2}B\big (t_1\big )\cdot \langle n|L\big (t_1\big )|m \rangle {\text {d}}t_1 -i\int _0^t E\big (t_1\big )\cdot \langle n|q\big (t_1\big )|m \rangle {\text {d}}t_1\nonumber \\&=- \frac{i}{2}\int _0^t \mathrm{e}^{it_1E(n,m)t_1}\left\{ B_1(t)\left( \sqrt{n_2m_3} {\mathcal {Z}}^{-1}_2{\mathcal {Z}}_3-\sqrt{n_3m_2} {\mathcal {Z}}_2{\mathcal {Z}}^{-1}_3\right) \right. \nonumber \\&\quad +\,B_2(t)\left( \sqrt{n_3m_1}{\mathcal {Z}}^{-1}_3 {\mathcal {Z}}_1-\sqrt{n_1m_3}{\mathcal {Z}}_3 {\mathcal {Z}}^{-1}_1\right) \nonumber \\&\quad \left. +\,B_3(t)\left( \sqrt{n_1m_2}{\mathcal {Z}}^{-1}_1{\mathcal {Z}}_2- \sqrt{n_2m_1}{\mathcal {Z}}_1{\mathcal {Z}}^{-1}_2\right) \right\} \delta [n-m]{\text {d}}t_1\nonumber \\&\quad -\,i\int _0^t \mathrm{e}^{it_1E(n,m)_1}\left\{ E_1\big (t_1\big )\left( \sqrt{\frac{m_1}{2}}{\mathcal {Z}}_1+ \sqrt{\frac{n_1}{2}}{\mathcal {Z}}^{-1}_1\right) \right. \nonumber \\&\quad +\,E_2\big (t_1\big )\left( \sqrt{\frac{m_2}{2}}{\mathcal {Z}}_2+ \sqrt{\frac{n_2}{2}}{\mathcal {Z}}^{-1}_2\right) \nonumber \\&\quad \left. +\,E_3\big (t_1\big )\left( \sqrt{\frac{m_3}{2}}{\mathcal {Z}}_3+ \sqrt{\frac{n_3}{2}}{\mathcal {Z}}^{-1}_3\right) \right\} \delta [n-m]{\text {d}}t_1 \end{aligned}$$
(43)
or on equivalently, defining
$$\begin{aligned} \hat{B_k}(n,t)=\int _0^t B_k\big (t_1\big )\text {e}^{\iota n t_1}{\text {d}}t_1 \end{aligned}$$
and
$$\begin{aligned} \hat{E_k}(n,t)=\int _0^t E_k\big (t_1\big )\text {e}^{\iota n t_1}{\text {d}}t_1 \end{aligned}$$
We get
$$\begin{aligned} \langle n|W_1(t)|m \rangle&= -\frac{i}{2}\left\{ \hat{B_1}(n-m,t)\left( \sqrt{n_2m_3}{\mathcal {Z}}^{-1}_2{\mathcal {Z}}_3- \sqrt{n_3m_2}{\mathcal {Z}}_2{\mathcal {Z}}^{-1}_3\right) \right. \nonumber \\&\quad +\,\hat{B_2}(n-m,t)\left( \sqrt{n_3m_1}{\mathcal {Z}}^{-1}_3{\mathcal {Z}}_1- \sqrt{n_1m_3}{\mathcal {Z}}_3{\mathcal {Z}}^{-1}_1\right) \nonumber \\&\quad \left. +\,\hat{B_3}(n-m,t)\left( \sqrt{n_1m_2}{\mathcal {Z}}^{-1}_1{\mathcal {Z}}_2- \sqrt{n_2m_1}{\mathcal {Z}}_1{\mathcal {Z}}^{-1}_2\right) \right\} \delta [n-m]\nonumber \\&\quad -\,i\sum _{k=1}^3 \hat{E_k}(n-m,t)\left( \sqrt{\frac{m_k}{2}}{\mathcal {Z}}_k+ \sqrt{\frac{n_k}{2}}{\mathcal {Z}}^{-1}_k\right) \delta [n-m], \end{aligned}$$
(44)
and also \(W_2(t)\) is given by
$$\begin{aligned} W_2(t)&=-\frac{i}{8}\int _0^t(B\big (t_1\big )\times q\big (t_1\big ))^2{\text {d}}t_1\\&-\int _{0<t_2<t_1<t}\bigg (\frac{1}{2}B\big (t_1\big )\cdot L\big (t_1\big )+ E\big (t_1\big )\cdot q\big (t_1\big )\bigg )\\&\bigg (\frac{1}{2}B\big (t_2\big )\cdot L\big (t_2\big )+E\big (t_2\big )\cdot q\big (t_2\big )\bigg ){\text {d}}t_1{\text {d}}t_2 \end{aligned}$$
On using Eq. (37), we get
$$\begin{aligned}&\langle n|(B(t)\times q)^2|m \rangle =\text {e}^{\iota t E(n,m)}\bigg \{B^2(t)\sum _{k=1}^3 \langle n|q_k^2|m \rangle -\sum _{k=1}^3 B_k(t)B_r(t)\langle n|q_kq_r|m \rangle \bigg \}\\&\quad =\bigg \{(B_2^2(t)+B_3^2(t)) \langle n|q_1^2|m \rangle + (B_3^2(t)+B_1^2(t)) \langle n|q_2^2|m \rangle \\&\qquad +\, (B_1^2(t)+B_2^2(t)) \langle n|q_3^2|m \rangle -2\sum _{1<k<r<3} B_k(t)B_r(t)\langle n|q_kq_r|m \rangle \bigg \}\text {e}^{\iota t E(n,m)t} \end{aligned}$$
To calculate \(\langle n|W_2(t)|m \rangle \), we also need \(\langle n|L_k\big (t_1\big )L_r\big (t_2\big )|m \rangle \) and \(\langle n|L_k\big (t_1\big )q_r\big (t_2\big )|m \rangle \)
$$\begin{aligned} L_k\big (t_1\big )L_r\big (t_2\big )=\text {e}^{\iota t_1H_0}L_kexp(-\iota (t_1-t_2)H_0)L_r\mathrm{e}^{-\iota t_2H_0} \end{aligned}$$
So ,
$$\begin{aligned} \langle n|L_k\big (t_1\big )L_r\big (t_2\big )|m \rangle =\text {e}^{\iota (E(n)t_1-E(m)t_2)}\sum _s \langle n|L_k|s \rangle \langle s|L_r|m \rangle \mathrm{e}^{-\iota (t_1-t_2)E(s)} \end{aligned}$$
The other quantity required is
$$\begin{aligned} T_1&= \langle n|\int _{0<t_1<t_2<t} B\big (t_1\big )\cdot L\big (t_1\big )B\big (t_2\big )\cdot L\big (t_2\big ){\text {d}}t_1{\text {d}}t_2|m \rangle \\&=\sum _{k,r=1}^3 \int _{0<t_1<t_2<t} B_k\big (t_1\big )\cdot B_r\big (t_2\big )\langle n|L_k\big (t_1\big )L_r\big (t_2\big )|m \rangle {\text {d}}t_1{\text {d}}t_2\\&=\sum _{k,r,s} \int _{0<t_1<t_2<t} B_k\big (t_1\big )\cdot B_r\big (t_2\big )\mathrm{e}^{-\iota (t_1-t_2)E(s)} \langle n|L_k|s \rangle \langle s|L_r|m \rangle \text {e}^{\iota (E(n)t_1-E(m)t_2}{\text {d}}t_1{\text {d}}t_2 \end{aligned}$$
We now substitute for \(\langle n|L_r|m \rangle \) in the above equation but firstly, let us define \(e_r=(\delta _{r1},\delta _{r2},\delta _{r3})\), \(r=1,2,3\) ,that is,
$$\begin{aligned} e_1=(100),e_2=(010),e_3=(001) \end{aligned}$$
As can be seen clearly,
$$\begin{aligned} {\mathcal {Z}}^{-1}_r{\mathcal {Z}}_s \delta [n-m]=\delta [n-m-e_r+e_s] \end{aligned}$$
So, on incorporating all the above we have \(T_1\) as
$$\begin{aligned} T_1&=-\sum _{k,r,s,\alpha ,\beta ,\mu ,\nu }\bigg (\int _{0<t_2<t_1<t}B_k\big (t_1\big )\cdot B_r\big (t_2\big )\mathrm{e}^{-\iota (t_1-t_2)E(s)} \epsilon (k\alpha \beta )\epsilon (r\mu \nu )\\&\quad \sqrt{n_\alpha s_\beta } \sqrt{s_\mu m_\nu }\delta [n-s-e_\alpha +e_\beta ] \delta \big [s-m-e_\mu +e_\nu \big ]\text {e}^{\iota (E(n)t_1-E(m)t_2)}\bigg ){\text {d}}t_1{\text {d}}t_2\\&=-\sum _{k,r,\alpha ,\beta ,\mu ,\nu }\bigg (\int _{0<t_2<t_1<t}B_k\big (t_1\big )\cdot B_r\big (t_2\big )\mathrm{e}^{-\iota (t_1-t_2)E(n-e_\alpha +e_\beta )} \epsilon (k\alpha \beta )\epsilon (r\mu \nu )\\&\quad \times \, \sqrt{n_\alpha \big (n_\beta -\delta _{\alpha \beta }+1\big ) m_\nu \big (m_\mu -1+\delta _{\mu \nu }\big )}\delta \big [n-m-e_\alpha +e_\beta +e_\mu -e_\nu \big ]\\&\quad \times \, \text {e}^{\iota \big (E(n)t_1-E(m)t_2\big )}\bigg ){\text {d}}t_1{\text {d}}t_2 \end{aligned}$$
Defining
$$\begin{aligned} K_{k,r}^{(n,m)}\big (t_1,t_2\big )&=\sum _s \mathrm{e}^{-\iota (t_1-t_2)E(s)}\langle n|L_k|s \rangle \langle s|L_r|m \rangle \mathrm{e}^{\iota \big (E(n)t_1-E(m)t_2\big )}\\&=-\sum _{\alpha \beta \mu \nu } \mathrm{e}^{-\iota \big (t_1-t_2\big ) E\big (n-e_\alpha +e_\beta \big )}\epsilon (k\alpha \beta )\epsilon (r\mu \nu )\\&\quad \times \,\sqrt{n_\alpha \big (n_\beta -\delta _{\alpha \beta }+1\big ) m_\nu \big (m_\mu -1+\delta _{\mu \nu }\big )}\delta \left[ n-m-e_\alpha +e_\beta +e_\mu -e_\nu \right] \\&\quad \times \, \text {e}^{\iota \big (E(n)t_1-E(m)t_2\big )} \end{aligned}$$
Thus,
$$\begin{aligned} \langle n|W_2(t)|m \rangle&=-\frac{i}{8}\int _0^t\langle n| (B\big (t_1\big )\times q\big (t_1\big ))^2|m\rangle {\text {d}}t_1\\&\quad -\frac{1}{8}\sum _{k,r}\int _{0<t_2<t_1<t}K_{k,r}^{(n,m)}\big (t_1,t_2\big )B_k\big (t_1\big )B_r\big (t_2\big ){\text {d}}t_1{\text {d}}t _2\\&\quad -\,\frac{1}{4}\int _{0<t_2<t_1<t}\langle n|B\big (t_1\big )\cdot L\big (t_1\big )E\big (t_2\big )\cdot q\big (t_2\big )|m \rangle {\text {d}}t_1{\text {d}}t_2\\&\quad -\,\frac{1}{4}\int _{0<t_2<t_1<t}\langle n|E\big (t_1\big )\cdot q\big (t_1\big )B\big (t_2\big )\cdot L\big (t_2\big )|m \rangle {\text {d}}t_1{\text {d}}t_2\\&\quad -\,\frac{1}{2}\int _{0<t_2<t_1<t}\langle n|E\big (t_1\big )\cdot q\big (t_1\big )E\big (t_2\big )\cdot q\big (t_2\big )|m \rangle {\text {d}}t_1{\text {d}}t_2 \end{aligned}$$
We need to calculate the general matrix element of \((B(t)\times q(t))^2\), so
$$\begin{aligned} (B(t)\times q(t))^2= & {} \sum _{k\alpha \beta \mu \nu }\{\epsilon (k\alpha \beta )\epsilon (k\mu \nu )B_\alpha (t)B\mu (t) q_\beta (t)q_\nu (t)\}\\ \langle n|(B(t)\times q(t))^2|m \rangle= & {} \sum _{k\alpha \beta \mu \nu }\{\epsilon (k\alpha \beta )\epsilon (k\mu \nu )B_\alpha (t)B\mu (t) \langle n|q_\beta (t)q_\nu (t)|m \rangle \} \end{aligned}$$
Define
$$\begin{aligned} G_{\alpha \mu }^{(n,m)}(t)=\sum _{k\beta \nu }\{\epsilon (k\alpha \beta )\epsilon (k\mu \nu ) \langle n|q_\beta q_\nu |m \rangle \}\mathrm{e}^{itE(n,m)} \end{aligned}$$
Thus,
$$\begin{aligned} T_2&\triangleq \int _o^t \langle n|(B(t)\times q(t))^2|m \rangle {\text {d}}t_1=\sum _{k,r} \int _o^t G_{kr}^{(n,m)}\big (t_1\big )B_k\big (t_1\big )B_r\big (t_1\big ){\text {d}}t_1\\ T_3&\triangleq \int _{0<t_2<t_1<t}\langle n|B\big (t_1\big )\cdot L\big (t_1\big )E\big (t_2\big )\cdot q\big (t_2\big )|m \rangle {\text {d}}t_1{\text {d}}t_2 \end{aligned}$$
where
$$\begin{aligned} \langle n|B\big (t_1\big )\cdot L\big (t_1\big )E\big (t_2\big )\cdot q\big (t_2\big )|m \rangle =\sum _{k,r}B_k\big (t_1\big )E_r\big (t_2\big )\langle n|L_k\big (t_1\big )q_r\big (t_2\big )|m \rangle \end{aligned}$$
and
$$\begin{aligned} \langle n|L_k\big (t_1\big )q_r\big (t_2\big )|m \rangle&= \langle n|\mathrm{e}^{\iota t_1H_0}L_k\mathrm{e}^{-\iota (t_1-t_2)H_0}q_r\mathrm{e}^{-\iota t_2H_0}|m \rangle \\&=\text {e}^{\iota (E(n)t_1-E(m)t_2)}\sum _s \langle n|L_k|s \rangle \langle s|q_r|m \rangle \mathrm{e}^{-\iota (t_1-t_2)}E(s)) \end{aligned}$$
Let the above equation be equal to \(H_{k,r}^{(n,m)}\big (t_1,t_2\big )\). Thus,
$$\begin{aligned} T_3= & {} \sum _{k,r}\int _{o<t_2<t_1<t}H_{k,r}^{(n,m)}\big (t_1,t_2\big )B_k\big (t_1\big )E_r\big (t_2\big ){\text {d}}t_1{\text {d}}t_2\\ T_4\triangleq & {} \int _{o<t_2<t_1<t} \langle n|E\big (t_1\big )\cdot q\big (t_1\big )B\big (t_2\big )\cdot L\big (t_2\big )|m \rangle {\text {d}}t_1{\text {d}}t_2\\= & {} \sum _{k,r} \int _{o<t_2<t_1<t} E_k\big (t_1\big )B_r\big (t_2\big ) \langle n|q_k\big (t_1\big )L_r\big (t_2\big )|m \rangle {\text {d}}t_1{\text {d}}t_2\\ \langle n|q_k\big (t_1\big )L_r\big (t_2\big )|m \rangle= & {} \text {e}^{\iota (E(n)t_1-E(m)t_2)} \langle n|q_k\mathrm{e}^{-\iota (t_1-t_2)H_0}L_r|m \rangle \\= & {} \text {e}^{\iota (E(n)t_1-E(m)t_2)} \sum _s \langle n|q_k|s \rangle \langle s|L_r|m \rangle \mathrm{e}^{-\iota (t_1-t_2)E(s)}\\\triangleq & {} L_{k,r}^{(n,m)}\big (t_1,t_2\big ) \end{aligned}$$
So \(T_4\) finally becomes
$$\begin{aligned} T_4=\sum _{k,r} \int _o^t L_{kr}^{(n,m)}\big (t_1,t_2\big )E_k\big (t_1\big )B_r\big (t_2\big ){\text {d}}t_1{\text {d}}t_2 \end{aligned}$$
Similarly,
$$\begin{aligned} T_5&\triangleq \int _{o<t_2<t_1<t}\langle n|E\big (t_1\big )\cdot q\big (t_1\big )E\big (t_2\big )\cdot q\big (t_2\big )|m \rangle {\text {d}}t_1{\text {d}}t_2\\&=\sum _{k,r} \int _{o<t_2<t_1<t} E_k\big (t_1\big )E_r\big (t_2\big ) \langle n|q_k\big (t_1\big )q_r\big (t_2\big )|m \rangle {\text {d}}t_1{\text {d}}t_2 \end{aligned}$$
we substitute the following in the above equation
$$\begin{aligned} \langle n|q_k\big (t_1\big )q_r\big (t_2\big )|m \rangle&=\text {e}^{\iota (E(n)t_1-E(m)t_2)} \sum _s \langle n|q_k|s \rangle \langle s|q_r|m \rangle \mathrm{e}^{-\iota (t_1-t_2)E(s)}\\&\triangleq \mu _{k,r}^{(n,m)}\big (t_1,t_2\big ) \end{aligned}$$
So \(T_5\) finally becomes
$$\begin{aligned} T_5=\sum _{k,r} \int _o^t \mu _{kr}^{(n,m)}\big (t_1,t_2\big )E_k\big (t_1\big )E_r\big (t_2\big ){\text {d}}t_1{\text {d}}t_2 \end{aligned}$$
Thus, we have \(W_2(t)\) as
$$\begin{aligned} \langle n|W_2(t)|m \rangle&=-\frac{i}{2}\sum _{k,r} \int _o^t G_{kr}^{(n,m)}\big (t_1\big )B_k\big (t_1\big )B_r\big (t_1\big ){\text {d}}t_1\\&\quad -\,\frac{1}{8}\sum _{k,r} \int _o^t K_{kr}^{(n,m)}\big (t_1,t_2\big )B_k\big (t_1\big )B_r\big (t_2\big ){\text {d}}t_1{\text {d}}t_2\\&\quad -\,\frac{1}{4}\sum _{k,r} \int _o^t H_{kr}^{(n,m)}\big (t_1,t_2\big )B_k\big (t_1\big )E_r\big (t_2\big ){\text {d}}t_1{\text {d}}t_2\\&\quad -\,\frac{1}{4}\sum _{k,r} \int _o^t L_{kr}^{(n,m)}\big (t_1,t_2\big )E_k\big (t_1\big )B_r\big (t_2\big ){\text {d}}t_1{\text {d}}t_2\\&\quad -\,\frac{1}{2}\sum _{k,r} \int _o^t \mu _{kr}^{(n,m)}\big (t_1,t_2\big )E_k\big (t_1\big )E_r\big (t_2\big ){\text {d}}t_1{\text {d}}t_2 \end{aligned}$$
Define the kernels
$$\begin{aligned} g_{11}\big (t_1,t_2;n,m,k,r\big )= & {} -\frac{i}{2} G_{kr}^{(n,m)}\big (t_1,t_2\big )\big (t_1\big )\delta \big (t_1-t_2\big )-\frac{1}{8} K_{kr}^{(n,m)}\theta \big (t_1-t_2\big )\\ g_{12}\big (t_1,t_2;n,m,k,r\big )= & {} -\frac{1}{4} H_{kr}^{(n,m)}\big (t_1,t_2\big )\big (t_1\big )\theta \big (t_1-t_2\big )-\frac{1}{4} L_{kr}^{(n,m)}\theta \big (t_2-t_1\big )\\ g_{22}\big (t_1,t_2;n,m,k,r\big )= & {} -\frac{1}{2} \mu _{kr}^{(n,m)}\big (t_1,t_2\big )\theta \big (t_1-t_2\big ) \end{aligned}$$
The general matrix element of \(W_2(t)\) becomes
$$\begin{aligned} \langle n|W_2(t)|m \rangle&=\sum _{k,r} \int _{[0,t]^2} g_{11}\big (t_1,t_2;n,m,k,r\big )B_k\big (t_1\big )B_r\big (t_2\big ){\text {d}}t_1{\text {d}}t_2\nonumber \\&\quad +\,\sum _{k,r} \int _{[0,t]^2} g_{12}\big (t_1,t_2;n,m,k,r\big )B_k\big (t_1\big )E_r\big (t_2\big ){\text {d}}t_1{\text {d}}t_2\nonumber \\&\quad +\,\sum _{k,r} \int _{[0,t]^2} g_{22}\big (t_1,t_2;n,m,k,r\big )E_k\big (t_1\big )E_r\big (t_2\big ){\text {d}}t_1{\text {d}}t_2 \end{aligned}$$
(45)
Likewise,
$$\begin{aligned} \langle n|W_1(t)|m \rangle&= -\frac{i}{2}\sum _{k} \int _o^t B_k\big (t_1\big ) \langle n|L_k\big (t_1\big )|m \rangle {\text {d}}t_1-i \sum _{k} \int _o^t E_k\big (t_1\big )\langle n|q_k\big (t_1\big )|m \rangle {\text {d}}t_1\\ \langle n|L_k(t)|m \rangle&=\mathrm{e}^{\iota E(n,m)t} \langle n|L_k|m \rangle \triangleq f(t,n,m,k)\\ \langle n|q_k(t)|m \rangle&=\mathrm{e}^{\iota E(n,m)t} \langle n|q_k|m \rangle \triangleq g(t,n,m,k) \end{aligned}$$
Then, define
$$\begin{aligned} h_1(t,n,m,k)= & {} -\frac{i}{2}f(t,n,m,k)\\ h_2(t,n,m,k)= & {} -i g(t,n,m,k) \end{aligned}$$