An efficient spectral method for numerical time-dependent perturbation theory

TL;DR

The paper introduces FLIPT, a fast and exact spectral method for calculating time-dependent perturbation expansions of the density matrix, optimized for complex laser pulse simulations.

Contribution

It presents FLIPT, a novel spectral approach that efficiently computes perturbative expansions with rigorous convergence, reducing computational complexity for high-dimensional integrals.

Findings

FLIPT achieves $O(N^2)$ complexity for n-dimensional integrals.

The method is highly suitable for simulating complex multiphoton laser experiments.

FLIPT provides numerically exact results with rigorous convergence guarantees.

Abstract

We develop the Fourier-Laplace Inversion of the Perturbation Theory (FLIPT), a novel numerically exact "black box" method to compute perturbative expansions of the density matrix with rigorous convergence conditions. Specifically, the FLIPT method is extremely well-suited to simulate multiphoton pulsed laser experiments with complex pulse shapes. The $n$-dimensional frequency integrals of the $n$-th order perturbative expansion are evaluated numerically using tensor products. The $N$ points discretized integrals are computed in $O(N^2)$ operations, a significant improvement over the $O(N^n)$ scaling of standard quadrature methods.

Tables

Table 1. Table 1 : Parameters used for the 2D models of pyrazine obtained from Ref.  Sukharev and Seideman, 2004 . All values have units of eV.

${\omega }_t$	${\omega }_c$	${\kappa }_1$	${\kappa }_2$	${ϵ}_1$	${ϵ}_2$	$\lambda$
0.0739	0.1139	-0.0981	0.1355	3.94	4.89	0.1830

Equations

$$H(t)=H_0+\lambda E(t)V,$$

$$\lambda E(t)V=-\sum _{\alpha }{\epsilon }_{\alpha }(\mathbf{r},t){\mu }_{\alpha },$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\rho (t)$$

$${\mathcal{L}}_0\rho$$

$$\mathcal{V}\rho$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\rho (t)$$

$$(z-{\mathcal{L}}_0)\rho (z)$$

$$\rho (z)={\int }_0^{\infty }\mathrm{d}t{e}^{-zt}\rho (t),$$

$$\rho (z)$$

$$+{\left(\frac{\lambda }{2\pi }\right)}^2{\iint }_{-\infty }^{\infty }\mathrm{d}{\omega }^{\prime }\mathrm{d}{\omega }^{\prime \prime }E({\omega }^{\prime \prime })E({\omega }^{\prime }){\mathcal{G}}_0(z)\mathcal{V}{\mathcal{G}}_0(z-i{\omega }^{\prime })\mathcal{V}{\mathcal{G}}_0(z-i{\omega }^{\prime }-i{\omega }^{\prime \prime }){\rho }_0+O({\lambda }^3)$$

$$={\rho }_0(z)+{\rho }_1(z)+{\rho }_2(z)+\cdots$$

$${\rho }_0(i\omega +\eta )$$

$${\rho }_n(i\omega +\eta )$$

$${\rho }_1(i\omega +\eta )$$

$$=\frac{\lambda }{2\pi }{\mathcal{G}}_0(i\omega +\eta )\mathcal{V}{\rho }_0{\int }_{-\infty }^{\infty }\mathrm{d}{\omega }^{\prime }\frac{E({\omega }^{\prime })}{i\omega -i{\omega }^{\prime }+\eta }.$$

$${\int }_{-\infty }^{\infty }\mathrm{d}t{e}^{-i\omega t-\eta t}\Theta (t)=\frac{1}{i\omega +\eta }.$$

$${\int }_{-\infty }^{\infty }\mathrm{d}{\omega }^{\prime }\frac{E({\omega }^{\prime })}{i\omega -i{\omega }^{\prime }+\eta }$$

$$=2\pi {\int }_{-\infty }^{\infty }\mathrm{d}t{e}^{-i\omega t-\eta t}\Theta (t)E(t).$$

$$E_{\eta }(t)=E(t)e^{-\eta t}.$$

$${\rho }_1(i\omega +\eta )$$

$${\int }_{-\infty }^{\infty }\mathrm{d}\Omega \delta (\Omega -\omega ),$$

$${\rho }_n(i\omega +\eta )$$

$$\times {\mathcal{G}}_0(i{\omega }_n+\cdots +i{\omega }_1+\eta )\mathcal{V}{\mathcal{G}}_0(i{\omega }_{n-1}+\cdots +i{\omega }_1+\eta )\cdots$$

$$\times \mathcal{V}{\mathcal{G}}_0(i{\omega }_1+\eta )\mathcal{V}{\rho }_0.$$

$${\rho }_n(t)=\frac{1}{2\pi i}{\int }_{\eta -i\infty }^{\eta +i\infty }\mathrm{d}z{e}^{zt}{\rho }_n(z)=\frac{e^{\eta t}}{2\pi }{\int }_{-\infty }^{\infty }\mathrm{d}\omega e^{i\omega t}{\rho }_n(i\omega +\eta ),$$

$${\rho }_n(t)$$

$$\times {\mathcal{G}}_0(i{\omega }_n+\cdots +i{\omega }_1+\eta )\mathcal{V}{\mathcal{G}}_0(i{\omega }_{n-1}+\cdots +i{\omega }_1+\eta )\cdots \mathcal{V}{\mathcal{G}}_0(i{\omega }_1+\eta )\mathcal{V}{\rho }_0.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}|\psi (t)⟩$$

$${\mathcal{L}}_0$$

$$\mathcal{V}$$

$${\rho }_n(t)$$

$$|{\psi }_n(t)⟩$$

$$\times E({\omega }_n)\cdots E({\omega }_2)E_{\eta }({\omega }_1-{ϵ}_0/\hslash )$$

$$\times G_0(i{\omega }_n+\cdots +i{\omega }_1+\eta )V{G}_0(i{\omega }_{n-1}+\cdots +i{\omega }_1+\eta )\cdots V{G}_0(i{\omega }_1+\eta )V|{\psi }_0⟩,$$

$$G_0(i\omega +\eta )$$

$$V$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\rho (t)$$

Full text

I Introduction

At the microscopic level, the interaction between semiclassical light and matter is well-described by perturbation theory.Shapiro and Brumer (2012) Indeed, optical processes such as two-photon absorption and Raman scattering are often described and classified in terms of discrete interactions with radiation, a picture based on perturbation theory.Cohen-Tannoudji, Dupont-Roc, and Grynberg (1992) Similarly, analytical results from perturbation theory are widely used to analyse nonlinear and ultrafast spectroscopic experiments,Gallagher Faeder and Jonas (1999); Zhuang et al. (2006); Quesada and Sipe (2014); Brühl, Buckup, and Motzkus (2018); Reppert and Brumer (2018a) the light-induced control of quantum systems,Pachón and Brumer (2013); Mukamel (2013); Am-Shallem and Kosloff (2014); Lavigne and Brumer (2017, 2019a) and the dynamics of photoactivated natural processes.Tscherbul and Brumer (2014, 2015); Brumer (2018); Reppert and Brumer (2018b)

