Starting from a set of orbitals localized in \(A\) at \(t=0\) it is possible to derive a time-propagation scheme with time step \(\Delta t\) by recursively applying the discrete time-evolution operator \(\hat{U}(\Delta t)\equiv \hat{U}(t+\Delta t,t)\) and splitting the components with eq. \eqref{eq:mask_split}. The result can be written in a closed form for \(\phi^A_i(\vec{r},t)\), represented in real space, and \(\phi^B_i(\vec{k},t)\), in momentum space, with the following structure: \[\begin{aligned} \label{eq:FMM_prop} \begin{array}{l} \phi^A_i(\vec{r},t+\Delta t) = \varphi^A_i(\vec{r},t+\Delta t) +\varphi^B_i(\vec{r},t+\Delta t)\ ,\\ \phi^B_i(\vec{k},t+\Delta t) = \vartheta^A_i(\vec{k},t+\Delta t)+\vartheta^B_i(\vec{k},t+\Delta t)\ , \end{array}\end{aligned}\] and the additional set of equations, \[\begin{aligned} \label{eq:FMM_prop_aux} \varphi^A_i(\vec{r},t+\Delta t)& = M \hat{U}(\Delta t) \phi^A_i(\vec{r},t)\ ,\\ \varphi^B_i(\vec{r},t+\Delta t)& = \frac{M}{(2\pi)^{3/2}}\int {\rm d}\vec{k} \mathrm{e}^{\mathrm{i}\vec{k}\cdot\vec{r}} \hat{U}_{\rm v}(\Delta t) \phi^B_i(\vec{k},t) \ ,\\ \vartheta^A_i(\vec{k},t+\Delta t)& = \frac1{(2\pi)^{3/2}} \int {\rm d}\vec{r} \mathrm{e}^{-\mathrm{i}\vec{k}\cdot\vec{r}} (1-M) \hat{U}(\Delta t) \phi^A_i(\vec{r},t) \ ,\\ \vartheta^B_i(\vec{k},t+\Delta t)& = \hat{U}_{\rm v}(\Delta t) \phi^B_i(\vec{k},t) \nonumber \\ &- \frac1{(2\pi)^{3/2}} \int {\rm d}\vec{r} \mathrm{e}^{-\mathrm{i}\vec{k}\cdot\vec{r}} \varphi^B_i(\vec{r},t+\Delta t)\ .\end{aligned}\] The momentum-resolved photoelectron probability is then obtained directly from the momentum components as \cite{DeGiovannini_2012} \[P(\vec{k})=\lim_{t\rightarrow \infty}\sum_i^N |\phi^B_i(\vec{k},t)|^2\,,\] while the energy-resolved probability follows by direct integration, \(P(E)=\int_{E=|\vec{k}|^2/2}{\rm d}\vec{k}P(\vec{k})\).

In eq. \eqref{eq:FMM_prop_aux} we introduced the Volkov propagator \(\hat{U}_{\rm v}(\Delta t)\) for the wavefunctions in \(B\). It is the time-evolution operator associated with the Hamiltonian \(\hat{H}_{\rm v}\) describing free electrons in an oscillating field. Given a time dependent vector field \(\vec{{A}}(t)\), the Hamiltonian \(\hat{H}_{\rm v}=\frac{1}{2}\left(-i\vec{\nabla}-\frac{\vec{{A}}(t)}{c}\right)^2\) expressed in the velocity gauge is diagonal in momentum and can be naturally applied to \(\phi^B_i(\vec{k},t)\).

For all systems that can be described by a Hamiltonian such that \(\hat{H}(\vec{r},t)=\hat{H}_{\rm v}(\vec{r},t)\) for \(\vec{r} \in B\) and all time \(t\), eqs. \eqref{eq:FMM_prop} and \eqref{eq:FMM_prop_aux} are equivalent to a time propagation in the entire space \(A\cup B\). In particular, it exactly describes situations where the electrons follow trajectories crossing the boundary separating \(A\) and \(B\) as illustrated in Fig. \ref{fig:pes_sheme}(b).

In Octopus we discretize eq. \eqref{eq:FMM_prop_aux} in real and momentum space and co-propagate the complete set of orbitals \(\phi^A_i(\vec{r},t)\) and \(\phi^B_i(\vec{k},t)\). The propagation has to take care of additional details since the discretization can introduce numerical instability. In fact, substituting the Fourier integrals in \eqref{eq:FMM_prop_aux} with Fourier sums (usually evaluated with FFTs) imposes periodic boundary conditions that spuriously reintroduces charge that was supposed to disappear. This is illustrated with a one-dimensional example in Fig. \ref{fig:pes_nfft}(a) where a wavepacket launched towards the left edge of the simulation box reappears from the other edge.

