# Methods for integration of oscillatory complex vectors as a function of time

I'm attempting to solve a problem of the form:

$$\mathbf{a}^{(n+1)}(t) = \int_{0}^{t}d\tau e^{i\mathbf{H}\tau} \mathbf{D}(\tau)e^{-i\mathbf{H}\tau}\mathbf{a}^{(n)}(\tau)$$ Where $\mathbf{D}(\tau)$ is a oscillatory time varying matrix ($\mathbf{D}\cdot E(\tau)$), $\mathbf{H}$ is a diagonal matrix, and $\mathbf{a}^{(n)}(t)$ is a vector function of time. I would like to obtain $\mathbf{a}^{(n+1)}(t)$. $\mathbf{D}$ oscillates with frequency significantly less than the oscillations due to the exponents, in general, though later iterations will change that. The frequencies for different terms of the vector are likely to be very different (by a factor of $10$ is likely). In later iterations $\mathbf{a}^{(n)}(t)$ will also oscillate.

Likely important to note: the first step in this ($a^{(0)}$) is constant, and thus the integrand of the first part is known for all $t$.

As you might be able to guess by the $(n+1)$ and $(n)$, this is an iterative process, and I would like to know both $\mathbf{a}^{(n)}(t)$ and $\mathbf{a}^{(n+1)}(t)$ for a number of evenly spaced points per period of $\mathbf{D}(t)$.

I'm not sure how I should be going about doing this. Since I would like to find the integral for a number of different values of $t$, something like a Riemann sum (like the trapezoid rule or similar) seems like a better idea than something like Gaussian quadrature, however, I'm not sure what is a good method when dealing with complex oscillatory functions. Is an adaptive method a good option? If so, what are good methods for finding the adaptive step for vectors?

The vector nature of this is only moderately important: the vectors are approximately 3000 values long, so I would rather not do a different calculation for each vector, instead preferring to find them all at once.

• Should the $D(t)$ under the integrand be $D(\tau)$, or is it correct like this? – davidhigh Apr 28 '17 at 0:34
• @davidhigh : good catch, I've fixed it. – Andrew Spott Apr 28 '17 at 0:36

That is a Volterra integral equation in principle, and there are several methods to solve it.

I would first differentiate the whole expression to obtain

$$\tag 1 \partial_t \mathbf{a}^{(n+1)}(t) = e^{i\mathbf{H}t} \mathbf{D}(t)e^{-i\mathbf{H}t}\mathbf{a}^{(n)}(\tau)$$

Define $\mathbf b^{(i)}(t) = e^{-i\mathbf{H}t} \mathbf a^{(i)}$, then $$\tag 2 i\partial_t \mathbf{b}^{(n+1)}(t) = \mathbf{H}(t) \mathbf{b}^{(n)}(t) + \mathbf{D}(t) \mathbf{b}^{(n)}(t),\\ \mathbf{b}^{(n+1)}(0) = 0\,.$$

Ok, this might look to you like going back, as this is probably where you came from: Schrödinger equation, electric field in dipole approximation, interaction picture, and so on -- you know the stuff. But the point is, you can solve that latter differential equation much better than the integral equation above.

So drop the upper iteration index, $$\tag 3 i\partial_t \mathbf{b}(t) = \Big(\mathbf{H}(t) + E(t)\mathbf{D} \Big)\mathbf{b}(t)$$ and just apply some basic Runge-Kutta for the first. Thereafter, I'd look at Crank-Nicolson, Lanczos schemes, and so on.

EDIT: ok, you want the iteration summands by itself. For that, one can use the boundary condition $\mathbf b^{(0)}(t) = 1$ (the first term in the Dyson series) and starting from this again solve the differential equation.

The underlying idea here is that the indexed version, Eq.(2), is some kind of fix-point problem (starting with a more or less irrelevant initial value), which, when it hopefully converges, also satisfies the Schrödinger equation, Eq.(3).

I know, this is probably again not what you've expected, but your integral and Eq.(2) are mathematically equivalent, and the latter is simply easier to solve, so I'd use it.

• This is definitely where I came from. :). Unfortunately, what I want is the iterative versions: I want $\mathbf{a}^{(n)}(t)$, not the sum of all of them... I already have a solver for the differential equation... :/ – Andrew Spott Apr 28 '17 at 1:08
• @AndrewSpott: out of curiosity: why? Otherwise, see my edit. – davidhigh Apr 28 '17 at 1:25
• Because I'm trying to look at the difference between the perturbative and the non-perturbative problems. In regards to your edit, I'm looking at a problem of the form y'(t) = f(t), so Runge-kutta, Crank-Nicolson or similar? It always feels a little weird using those to solve what is essentially an integral... – Andrew Spott Apr 28 '17 at 1:33
• So, Crank-Nicolson reduces to the trapezoid method for $f$ not a function of $y$, and Runge-Kutta reduces to Simpsons rule when $f$ isn't a function of $y$. Are these methods good enough for this kind of oscillatory integral? – Andrew Spott Apr 28 '17 at 2:01
• @AndrewSpott: even if it's weird from the first, note that time-propagation is also often called "integration". And your stated equivalences hold at least intuitively. For explicit methods (and under special circumstances) it's clearer: for example, as you said, the 4th-order Runge-Kutta scheme becomes the Simpson rule. Not quite sure at the moment about Crank-Nicolson. – davidhigh Apr 28 '17 at 9:47