Although perturbative analysis underlies much of the theory of optical processes, it is not generally used in numerical simulations without additional approximations. For example, pulsed electric fields are often approximated as an infinitely short $\delta(t)$ pulse,Smallwood, Autry, and Cundiff (2017); Perlík, Hauer, and Šanda (2017) and nearly monochromatic electric fields as infinitely long continuous wave (CW) oscillations, i.e. $\delta$ functions in frequency. Both approximations suffer from important drawbacks. It is well-known that the phase of an electric field affects the dynamics of matter interacting with said field.Shapiro and Brumer (2012) Such phase effects range from the trivial reduction in two-photon absorption probability when chirping ultrashort pulsesLavigne and Brumer (2019a) to the non-trivial quantum control of molecular dynamics.Brumer and Shapiro (1992); Mukamel (2013); Am-Shallem and Kosloff (2014); Lavigne and Brumer (2017); Brühl, Buckup, and Motzkus (2018); Lavigne and Brumer (2019a) However, these effects are lost in the limit of infinitely short pulses. Conversely, simulations of interaction with CW radiation are of limited applicability, especially in the presence of other time-dependent processes such as decoherence and dissipation.Spanner, Arango, and Brumer (2010) Furthermore, the resulting time-independent equations do not converge without renormalizationFaisal (1987) or the inclusion of broadening factors.Mukamel (1995); Plenio, Almeida, and Huelga (2013)

In cases where an experimentally realistic electric field is required, e.g., in the interpretation of control experiments,Lavigne and Brumer (2017) the relevant Schrodinger equation can be solved by means of a non-perturbative propagation of the time-dependent Schrodinger equation, with the radiation included in the Hamiltonian as a time-dependent potential.Katz, Ratner, and Kosloff (2010); Arango and Brumer (2013); Abe et al. (2005); Brühl, Buckup, and Motzkus (2018) A “one-photon” or “two-photon” result would be obtained in such a calculation by choosing a small intensity for the electric field.

Time-dependent propagation methods, unlike perturbative approaches, do not suffer from the lack of phase effects or convergence difficulties described above. However, direct propagation is non-perturbative, which has important drawbacks when computing weak-field processes. First, from a numerical analysis point of view, a high accuracy is required to adequately simulate the small population excited by the field.Arango and Brumer (2013) In addition, short timesteps (sub-fs in the case of visible radiation) are required to capture the fast oscillations of electric fields at optical frequencies. The rotating wave approximation (RWA) greatly speeds up convergence by allowing for much longer timesteps. However, the RWA breaks down for non-resonant processesMilonni and Eberly (2010) and multiphoton absorption.Faisal (1987)

Second, non-perturbative results are difficult to interpret and to compare with experimental results that are described in the framework of perturbation theory, e.g., described by $n$-wave mixing or $n$ photon absorption processes. For example, higher order processes can appear in “linear regime” experiments and simulations thereof. Han and Shapiro (2012, 2013); Brühl, Buckup, and Motzkus (2018); Lavigne and Brumer (2019b) In fact, non-perturbative microscopic simulations (i.e. where only the wavefunction of a single molecule is evolved) suffer in particular from uncontrolled contributions of unwanted optical signals as they are not subject to the phase matching conditions arising from the macroscopic Maxwell equations. The final computed result is then a sum of signals from processes that can only be unmixed with great difficulty in simulations.Meyer and Engel (2000) In contrast, the corresponding experimental signals are spatially separated.

In this article, a numerical algorithm is introduced to evaluate arbitrary order perturbative expansions of the density matrix using the Fourier-Laplace inversion.Dubner and Abate (1968); Veillon (1974); Crump (1976); de Hoog, Knight, and Stokes (1982); Piessens and Huysmans (1984) The Laplace-transformed time-dependent perturbation theory is shown to be convergent for experimentally relevant light-matter processes. Perturbative contributions are obtained as iterated frequency integrals. A Fourier series discretization is used to evaluate the perturbative integrals exponentially faster than standard quadrature using a tensor product technique.Gerstner and Griebel (1998) The resultant Fourier-Laplace Inversion of the Perturbation Theory (FLIPT) method can be expressed succinctly using tensor notation and the associated tensor product algebra is easily implemented using multidimensional arrays. The resultant, highly efficient implementation is made freely available by the authors.Lavigne (2019)

The FLIPT algorithm introduced here can be used as a “black box” method to simulate light-matter interactions from coherent,*[Incoherentradiation; suchassunlight; isbettertreatedusingothermethodssuchasthatdescribedin][]axelrod_efficient_2018 pulsed laser fields of the type used in ultrafast spectroscopy. The implementation is fully automatic and contains no free parameters. Since the computation is performed in the frequency domain and not in the time domain, fast and slow dynamical observables are equally resolved. This property is particularly useful for chemical processes where coherent excitation dynamics on the fs timescale lead to ps to ns reaction dynamics. Both are obtained in the FLIPT method in a numerically exact manner, with uniform convergence and without the RWA.

Calculations below and in Ref. Lavigne and Brumer, 2019a show that the FLIPT method can be applied to moderately large multilevel systems (N $\approx$ 300-600) with timescales ranging from the sub-fs dynamics of electronic coherences to the ps dynamics of molecular vibrations. Furthermore, the method is not limited to specific form for the exciting light and can equally treat resonant and non-resonant multiphoton processes. Extensions of the algorithm to the perturbative analysis of wavefunction and non-Markovian dynamics are discussed. Convergence is rigorously established for time-limited excitations of varying duration and is independent of the underlying system dynamics.Dubner and Abate (1968); Crump (1976)

II Theory

The theory underlying the efficient numerical method alluded to above is developed here. The quantity to be computed is the perturbative expansion of the density matrix under the action of time-dependent potentials of the type relevant in ultrafast laser experiments. The usual time-dependent perturbative expansion of the density matrix is re-derived to show that the Laplace transform formalism used here is convergent and well-behaved. Below, the perturbative expansion and its computation through a Laplace inversion is demonstrated for the density matrix obeying the Liouville-von Neumann equation. Possible extensions to wavefunction dynamics and open system dynamics are briefly described.

The Hamiltonian of interest is given by,

[TABLE]

where $H_{0}$ is the zeroth-order Hamiltonian, $\lambda$ is a small dimensionless scalar and $V$ is the coupling operator for the perturbation that evolves under $E(t)$, a scalar function of time.111The algorithm can easily be extended to the case where the perturbation is composed of multiple components $\sum_{\alpha}E_{\alpha}(t)V_{\alpha}(t)$ by computing and summing over all unique combinations of $\alpha$. Significantly, such an Hamiltonian describes the interaction of a molecule with a time-varying classical electric field in the dipole approximation,Shapiro and Brumer (2012) with

[TABLE]

where $\bm{r}$ is the position of the molecule, $\varepsilon_{\alpha}(\bm{r},t)$ is the $\alpha=x,y,z$ component of the electric field and $\mu_{\alpha}$ is the $\alpha$ component of the dipole transition operator for the molecule.

The FLIPT method depends crucially on the convergence of the Laplace transform applied to the perturbation, which is guaranteed in the experimentally relevant case of a perturbation of finite duration, as shown in this section. Specifically, the perturbation $E(t)$ is taken to be bounded and time-limited, i.e., $E(t)=0$ for all $t$ less than some “turn-on time” $t_{\text{on}}$ and larger than some “turn-off time” $t_{\text{off}}$. This condition guarantees convergence, both analytically and numerically, as shown below.

II.1 Perturbative expansion and Laplace inversion

The Liouville-von Neumann equation of motion for the density matrix in the superoperator formalism is given by,Lendi (1977)

[TABLE]

where Liouvillian and coupling superoperators are so defined. A Fourier inversion of $E(t)$ yields the following,

[TABLE]