An alternative discretization strategy is zero padding. This is done by embedding the system into a simulation box enlarged by a factor \(\alpha>1\), extending the orbitals with zeros in the outer region as shown in Fig. \ref{fig:pes_nfft}(b). In this way, the periodic boundaries are pushed away from the simulation box and the wavepackets have to travel an additional distance \(2(\alpha -1)L\) before reappearing from the other side. In doing so, the computational cost is increased by adding \((\alpha -1)n\) points for each orbital.

This cost can be greatly reduced using a special grid with only two additional points placed at \(\pm \alpha L\) as shown in Fig. \ref{fig:pes_nfft}(c). Since the new grid has non uniform spacing a non-equispaced FFT (NFFT) is used \cite{Kunis_2006,Keiner_2009}. With this strategy, a price is paid in momentum space where the maximum momentum \(k_{\rm max}\) is reduced by a factor \(\alpha\) compared to ordinary FFT. In Octopus we implemented all three strategies: bare FFT, zero padding with FFT and zero padding with NFFT.

All these discretization strategies are numerically stable for a propagation time approximately equivalent to the time that it takes for a wavepacket with the highest momentum considered to be reintroduced in the simulation box. For longer times we can employ a modified set of equations. It can be derived from \eqref{eq:MM_prop_aux} under the assumption that the electron flow is only outgoing. In this case we can drop the equation for \(\varphi^B_i\) responsible for the ingoing flow and obtain the set \[\begin{aligned} \label{eq:MM_prop_aux} \begin{array}{l} \varphi^A_i(\vec{r},t+\Delta t) = M \hat{U}(\Delta t) \phi^A_i(\vec{r},t)\ ,\\ \varphi^B_i(\vec{r},t+\Delta t) = 0 \ ,\\ \vartheta^A_i(\vec{k},t+\Delta t) = \frac1{ (2\pi)^{{3}/{2}}} \int {\rm d}\vec{r} e^{-\mathrm{i}\vec{k}\cdot\vec{r}} (1-M) \hat{U}(\Delta t) \phi^A_i(\vec{r},t) \ ,\\ \vartheta^B_i(\vec{k},t+\Delta t) = \hat{U}_{\rm v}(\Delta t) \phi^B_i(\vec{k},t)\ . \end{array}\end{aligned}\] This new set of equations together with \eqref{eq:FMM_prop} lifts the periodic conditions at the boundaries and secures numerical stability for arbitrary long time propagations. A consequence of this approximation is the fact that the removal of charge is performed only in the equation for \(\varphi^A_i\) by means of a multiplication by \(M(\vec{r})\). This is equivalent to the use of a mask function boundary absorber that is known to prevent reflections in an energy range that depends on \(M(\vec{r})\) \cite{DeGiovannini:2014wo}. Carefully choosing the most appropriate mask function thus becomes of key importance in order to obtain accurate results.

We conclude briefly summarizing some of the most important features and applications of our approach. The method allows us to retrieve \(P(\vec{k})\), the most resolved quantity available in experiments nowadays. In addition, it is very flexible with respect to the definition of the external field and can operate in a wide range of situations. In the strong-field regime, it can handle interesting situations, for instance, when the electrons follow trajectories extending beyond the simulation box, or when the target system is a large molecule. This constitutes a step forward compared to the standard theoretical tools employed in the field which, in the large majority of cases, invoke the single-active-electron approximation. In Ref. \cite{DeGiovannini_2012} the code was successfully employed to study the photoelectron angular distributions of nitrogen dimers under a strong infrared laser field. The method can efficiently describe situations where more than one laser pulse is involved. This includes, for instance, time-resolved measurements where pump and probe setups are employed. In Ref. \cite{DeGiovannini_2013} Octopus was used to monitor the time evolution of the \(\pi\rightarrow\pi^*\) transition in ethylene molecules with photoelectrons. The study was later extended to include the effect of moving ions at the classical level \cite{CrawfordUranga_2014}. Finally, we point out that our method is by no means restricted to the study of light-induced ionization but can be applied to characterize ionization induced by other processes, for example, ionization taking place after a proton collision.