Without loss of generality, the origin $t=0$ is chosen before the “turn-on time” $t_{\text{on}}$. A Laplace transform yields,Boas (2005)

[TABLE]

where

[TABLE]

and $\rho_{0}$ is the initial state $\rho(t=0)$. Multiplying through by the Green’s function $\mathcal{G}_{0}(z)=(z-\mathcal{L}_{0})^{-1}$, an implicit integral equation for $\rho(z)$ is obtained,Löwdin (1982)

[TABLE]

The perturbative expansion results from an iteration over $\rho(z)$:

[TABLE]

The Laplace-transformed $\rho_{n}(z)$ is the $n$-th order perturbative contribution to $\rho(z)$. The iteration procedure can be written explicitly and succinctly as,

[TABLE]

Thus, $\rho_{n}(i\omega+\eta)$ is obtained from a convolution of $E(\omega)\mathcal{V}$ with $\rho_{n-1}(i\omega-i\omega^{\prime}+\eta)$ followed by an application of $\mathcal{G}_{0}(i\omega+\eta)$.

An important special case is where the system is initially in a steady state, such that $\mathcal{L}_{0}\rho_{0}=0$. Then, the first order contribution becomes,

[TABLE]

The integral over $\omega^{\prime}$ is analytically solvable. The one-sided Fourier transform of an exponentially decaying function is given by,

[TABLE]

where $\Theta(t)$ is the Heaviside step function. The integral in eq. (15) is thus given by,

[TABLE]

Since the perturbation $E(t)$ is zero for $t<0$, eq. (18) is the Fourier transform of the following function,

[TABLE]

Therefore, the first order term in eq. (13) becomes,

[TABLE]

A convenient form for the $n$-th order term can be obtained by introducing an additional frequency integration,

[TABLE]

and performing the change of variables $\Omega=\omega_{1}+\omega_{2}\cdots$, where $\omega_{i}$ is the frequency variable for the $i$-th perturbation. The $n$-th perturbative contribution is then given by a $n$-dimensional frequency integral, obtained from $n$ perturbations $E(\omega_{1}),E(\omega_{2})\cdots$ of the initial state $\rho_{0}$,

[TABLE]

The case of a system initially in a steady state is particularly important, as that describes most spectroscopic and control experiments.Lavigne and Brumer (2019b) It is also significantly simpler as it yields time-translationally invariant dynamics, that is dynamics that do not depend on the absolute value of the initial time $t_{\text{on}}$.222 The case where the initial state is not a steady state can be computed directly from eq. (12) and (13) above. This is numerically more expensive as it lead to an $n+1$ dimensional integral for the $n$-th perturbative term instead of an $n$ dimensional integral.

The time-dependent $n$-th order contribution $\rho_{n}(t)$ can be obtained by the Laplace inversion integral,

[TABLE]

where $z=i\omega+\eta$ and $\eta$ is a real number greater than the real part of all the poles of $\rho_{n}(z)$. Substituting eq. (22) into eq. (23) yields the following $n$ multidimensional inverse Fourier transform solution for $\rho_{n}(t)$,

[TABLE]

For a closed multilevel system, all eigenvalues of the Liouvillian superoperator $\mathcal{L}_{0}$ have a zero real part; thus all poles of $\mathcal{G}_{0}(z)$ lie on the real line.Löwdin (1982) As $E(\omega)$ is an entire function, this integral converges for any positive value of $\eta$. Numerically, as described below, a value of $\eta$ that minimizes numerical error is used.

II.2 Extensions to open system and wavefunction formalisms

The perturbative analysis described above is readily extended to general open system dynamics and to wavefunction calculations. Below, such extensions are briefly discussed.

Perturbative expansions of the wavefunction, which are significantly more efficient to compute for large closed systems,Rose and Krich (2019) also admit a Laplace solution. The Schrodinger equation describing the evolution of the wavefunction under the action of a time-dependent potential is given by,

[TABLE]

This equation is identical to eq. (3) when the following substitutions are performed,

[TABLE]

It is thus unsurprising that a near identical equation to eq. (24) is obtained upon performing the Laplace transform and inversion,

[TABLE]

where

[TABLE]

and $\ket{\psi_{0}}$ is an eigenstate of $H_{0}$ with eigenenergy $\epsilon_{0}$. The term $E_{\eta}(\omega_{1}-\epsilon_{0}/\hbar)$ is obtained from eq. (15) for the initial state with $G_{0}(i\omega+\eta)\ket{\psi_{0}}=(i\omega+\eta-\epsilon_{0}/i\hbar)^{-1}\ket{\psi_{0}}$. While the wavefunction approach can be computationally more efficient than the density matrix approach for large closed systems, the perturbative expansion suffers from well-known issues when it is used to compute expectation values.Mukamel (1995) For example, a second order expansion of the wavefunction $\ket{\psi(t)}=\ket{\psi_{0}}+\ket{\psi_{1}(t)}+\ket{\psi_{2}(t)}$ yields nine terms when taking the expectation value $\braket{\psi(t)}{O}{\psi(t)}$, with perturbative orders ranging from zero to four and with mixed time ordering. In contrast, the density matrix perturbation treat the state (the density matrix) and the observables on an equal footing.Lavigne and Brumer (2019a)

The FLIPT method can also be used to compute perturbative expansions of open systems. For example, open system dynamics of the Lindblad type have been used by the authors in a study of two-photon control, using the extension described here.Lavigne and Brumer (2019a) Consider for instance the generalized master equation,

[TABLE]

where $\mathcal{K}(t-t^{\prime})$ is the memory kernel superoperator.de Vega and Alonso (2017) A Laplace transform yields a similar expression to eq. (9),

[TABLE]

provided of course that the Laplace transform of the kernel exists,

[TABLE]

Then, a perturbative expansion leads to a set of equations identical to eqs. (12) and (13) above other than having a different Green’s function,

[TABLE]

For the case of a Markovian environment, the memory kernel is in fact memory-less,

[TABLE]

where $\mathcal{R}$ is the relaxation tensor from e.g., the Redfield equation.Weiss (2012) Then, an identical equation to the closed system case above is obtained, with the exception that the eigenvalues of $\mathcal{L}_{0}$ are no longer strictly imaginary but include negative real components.Albert and Jiang (2014) Evaluating $\mathcal{G}_{0}(i\omega+\eta)\rho$ is more difficult but those not affect the overall convergence; the additional broadening of the spectra from decay and decoherence processes makes the integral of eq. (22) better behaved than in the closed system case. The case of non-Markovian dynamics would follow the same approach but requires further consideration as to convergence and is beyond the scope of this paper.

III The FLIPT method

The FLIPT algorithm, introduced below, provides a highly efficient method to numerically evaluate terms of the perturbative series. The numerical inversion of the Laplace transform is performed using a well-known and well-understood Fourier series approach.Dubner and Abate (1968); Veillon (1974); Crump (1976); de Hoog, Knight, and Stokes (1982); Piessens and Huysmans (1984) The Fourier inversion corresponds to using a finite difference grid in the frequency domain, which is simple to implement and has well-known error properties. Importantly, this simple discretization scheme can be used to exploit the iterative structure of the multidimensional frequency domain integrals of eq. (22) and thereby greatly reduce the number of integrand evaluations required.

Below, we describe a Fourier series discretization of the Laplace inversion integral of eq. (24). The discretized integral is expressed as a product of tensors. In the second section below, we show how this tensor form is exploited to obtain the FLIPT algorithm. The numerical complexity and error properties of this algorithm are then described. Finally, we address the numerical evaluation of $\mathcal{V}\rho$ and $\mathcal{G}_{0}(\omega)\rho$ and related performance considerations.

III.1 Tensor product representation

The Laplace transform can be inverted using a Fourier series decomposition in the time domain.Dubner and Abate (1968) Applying a finite difference discretization to eq. (24) gives the following approximation to $\rho_{n}(t)$,

[TABLE]

That is, the Laplace inversion has been approximated by its discrete Fourier series over an interval of size $T$, with a corresponding frequency $\Omega=2\pi/T$.Dubner and Abate (1968); Crump (1976) If $T$ is longer than the length of the field $t_{\text{off}}-t_{\text{on}}$, $E_{\eta}(t)$ can be exactly represented by its Fourier series transform over the interval $T$. Without loss of generality (as described below) the interval is taken here to be $[0,T]$. Then, the Fourier series for $E_{\eta}(t)$ is given by

[TABLE]

with Fourier coefficients given by,Boas (2005)

[TABLE]

The subscript $\eta$ will only be included for those term where $\eta\neq 0$. Since the Fourier series approximation to $E(t)$ is periodic, it has the following Dirac comb as its Fourier transform:

[TABLE]

Substituting eq. (40) into eq. (22) yields a frequency-discretized expression,

[TABLE]

where $k_{i}$ is an integer index for the grid points of the discretized integral over $\omega_{i}$. This equation, the discretized integral that yields the $n$-th order perturbative contribution to $\rho(t)$, is at the core of the FLIPT algorithm. Formulas are given in the Appendix for the computation of spectral quantities.

As previously stated, the start of the interval over which the field is on does not need to be explicitly included. This property is a consequence of the steady initial state of the system. Consider the value of $\rho_{n}(t)$ due to the translated perturbation $E^{\prime}(t)=E(t-t_{0})$. The translated perturbation yields the following Fourier series,

[TABLE]

That is, the translation yields an additional factor of $e^{-i\Omega kt_{0}-\eta t_{0}}$ for each $E_{\eta}[k_{i}]$. These factors generate a corresponding translation of $\rho_{n}(t)$ to $\rho_{n}(t-t_{0})$ by acting on the exponential time-dependence of eq. (37). Hence, $\rho_{n}(t-t_{0})$ due to $E(t-t_{0})$ equals $\rho_{n}(t)$ due to $E(t)$, as is expected from the steady-state initial condition.

Importantly, this invariance to time-translations removes the need to explicitly provide the interval over which the Fourier series of $E(t)$ is computed — this information is entirely encoded in the function $E_{\eta}[k]$. Thus, the only parameters are $T$, the duration of the propagation of $\rho_{n}(t)$ after the field is on, and the convergence parameter $\eta>0$ which, as shown below, can be expressed in terms of $T$.

The specific structure of equation (41) that is responsible for the numerical efficiency of the FLIPT algorithm is exposed here by expressing the discretized integral as a product of tensors. Specifically, frequency indices $k_{1}\cdots k_{n}$ are expressed as an additional index (superscript $k$ below) which is summed over. In this notation, the density matrix at order $n$ from eq. (41) is given by,

[TABLE]

where $\hat{\rho}_{n}^{k}$ is the $n$-th order frequency-resolved (denoted by a caret) density matrix at frequency $\Omega k$,

[TABLE]

The frequency-resolved density matrix $\hat{\rho}_{n}$ is then a 3-tensor with two indices over the quantum mechanical states from the density matrix and one index over the frequency. (The parameter $\eta$ is implicit and is zero for order $n=0$). The zeroth-order term is time-independent and thus is nonzero only at the frequency index $k=0$. Thus, the initial state has the following tensor form,

[TABLE]

The next order of the perturbative expansion is given by first applying $\mathcal{V}$ to $\rho_{0}$, then multiplying the result by $E_{\eta}[k]$ and followed by applying $\mathcal{G}_{0}(i\Omega k+\eta)$. Each of those operations can be expressed by a tensor product with a frequency-resolved superoperator. The tensors $\hat{G}$ and $\hat{V}$ are diagonal in the frequency indices $k,k^{\prime}$ and are given by

[TABLE]

The tensor $\hat{E}_{\eta}$, diagonal in the system state indices, describe the action of $E_{\eta}(\omega)$ on the frequency of the system, with

[TABLE]

In effect, $\hat{E}_{\eta}$ raises or lowers the frequency indices of a state by $k$ for every nonzero $E_{\eta}[k]$. Then, the first order frequency-resolved density matrix is given succinctly by the following tensor product,

[TABLE]

The iterated structure of the perturbative expansion [Eqs. (12) and (13)] translates naturally into a repeated tensor product form. Indeed, higher order terms can be obtained simply by the repeated application of $\hat{G}\hat{E}\hat{V}$. In the tensor notation, the order $n+1$ term is given by

[TABLE]

Note that for $n>1$ $E[k]$ (i.e. with $\eta=0$) is used instead of $E_{\eta}[k]$. Using this notation, the $n$-th order contribution to the density matrix, from equation (41), is a product of the tensors defined above,

[TABLE]

where the square bracket with superscript $k$ denotes the $k$-th frequency element of the overall product of tensors. The tensor form described here not only makes for a convenient notation, it is also much more efficiently evaluated, as shown below.

III.2 Algorithm

The tensor structure of eq. (53) and the uniform discretization of the Fourier series are exploited in FLIPT method to reduce the amount of numerical operations performed to compute eq. (41). This is best described by an example, as done here. The discretized perturbation $E_{\eta}[k]$ from eq. (39) above is taken to be computed at specific values of $\Omega$ and $\eta$ and truncated at frequency indices $L<|k|<U$. The number of grid points of the discretized perturbation is 2$N_{d}$, where $N_{d}=U-L$. In this section, we focus on the frequency integral and thus on the action of the operator $\hat{E}$; the frequency-resolved Green’s function $\hat{G}$ and potential $\hat{V}$ operators are discussed in Sec. III.4.

The frequency-resolved density matrix at order zero consists of only one element, $\rho_{0}$, with frequency index $k=0$. The first order contribution is obtained by performing $\hat{G}\hat{E}_{\eta}\hat{V}\hat{\rho}_{0}$. Specifically, the $k$-th frequency index of $\hat{E}_{\eta}\hat{V}\hat{\rho}_{0}$ is given by the product of $E_{\eta}[k]$ (a scalar) and $\mathcal{V}\rho_{0}$ (the sole nonzero density matrix of $\hat{\rho}_{0}$),

[TABLE]

Then, $\hat{\rho}_{1}=\hat{G}\hat{\rho}_{1^{\prime}}$ yields the first order contribution. Hence, at first order, $2N_{d}$ density matrices are obtained, with frequency indices spanning $[L,U]$ and $[-U,-L]$.

The second order $\hat{\rho}_{2}$ is obtained from the first order contribution $\hat{\rho}_{1}$ by repeating this process. First, $\mathcal{V}_{0}$ is applied to each of the $2N_{d}$ density matrices of $\hat{\rho}_{1}$ to obtain $\hat{\rho}_{1^{\prime\prime}}$. Then the product with $\hat{E}$ is performed as follows,

[TABLE]

That is, every nonzero density matrix $\hat{\rho}_{1^{\prime\prime}}^{k}$ is multiplied with every nonzero values $E[k^{\prime}]$ and summed into $\hat{\rho}_{2^{\prime}}^{k+k^{\prime}}$. This yields the following nonzero frequency indices,

[TABLE]

Thus, $\hat{G}\rho_{2^{\prime}}=\hat{\rho}_{2}$ consists of 8 $N_{d}$ terms.

Importantly, there are $2N_{d}\times 2N_{d}=4N_{d}^{2}$ possible products of $E[k^{\prime}]$ and $\rho_{1^{\prime\prime}}^{k}$, but only $8N_{d}$ distinct values of $k+k^{\prime}$. That is, there are more than one “pathways” to a given frequency index. For example, the density matrix $\hat{\rho}_{1}^{k}$ multiplied by $E[-k^{\prime}]$ contributes to the second order density matrix at index $k-k^{\prime}$; the term $E[k]\hat{\rho}_{1}^{-k^{\prime}}$ also sums to the same frequency index. Contributions such as these to the same final frequency index are summed over as soon as they are available. Significantly, this summing is responsible for the exponential speedup of the FLIPT algorithm over conventional numerical quadrature.

Specifically, the number of nonzero frequency indices of $\hat{\rho}_{n}$ is linear in $N_{d}$ as a consequence of the large number of repeated indices. The frequency indices where $\hat{\rho}_{n}$ is nonzero are given by,

[TABLE]

At order $n$, there are $n+1$ such intervals, each with $nN_{d}$ nonzero frequency indices.333 The regions are taken here to be disjoint. When the regions overlap, fewer points are needed and the numerical effort is reduced. The number of distinct frequency indices is thus $(n+1)nN_{d}$, i.e. a linear function of $N_{d}$.

In contrast, if the frequency discretization is performed using an arbitrary set of grid points $\{\omega_{i}\}$, the number of nonzero elements of $\hat{\rho}_{n}$ would scale exponentially with $n$. Indeed, given an arbitrary set of frequency points $\{\omega_{k}\}$ of size $2N_{d}$, the set of points $\omega_{i}+\omega_{j}$ for all $i$ and $j$ (i.e. the number of frequency indices of $\hat{\rho}_{2}$ from a general, non-uniform grid) contains $4N_{d}^{2}$ unique frequencies. The third order result (i.e. the number of frequency indices of $\hat{\rho}_{3}$) contains $2N_{d}\times 4N_{d}^{2}=8N_{d}^{3}$ unique frequencies, the fourth order result contains $32N_{d}^{4}$ etc. That is, the number of frequency indices would increase exponentially. However, the Fourier series representation (and resulting uniform grid) used here bypasses the exponential increase without introducing artificial coarse-graining or filtering of any sort due to the summing of repeated indices.

III.3 Numerical analysis

Formally, the FLIPT algorithm scales as $O(N_{d}^{2})$ at all orders $n$ of the perturbation, instead of the exponential $O(N_{d}^{n})$ scaling of standard quadrature. Consider the standard quadrature of eq. (22) using a grid with the $2N_{d}$ of points. The integrand for the $n$-th perturbative contribution has $n$ applications of both $\mathcal{V}$ and $\mathcal{G}_{0}(i\omega+\eta)$. Denoting the number of operations required to perform $\mathcal{V}\rho$ and $\mathcal{G}_{0}(i\omega+\eta)\rho$ by $N_{V}$ and $N_{G}$ respectively, the number of operations required to compute the integrand once is $n(N_{v}+N_{G})$. The integrand has to be evaluated at every multidimensional quadrature points. At order $n$, computing the $n$-dimensional quadrature using a one-dimensional grid of $2N_{d}$ points require $(2N_{d})^{n}$ integrand evaluation. Thus, the number of operations required by standard quadrature obeys,

[TABLE]

Therefore, standard quadrature integration scales as $N_{d}^{n}$ — exponential in the perturbative order $n$.

In contrast, the FLIPT algorithm presented above requires a number of operation at most quadratic in the number of discretization points $N_{d}$ at all orders. The equivalent number of operations required to compute the $n$-th perturbative order can be derived from eq. (53) above,

[TABLE]

The first term corresponds to the number of operations required to obtain the $n-1$ order from which order $n$ is computed. The remaining terms are the number of operations required to compute tensor products with $\hat{V}$, $\hat{E}$ and $\hat{G}$ respectively. In most applications, the highest order of the perturbative expansion dominates the CPU time and other terms can be ignored, such that the following complexity is obtained,

[TABLE]

For a system of $d$ levels, $N_{E}\propto d^{2}$, as will be shown below. Depending on the number of operations required by $\hat{V}$ and $\hat{G}$, the FLIPT algorithm is linear to quadratic in the number of frequency points $N_{d}$, an enormous improvement over standard quadrature for perturbative order 3 and above.

The FLIPT method is well-behaved numerically as the Fourier series representation is numerically “exact” for continuous and time-limited $E(t)$.Crump (1976) Furthermore, the summation procedure over repeated indices described above is analytic and not the result of any numerical coarse-graining. The numerical error of a function $y(t)$ computed by Fourier-Laplace inversion, such as e.g. a matrix element of $\rho(t)$ computed using the FLIPT algorithm, obeys the following relation,Dubner and Abate (1968)

[TABLE]

where $C=\max y(t)$ over an interval from 0 to $T$. Hence, the error increases on the approach to $T$ in proportion to $\eta$. Error analysis can be used to compute an optimal value for $\eta$. Following Ref. Piessens and Huysmans, 1984, $\eta$ is taken here to be $2\pi/T=\Omega$, which provides a good balance between the Laplace inversion error and any floating point errors arising from evaluating the Green’s function near poles of $\mathcal{L}_{0}$. For $t<T/2$ the relative error is $\leq 10^{-3}$. It should be noted that eq. (60) is an upper bound; in practice, the recurrence of the Fourier series at $T$ dwarfs any numerical errors due to the inversion.

III.4 Tensor operations

The numerical performance of the FLIPT method depends crucially on a fast computation of the tensor products of eq. (41). Below, the numerical implementation and performance properties of those operations are discussed for a system of $d$ levels. It should be noted that the same operations also determine the performance of other numerical methods; non-perturbative propagation methods, for example, also rely on the fast evaluation of $\mathcal{V}\rho$ and $\mathcal{L}_{0}$.

Computing $\hat{\rho}^{\prime}=\hat{V}\hat{\rho}$ is done by a straightforward computation of $\mathcal{V}\rho$ at each frequency index $k$ (the number of which is denoted $N_{k}$) where $\hat{\rho}^{k}$ is nonzero,

[TABLE]

Standard matrix-matrix multiplication is used for this step. The computational complexity of $\hat{V}\hat{\rho}$ is given by the complexity of performing the $\mu\rho$ matrix-matrix multiplication ($N_{v}\propto d^{3}$ above) multiplied by $N_{k}$.

Similarly, $\hat{\rho}^{\prime}=\hat{G}\hat{\rho}$ is evaluated by computing $\hat{\rho}^{\prime k}=\mathcal{G}_{0}(i\Omega k+\eta)\hat{\rho}^{k}$ for all values of the frequency index $k$. For a closed $d$-level system with energy eigenvalues $E_{i}$, $\mathcal{G}_{0}(i\omega+\eta)\rho$ is given by the following analytical formula,Löwdin (1982)

[TABLE]

where $\rho_{ij}$ is the $i,j$ matrix element of $\rho$ in the eigenbasis of $H_{0}$ of eq. (4). Thus, for such a system, the number of operations per frequency index $N_{G}=O(d^{2})$ and the overall complexity of computing $\hat{\rho}^{\prime}$ is $O(N_{k}d^{2})$.

The Green’s function can still be evaluated even when $\mathcal{L}_{0}$ does not have the form of eq. (4) or when the Hamiltonian is not diagonal. This is the case in, e.g., Ref. Lavigne and Brumer, 2019a where the Redfield equation is used. The Green’s function applied to a density matrix $\rho$ yields

[TABLE]

Multiplying both sides by $(z-\mathcal{L}_{0})$, the following is obtained,

[TABLE]

That is, the product of the Green’s function on $\rho$ can be obtained by solving the above system of linear equations for $\rho^{\prime}$. Using an iterative method (such as GMRES,Baker, Jessup, and Manteuffel (2005) used in Ref.Lavigne and Brumer, 2019a) and an efficient algorithm for the product $\mathcal{L}_{0}\rho$ (such as the Pollard-Friesner methodPollard and Friesner (1994)), a numerical complexity of $O(d^{3})$ is obtained. For small systems, direct methods (such as the LU decompositionAnderson et al. (1999)) can also be used.

Finally, $\hat{\rho}^{\prime}=\hat{E}\hat{\rho}$ is computed using a matrix-matrix product over frequency indices,

[TABLE]

where $\hat{E}^{k^{\prime},k}$ is given by eq. (50). This operation has a numerical complexity of $O(N_{k}^{\prime}N_{k}d^{2})$ where $N_{k}$ and $N_{k}^{\prime}$ are the number of indices $k$ and $k^{\prime}$ where $\hat{\rho}^{k}$ and $\hat{\rho}^{\prime k^{\prime}}$ respectively are nonzero. This is done numerically by a matrix-matrix multiplication with $\hat{E}$ expressed as a $N_{k}^{\prime}\times N_{k}$ matrix and $\hat{\rho}$ expressed as a $N_{k}\times d^{2}$ matrix.

Importantly, it should be noted that the quadratic scaling described in this paper is for the case where the computation of $\hat{E}$ is the most expensive step. This is true when the number of levels is low and $\mathcal{V}\rho$ and $\mathcal{G}_{0}(\omega)\rho$ are relatively inexpensive. For larger systems, this is no longer the case and the FLIPT algorithm becomes dominated by the linear in $N_{d}$ terms of eq. (59); in those cases, the performance advantage of FLIPT over other quadrature method is significant even for second order processes such as linear absorption.

IV Implementation and example calculations

The FLIPT algorithm is easily implemented using tensor algebra. Such an implementation for closed and open secular dynamics is made freely available by the authors.Lavigne (2019) The resultant code is very short (less than 600 lines in the Julia programming language) and can thus be easily translated to other programming languages and environments.

The implementation can be used as a “black box” to compute arbitrary orders of the perturbative expansion as well as arbitrarily complex nonlinear spectroscopy diagrams. The implementation includes subroutines to build $\hat{V}$ and $\hat{G}$ from $\mu$ and $H_{0}$ for closed systems. The ability to optionally operate on only the ket or bra sides of the density matrix with either solely positive or negative frequency components of the field can be used to compute the response of arbitrary nonlinear spectroscopy diagrams including $\bm{k}$-vector phase matching.Mukamel (1995) The maximum propagation time $T$, the energy levels, the dipole transition operator and the electric field $E(\omega)$ are the only required parameters; the implementation is otherwise fully automatic.

We have recently applied the FLIPT algorithm to the simulation of quantum dynamics arising from two-photon absorption in retinal to understand the mechanism of a quantum control experiment on living brain cells.Paul et al. (2017); Lavigne and Brumer (2019a) Here, we focus on the performance properties of the FLIPT algorithm, and the application of the FLIPT method to spectroscopy. First, performance is studied using a small model system; results are compared with non-perturbative propagation. Then, a sample calculation on a model molecule, pyrazine, shows how the FLIPT method can be applied to the simulation of spectroscopy experiments. The model exhibits an ultrafast conical intersection mediated population transfer between two electronic states. A signature of the conical intersection is present in a simulated transient absorption spectrum obtained from a third-order perturbative calculation.

IV.1 Numerical scaling and performance

The numerical performance and numerical error of the FLIPT method are studied here using a four-level model system (Fig. 1(a)). This small model exhibits slow ($\approx 100$ fs) and fast dynamics ($\approx 1fs$) dynamics resulting from separated manifolds $\ket{g}$ and $\ket{e}$, each consisting of two closely spaced levels. The system interacts with an ultrafast, coherent pulse. The spectrum of this pulse is given by,

[TABLE]

where the function $g(\omega;\omega_{0},\sigma)$ is a dimensionless, normalized Gaussian function centered at frequency $\omega_{0}$,

[TABLE]

The standard deviation of the field $\sigma$ in the Fourier domain is obtained from the FWHM in the time domain as,

[TABLE]

The FLIPT results are compared to the result of a non-perturbative propagation of the density matrix. A perturbative decomposition of the non-perturbative result is approximated with a least-squares fit at multiple values of the field intensity. Computations are performed using one core of an Intel Xeon E5-2630 (2.20GHz) processor, with parameters given in the caption of Fig. 1.

It should be noted that no artificial decoherence or broadening parameters are introduced — the only parameter of the FLIPT method is the propagation time $T$. The value of $T$ sets the maximum duration of a valid FLIPT result; indeed, the Fourier series decomposition imply a recurrence of the dynamics at time $T$ after the excitation. Away from the recurrence at $T$, the FLIPT method is highly accurate even when compared to the non-perturbative method, as shown in Fig. 1(b). In the quadratic case, the relative error between the FLIPT result and non-perturbative propagation is less than $10^{-3}$ and stable until at least $T/2$. The quartic case is similar; the higher relative error is due to the order extrapolation of the non-perturbative propagation.

The linear to quadratic scaling with respect to the number of frequency points given in eq. (59) is demonstrated in Fig. 2, where the duration $T$ (and thus the number of points $N_{d}$) of the propagation is adjusted. The time required to compute each of the first four perturbative orders of $\rho(t)$ are shown. As described above, the complexity is at most quadratic for all perturbative orders, making the FLIPT method highly suited to the computation of higher-order processes.

The FLIPT method is uniquely suited to study processes where fast and slow timescales (e.g., electronic and nuclear motion) are both present. The spectral nature of the method leads to an important property: time-scale invariance. This sets it apart from non-perturbative propagation and other time-dependent methods.Rose and Krich (2019); Arango and Brumer (2013); Zhuang et al. (2006) Indeed, the FLIPT method is independent of the energy spacing of the system under investigation, since that spacing is fully subsumed in the smooth function $\mathcal{G}_{0}(\omega+i\eta)$. The only timescales of importance are those of the interaction and of the duration of propagation $T$. The number of points $N_{d}$ required is directly proportional to the bandwidth of $E(\omega)$ and inversely proportional to the duration of the simulation, a result similar to that obtained with the rotating wave approximation, but here without any approximations. Hence, computing dynamics over 300 ps arising from a 10 ps interaction with frequency $\omega_{0}$ requires the same computation time as computing 300 fs of dynamics arising from a 10 fs pulse. In contrast, a non-perturbative propagation of these cases will scale at least linearly as the minimal time step is set by the fast, optical transients of the electric field. This property is demonstrated in Fig. 3.

IV.2 Sample calculation: spectroscopy of pyrazine

Pyrazine, a small molecule, is a well-known model system to study ultrafast internal conversion.Raab et al. (1999); Sukharev and Seideman (2004); Christopher, Shapiro, and Brumer (2005, 2006); Ryabinkin, Joubert-Doriol, and Izmaylov (2014) A model for pyrazine with three electronic surfaces (denoted $S_{0}$, $S_{1}$ and $S_{2}$) and two vibrational modes (denoted $6a$ and $10a$) is used here to demonstrate the use of FLIPT with a sizable multilevel system (with $\approx 300$ levels). The Hamiltonian for this system is given by,

[TABLE]

The operators $p_{i}$ and $q_{i}$ are the momentum and position operators for the modes in mass-weighted coordinates, $\epsilon_{i}$ is the Franck-Condon energy of surface $i$, $\kappa_{i}$ is the displacement of the tuning mode $6a$ on surface $i$ and $\lambda$ is the non-adiabatic coupling between surfaces $S_{1}$ and $S_{2}$. The transition dipole operator $\mu$ couples the ground and excited surfaces,

[TABLE]

where $E_{0}\mu_{0}=1eV$ is a scaling parameter for the perturbative expansion. Parameters for this model are given in Table 1. Direct diagonalization yields $\approx 300$ basis states. Gaussian pulses, described by eqs. (66)-(68), are used throughout.

A measurement of the linear absorption is simulated using the FLIPT algorithm (Fig. 4). The absorption is computed from the heterodyne detection formula given by eq. (84) in the Appendix,

[TABLE]

where $\omega_{0}$ is the central frequency of the field $E(\omega)$, a Gaussian pulse with a FWHM of 300 fs. The central frequency of the pulse is swept to obtain the absorbance,

[TABLE]

The resolution of the absorbance spectrum is directly proportional to the length of the (minimum uncertainty) probe pulses. This is shown in the inset of Fig. 4, where the absorbance spectrum obtained with pulses with a FWHM of 1000 fs is compared with that obtained with the above shorter pulses. Using pulses that are longer in the time-domain and narrower in the frequency domain significantly increases the resolution of the spectral peaks, a finite pulse effect.Smallwood, Autry, and Cundiff (2017); Perlík, Hauer, and Šanda (2017) It should be noted that both results are computed using the same number of discrete frequency points per pulse and thus require the same amount of computing time. Furthermore, no additional decoherence processes, phenomenological broadening or system-bath interactions are added: the “smoothness” of the spectrum is entirely due to the limited resolution of the probe laser. This should be contrasted with other common methods of computing the spectrum, e.g., through the Fourier transform of the autocorrelation functionRaab et al. (1999) or the response of the system to a CW field, where ad-hoc broadening factors or signal windowing are required to obtain numerical convergence.

A pump-probe transient absorption (TA) spectrum is computed as an example of the kind of higher order process that can be studied with the FLIPT method. In a transient absorption experiment, a linear absorption spectra is measure with an ultrashort probe pulse after excitation with an ultrashort pump pulse and plotted as a function of the delay $\tau$ between the pump and probe pulses. This is a modeled four-wave mixing experiment and thus third order in the perturbative series; a three-dimensional integral is computed at every value of the pump-probe delay.

Here neither the pump nor the probe are approximated; a realistic pulse shape is used in both case. The pump pulse has a FWHM of 20 fs and is centered at 4.8 eV, the absorption maximum of $S_{2}$. The probe pulse $E_{\text{probe}}(\omega)$ has a FWHM of 5 fs and a central frequency of 4.3 eV; its bandwidth is sufficient to probe both the $S_{1}$ manifold between 3.7 and 4.2 eV and the $S_{1}$ manifold between 4.5 and 5.5 eV. This simulation require in this case a discretization of 180 points (positive and negative) for the pump and 360 points (positive only) for the probe. Using the same discretization with a standard quadrature calculation would require more than 10 million integrand evaluations, each computed using seven matrix-matrix multiplications of $\mu$ and $\rho$. In comparison, the FLIPT result is evaluated using only 1000 such multiplications — a factor of $10^{4}$ improvement. This integral is repeated for 300 different values of the pump-probe delay.

The absorbance measured with the probe pulse is given by,

[TABLE]

where the heterodyne intensity is computed from eq. (88) and the probe intensity is given by $I_{\text{probe}}(\omega_{p})=|E_{\text{probe}}(\omega_{p})|^{2}$. The absorbance measured with the probe pulse after excitation with the pump pulse is given by,

[TABLE]

where $I_{\text{het,pp}}(\omega_{p};\tau)$ is the heterodyne-detected absorption of the probe pulse following excitation with the pump pulse,444Due to the strict time-ordering of the perturbative expansion, $I_{\text{het,pp}}(\omega;\tau)$ is zero when the probe arrives before the pump, i.e. it is background-free. a fourth order perturbative term. The transient absorption spectrum is the difference of these two quantities,

[TABLE]

A simulated TA spectrum is shown in Fig. 5. Red and blue colors represented heightened and lowered absorption compared to the unpumped system. The pump pulse excites a coherent wavepacket on the $S_{1}$ and $S_{2}$ surfaces; wavepacket dynamics on those surfaces generate the time-dependent transient absorption. Fig. 5 shows the initial excitation at delay zero, followed by wavepacket oscillations on the $S_{2}$ surface at $\hbar\omega_{p}\approx$ 4.7 eV and on the $S_{1}$ surface at $\hbar\omega_{p}\approx$ 3.8 eV due to the tuning mode $6a$. The structured bands in the spectrum are due to vibronic states of the two electronic surfaces. As the model includes no decoherence, there is no decay of the signal.

The TA spectrum reveals details of the two excited electronic surface, including the presence of a conical intersection. The first 80 fs of the TA spectrum in the energy region of the $S_{2}$ surface are highlighted in Fig. 6(a). Qualitatively, a band of increasing energy (along the solid curve marked B) can be identified that splits off at $\tau\approx 40$ fs from the main energy oscillations around 4.8 eV (along dashed curves marked A and C).

These two features are produced by the wavepacket motion along the tuning mode $6a$, shown diagrammatically in Fig. 6(b). Initial excitation by the pump generates a coherent wavepacket at the Franck-Condon (FC) point. Then, the wavepacket moves on the $S_{2}$ surface towards the conical intersection (CI) between the $S_{1}$ and $S_{2}$ surface along the path marked by A in both figures. Upon encountering the CI, the wavepacket bifurcates. Part of the wavepacket moves past the CI and continues along the $S_{2}$ surface, producing the energy increasing feature along B, while some of the wavepacket transits to $S_{1}$ and reflects back onto $S_{2}$, which yields the energy oscillation shown by C. Finally, some of the population transits onto $S_{1}$ and moves to the $S_{1}$ minimum. This latter contribution is mostly dark at those energies, but is responsible for the transient absorption signal below the $S_{2}$ minimum.

The paths A, B and C in Fig. 6(a) are obtained from the $S_{2}$ potential energy $V_{2}(q(\tau))$ as a function of the normal mode $q(\tau)=\omega\tau$, where $q(0)=0$ is the FC point. That is, the paths are approximations of the potential energy of a ballistic particle moving on $S_{2}$. Specifically, this qualitative treatment shows how a conical intersection could be experimentally identified from the bifurcation of a wavepacket as measured in a transient absorption experiment. Importantly, all parameters leading to this result have clear physical origins; the fully automatic FLIPT computation requires no broadening or convergence factors not directly related to physical properties of the molecule and radiation.

V Conclusion

Time-dependent perturbation theory plays a central role in applications of quantum mechanics. However, higher order perturbative contributions are difficult to evaluate numerically. Inversion of the Laplace transform by Fourier seriesDubner and Abate (1968); Crump (1976); Piessens and Huysmans (1984) is used to compute arbitrary order time-dependent perturbative expansions of the Liouville-von Neumann equation with analytically demonstrated convergence for perturbations of finite duration. Importantly, nonlinear spectroscopyCohen-Tannoudji, Dupont-Roc, and Grynberg (1992); Mukamel (1995) and quantum control experimentsShapiro and Brumer (2012) performed with pulsed lasers are well described by perturbation theory with finite duration fields. Here we have introduced the FLIPT algorithm that uses the particular structure of the perturbative expansion to efficiently compute the Fourier-Laplace inversion. It is a numerically exact scheme to compute the Fourier-Laplace inversion of the perturbative series, that has significant advantages over other integration methods and over non-perturbative propagation methods. Compared with a propagation of the density matrix, the FLIPT method yields perturbative results that are more readily understood and experimentally applicable. Since it is a spectral method, its computational complexity is independent of the fast and slow inherent timescales of the system. In contrast, time-dependent propagation methods require fine timesteps when the RWA breaks down,Arango and Brumer (2013); Rose and Krich (2019) as is the case in the presence of significant nonresonant processes.

The iterated structure of the perturbative expansion can be expressed succinctly using tensor algebra; exploiting this structure, the FLIPT method achieves an exponential speedup with respect to the perturbation order over standard multidimensional quadrature.Novak and Ritter (1997) Indeed, the computational complexity of the algorithm is at most quadratic at any order $n$ of the perturbative expansion, i.e. the $N$-point discretized multidimensional integration can be computed in $O(N^{\alpha})$ operation with $\alpha\lesssim 2$. That is significantly better than the $O(N_{d}^{n})$ scaling of standard quadrature methods. Furthermore, the obtained fixed-grid spectral representation is easily Fourier transformed to the time-domain, a property not shared by quadrature and Monte Carlo integration.

Other methods for simulating the interaction of light with matter include the SPECTRON programZhuang et al. (2006) and the NISE method,Torii (2006) as well as symbolic or analytical approaches.Perlík, Hauer, and Šanda (2017); Smallwood, Autry, and Cundiff (2017) SPECTRON and NISE are significantly less general than the FLIPT method, since they are limited to computing responses of specific types of linear to quartic order spectroscopies. Analytical approaches are limited to specific idealized pulse shapes and must be re-derived, with significant effort, for every perturbation order. In contrast, the FLIPT algorithm can be used to compute any observable quantities (i.e. not only responses but also electronic populations, vibrational displacements, etc.), at arbitrary order of the perturbation theory and for arbitrary pulse shapes.

The recently proposed UFF method of Ref. Rose and Krich, 2019 is in many respects similar to the FLIPT algorithm. UFF is an arbitrary order, arbitrary interaction time-domain approach developed for wavefunctions. Hence, UFF scales significantly better than the FLIPT method with respect to the size of the perturbed system. However, it has convergence issues with respect to fast oscillations in the absence of the RWA and is limited to energy-diagonalized Hamiltonian systems. The choice between the UFF and FLIPT algorithm should be made based on considerations such as the size of the system under study, the importance of nonresonant processes and the order of the perturbative expansion.

Future work should focus on extending the applicability of the FLIPT algorithm. Here, the high performance of the method is achieved at the cost of a high memory usage. Indeed, for large systems, a copy of the density matrix must be stored at each discretized frequency index, with a correspondingly large memory cost. The use of non-Markovian equations can significantly reduce memory usage, as can the propagation of wavefunctions instead of density matrices. Such extensions are under development. Finally, although the primary focus of this paper has been on light-matter interaction, the FLIPT method can conceivably be applied to other kinds of perturbation theory. In particular, it is closely related to some path integral methods with similar recursive integral structures.Thirumalai, Bruskin, and Berne (1983); Shao and Makri (2001); Jadhao and Makri (2008) Applications to other systems are being investigated.

Acknowledgments: This work was supported by the U.S. Air Force Office of Scientific Research under Contract No. FA9550-17-1-0310, and by the Natural Sciences and Engineering Research Council of Canada.

Appendix A Spectral quantities from FLIPT result

The Fourier-Laplace inversion can be used to compute spectrally resolved quantities without introducing ad-hoc broadening or decoherence factors, but care must be taken to ensure convergence. Indeed, consider the expectation value of an operator $O$ evaluated using the frequency-resolved, $n$-th order density matrix at index $k$ [eq. (46) above],

[TABLE]

This quantity is not the Fourier transform of $O_{n}(t)=\text{Tr}O\rho_{n}(t)$ at the frequency $\Omega k$. Specifically, $\hat{O}_{n}^{k}$ is given by a convergent Fourier series over a finite interval while the Fourier transform of $O_{n}(t)$ may not even exist. For example, in the absence of an environment, oscillatory coherences do not decay and the Fourier transform of an oscillatory expectation value does not converge.Lavigne and Brumer (2017) However, approximate Fourier transforms can be computed using the FLIPT method as shown below.

First, the spectrum of an observable $O$ at perturbative order $n$ can be approximated by taking the Fourier transform of the time-dependent value $O_{n}(t)$ over a finite interval. The obtained spectrally resolved expectation value $O_{T,n}(\omega)$, defined below, is an approximation to the true Fourier transform.Wiener (1930) Using the Fast Fourier Transform to compute $O_{T,n}(\omega)$ from $O_{n}(t)$ is expensive in the presence of high-frequency components since it requires $O_{n}(t)$ to be meshed over a fine grid. Fortunately, $O_{T,n}(\omega)$ can be computed directly from the FLIPT result at the grid points $\omega=\Omega k$ using the Fourier series,

[TABLE]

The expectation value $O_{T,n}(\Omega k)$ converges as a distribution to the Fourier transform when $T\rightarrow\infty$.Rudin (1987) For the case where $\eta=\Omega$, as is done in the present implementation, the above further simplifies to,

[TABLE]

Second, specific formulas can be obtained for spectroscopic signals detected through heterodyning. This is the case in many ultrafast spectroscopy experiments, e.g., transient absorption, pump-probe and 2D spectroscopy. A heterodyne signal $I^{(n)}_{\text{het}}(t)$ is obtained by mixing the response of the system, given by an observable $\mu_{n}(t)$, with a probe electric field $E(t)$ and detected in the direction of the probe,Mukamel (1995)

[TABLE]

where the subscripts $+$ denote that only positive frequency components (i.e. positive phase-matched $\bm{k}$ components) are detected. The heterodyne signal can be computed from the Fourier series representation as follows,

[TABLE]

The integrated heterodyne signal can then be approximated as above by a finite time integral,

[TABLE]

This is the signal as measured in Fig. 4 or in a pump-probe experiment. The heterodyne signal can also be dispersed through a monochromator to obtain a frequency-resolved measurement, as is done in transient absorption spectroscopy. The monochromated quantities at output frequency $\Omega k_{\text{out}}$ are given by,

[TABLE]

where $k_{\text{out}}$ is the frequency of the monochromator. Then, the integrated signal is given by

[TABLE]

This formula is used to obtain the transient absorption spectra above.