# Phase delays in $\omega-2\omega$ above-threshold ionization

S. D. López1, S. Donsa2, S. Nagele2, D. G. Arbó1,3 and J. Burgdörfer2 1
Institute for Astronomy and Space Physics - IAFE (CONICET-UBA), CC 67, Suc.
28, C1428ZAA, Buenos Aires, Argentina  2Institute for Theoretical Physics,
Vienna University of Technology, Wiedner Hauptstr. 8-10/E136, A-1040 Vienna,
Austria, EU  3Facultad de Ciencias Exactas y Naturales y Ciclo Básico Común,
Universidad de Buenos Aires, Argentina

###### Abstract

Relative phases of atomic above-threshold ionization wavepackets have been
investigated in a recent experiment [L. J. Zipp, A. Natan, and P. H.
Bucksbaum, Optica 1, 361-364 (2014)] exploiting interferences between
different pathways in a weak probe field at half the frequency of the strong
ionization pulse. In this work we theoretically explore the extraction of
phase delays and time delays of attosecond wavepackets formed in strong-field
ionization. We perform simulations solving the time-dependent Schrödinger
equation and compare these results with the strong-field and Coulomb-Volkov
approximations. In order to disentangle short- from long- ranged effects of
the atomic potential we also perform simulations for atomic model potentials
featuring a Yukawa-type short-range potential. We find significant deviations
of the ab-initio phase delays between different photoelectron pathways from
the predictions by the strong-field approximation even at energies well above
the ionization threshold. We identify similarities but also profound
differences to the well-known interferometric extraction of phase- and time
delays in one-photon ionization.

###### pacs:

32.80.Rm,32.80.Fb,03.65.Sq

##  I Introduction

Measuring and analyzing ionization phases and timing information on electron
wavepackets ionized by absorption of an XUV photon represents one of the major
advances attosecond pulses and phase-controlled femtosecond laser pulses have
enabled during the last decade [2, 1]. Such XUV pulses in combination with
near-infrared or visible (NIR/V) laser light permit the control of electronic
motion on the shortest accessible timescales [3, 4, 5, 6]. Pump-probe
techniques such as attosecond streaking [7, 8, 9] and reconstruction of
attosecond harmonic beating by interference of two-photon transitions (RABBIT)
[10, 11] have been employed to measure attosecond time-resolved electron
emission from noble gas atoms [12, 13, 14, 15, 16, 17], molecules [18, 19],
and solids [20, 21, 22]. Whereas attosecond streaking of electrons ionized by
an XUV pulse can be understood in terms of a classical time-resolved shift in
momentum and energy by the probing IR field [7, 23, 24, 25, 26, 27], RABBIT
employs two interfering quantum paths to the same final state in the continuum
called sideband [13, 27]. This sideband energy can be reached through a two-
photon process involving absorption of photons from one of two adjacent
harmonic orders of a high-order harmonic generation (HHG) radiation followed
by absorption or emission of an IR photon of the fundamental driving frequency
$\omega$ [28, 29, 30].

Two-color ($\omega-2\omega$) laser fields with well-controlled relative phases
between both colors have been experimentally and theoretically studied since
the last decade of the last century [31, 32, 33, 34, 35]. Recently, they have
also been employed as alternative tool to extract information on ionization
phases and time delays [36, 37, 38, 39]. One key feature is that the broken
inversion symmetry of the $\omega-2\omega$ field allows for interference
between odd and even partial waves of the outgoing photoelectron which leads
to a $(\theta\leftrightarrow\pi-\theta)$ asymmetry of the emission signal.

Recently, Zipp et al. [40] extended the measurement of ionization phases and
attosecond time delays to the strong-field multiphoton regime, providing new
perspectives on time-resolved strong-field ionization. In this novel
$\omega-2\omega$ interference protocol the role of electron wavepackets
emitted by absorption of a single photon from either one or two subsequent
harmonics in the RABBIT protocol is replaced by adjacent ATI peaks generated
by a strong driving field of frequency $2\omega$. The concomitant weaker
$\omega$ field opens up interfering pathways to side bands in between
neighboring ATI peaks by absorbing or emitting one $\omega$ photon. Measuring
the photoelectron angular distribution as a function of the relative phase
$\phi$ between the $\omega$ and the $2\omega$ fields provides information on
the ATI amplitudes. This interferometric approach to multi-photon ionization
(Fig. 1a) resembles the original RABBIT protocol for the extraction of the
ionization phase in one-photon ionization (Fig. 1b). It promises new insights
into relative phases and, possibly, attosecond-scale timing information of
multi-photon strong-field processes. Somewhat simplified, it can address the
question what additional phase delays incur or how much longer it takes
forming a wavepacket by absorbing $N+1$ rather than by $N$ photons. Some works
based on the strong field approximation were recently reported in this
direction [42, 41]. Indeed, first simulations employing semiclassical
trajectory methods [43, 44, 41] highlighted the role of transient trapping of
the wavepacket for the phase shift of the ATI peaks close to or even below the
threshold. A detailed analysis of the information encoded in the ionization
phases, their dependence on the intensities of the driving ($I_{2\omega}$) and
probing ($I_{\omega}$) field, and on the properties of the atomic potential
appears to be still missing.

Figure 1: Comparison between (a) multi-photon strong-field interference
(MPSFI) and the standard RABBIT protocol (b) for two interfering pathways from
the initial bound state $\left|i\right\rangle$ to final states
$\left|f\right\rangle$ in the continuum, schematically. While RABBIT applies
to two pathways involving ionization by one photon with energies
$\left(2m-1\right)\omega$ and $\left(2m+1\right)\omega$ generated by HHG,
MPSFI involves (at least) two ATI peaks generated by absorbing $N$ or
$\left(N+1\right)$ photons from the strong pump field with frequency
$2\omega$. The final state $\left|f\right\rangle$ is reached in either case by
the absorption ($\omega$) or emission ($-\omega$) of one photon of the weak
probe field. Each arrow denotes a one-photon transition.

As will be shown in the following, multi-photon strong-field interference
(MPSFI, Fig. 1a) substantially differs from the standard RABBIT protocol as a
multitude of pathways with different number of photons and a broad range of
partial waves of the emerging electronic wavepacket contribute. Phase delays
can be extracted by this photoelectron interferometry not only at energies
near the so-called sidebands but also near the ATI peaks. Moreover, in the
strong-field setting phase delays are, unlike in the RABBIT protocol, found to
be remarkably sensitive to the probe field strengths, rendering the separation
of the atomic field and laser field influence on the resulting phase and time
delay more challenging.

In this work, we theoretically investigate the phase delays in the multi-
photon regime accessible by such a $\omega-2\omega$ interference protocol for
two collinearly polarized laser fields. We find strong deviations of the time-
dependent Schrödinger equation (TDSE) results from SFA predictions clearly
indicating that the atomic potential has a crucial influence on the ionization
phase of ATI peaks in this strong-field regime even at energies well above the
ionization threshold. We also present a simplified analytical description of
the MPSFI phase delays and discuss their potential to access timing
information.

In Sec. II we briefly introduce the simulation methods employed. In Sec. III
we present numerical results for quantum path interferences in multi-photon
ionization. An approximate analytical approach to the extraction of the
information on ionization phases, phase delays and time delays from such a
$\omega-2\omega$ protocol as well as numerical results for a model atom with a
short-ranged Yukawa-type atomic binding potential are discussed in Sec. IV.
The comparison with experimental data for argon described by a suitable model
potential [45] in single active electron (SAE) approximation [32, 46, 47] is
presented in Sec. V. Concluding remarks are given in Sec. VI. Atomic units are
used unless stated otherwise.

##  II Methods

We consider a multi-femtosecond laser pulse with frequency $\omega$ and its
second harmonic $2\omega$ with electric field amplitude

$F(t)=f(t)\left[F_{2\omega}\sin(2\omega t+\phi)+F_{\omega}\sin(\omega
t)\right]\hat{\bm{z}},$ (1)

where $f(t)$ is the overall pulse envelope and $\hat{\bm{z}}$ is the
polarization direction of both fields. In the present $\omega-2\omega$
scenario $F_{2\omega}$ is the amplitude of the strong pump field giving rise
to ATI peaks and $F_{\omega}$ is the amplitude of the weak probe field, i.e.
$F_{\omega}\ll F_{2\omega}$. The relative phase $\phi$ between the $\omega$
and $2\omega$ fields is the experimentally accessible knob to control the
interference between different multi-photon pathways. In the following, we
will present results for the integral and angular differential photoelectron
spectra as a function of $\phi$. For the envelope function we choose the form
$f\left(t\right)=\sin^{2}\left(\frac{\pi t}{\tau}\right)$, where $\tau$ is the
pulse duration covering 16 cycles in the strong pump field or eight cycles of
the probe field, i.e., $\tau=16\pi/\omega$.

We solve the time-dependent Schrödinger equation (TDSE) in the single-active
electron (SAE) approximation in the length gauge [32, 46, 47],

$i\frac{\partial\psi(\vec{r},t)}{\partial
t}=\left(\frac{p^{2}}{2}+V_{a}(r)+\vec{r}\cdot\vec{F}(t)\right)\psi(\vec{r},t).$
(2)

In our simulation for argon to be compared with the experiment [40] we employ
as atomic potential $V_{a}$ in Eq. (2) the Muller model potential [45]. In
order to delineate the role of short-ranged and long-ranged potentials we
alternatively use a Yukawa-type atomic potential

$V_{a}\left(r\right)=-\frac{b}{r}e^{-r/a},$ (3)

with charge parameter $b$ and the screening length $a$.

In addition to full solutions of the TDSE, we employ two popular versions of
the distorted-wave Born approximation (DWBA) that allow to account for multi-
photon and strong-field processes, namely the strong-field approximation (SFA)
[48, 49, 50] and the Coulomb-Volkov approximation (CVA) [51]. Accordingly, the
transition amplitude from an initial atomic state
$\left|\phi_{i}(t)\right\rangle$ to a final state
$\left|\varphi_{\vec{k}}\right\rangle$ with asymptotic momentum $\vec{k}$ in
the continuum, i.e.,
$a_{\vec{k}}(\phi)=\lim_{t\rightarrow\infty}\left\langle\varphi_{\vec{k}}\right|\left.\psi(t)\right\rangle$
in the DWBA, is given by

$a(\vec{k},\varphi)=-i\int\limits_{-\infty}^{+\infty}dt\
\langle\chi_{\vec{k}}^{\textsc{DW}}(t)|z\,F\,(t)\left|\phi_{i}(t)\right\rangle.$
(4)

From Eq. (4), the SFA follows when the Volkov state is used as the distorted
wave [48, 49, 50]

$\chi_{\vec{k}}^{(\textsc{DW})-}(\vec{r},t)=\chi_{\vec{k}}^{(\textsc{V})-}(\vec{r},t)=\frac{\exp\mathbf{[}i(\vec{k}+\vec{A})\cdot\vec{r}\mathbf{]}}{\left(2\pi\right)^{3/2}}\exp\left[-i\int_{t}^{+\infty}dt^{\prime}\frac{(\vec{k}+\vec{A}(t^{\prime}))^{2}}{2}\right]\
.$ (5)

The CVA results when approximating the distorted wave by a product of the
Volkov solution and the Coulomb wave [51]

$\chi_{\vec{k}}^{(\textsc{DW})-}(\vec{r},t)=\chi_{\vec{k}}^{(CV)-}(\vec{r},t)=\chi_{\vec{k}}^{(V)-}(\vec{r},t)\;\mathcal{D}_{C}(Z_{T},\vec{k},t),$
(6)

where $\mathcal{D}_{C}(Z_{T},\vec{k},t)=N_{T}^{-}(k)\
_{1}F_{1}(-iZ_{T}/k,1,-ik\ r-i\vec{k}\cdot\vec{r})$ for a hydrogenic atom. The
Coulomb normalization factor $N_{T}^{-}(k)=\exp(\pi
Z_{T}/2k)\Gamma(1+iZ_{T}/k)$ coincides with the amplitude of the Coulomb wave
function at the origin, ${}_{1}F_{1}$ denotes the confluent hypergeometric
function, and $Z_{T}$ is the electric charge of the parent ion. Eq. (5)
describes the final state of a free electron wave in the strong laser field
while completely neglecting the atomic potential. The CVA in Eq. (6) includes
also the Coulomb scattering of the free electron but neglects the effect of
binding and of dynamical Stark shifts. These two DWBA approximations provide
points of reference for identifying dynamical multi-photon effects on
ionization phases.

Because of the azimuthal symmetry, the electron probability distribution
$P(\vec{k})=\left|a_{\vec{k}}\right|^{2}$ depends only on the electron
momentum parallel ($k_{z}$) and transverse ($k_{\perp}$) to the field
polarization direction or, alternatively, on the kinetic energy $E$ and the
polar emission angle $\theta$, i.e.,
$P\left({k_{\perp},k_{z},\phi}\right)=(2E)^{-1/2}P\left({E,\cos\theta,\phi}\right)$.
In the multiphoton regime, the photoelectron spectrum is composed of a series
of peaks positioned at energies $E_{n}$

$E_{n}=n\omega-(I_{p}+U_{p}),$ (7)

corresponding to absorption of a given number of $n_{\omega}$ photons of
frequency $\omega$ and $N$ photons of frequency $2\omega$, such that
$n\omega=n_{\omega}\omega+N\left(2\omega\right)$. In Eq. (7), $I_{p}$ and
$U_{p}$ denote the ionization potential and the ponderomotive energy,
respectively. As a given peak $E_{n}$ can be reached by different combinations
of photon numbers $n_{\omega}$ and $N$, photo-electron interferometry in this
strong-field setting is characterized by multi-path interferences of partial
waves with opposite parity. Consequently, an important quantity for
characterizing interferences between partial waves of opposite parity and,
thus, to map out ionization phases in the $\omega-2\omega$ protocol is the
forward-backward ($\theta\leftrightarrow\pi-\theta$) asymmetry of the
photoelectron emission probability

$A(E,\phi)=\frac{S_{+}(E,\phi)-S_{-}(E,\phi)}{S_{+}(E,\phi)+S_{-}(E,\phi)},$
(8)

where the forward (backward) emission spectra $S_{+}$ ($S_{-}$) are obtained
by integrating the momentum distribution over the +z (-z) hemisphere

$S_{+(-)}(E,\phi)=\int_{0(-1)}^{1(0)}\mathrm{d}\cos\theta\>P(E,\cos\theta,\phi).$
(9)

The calculated (or measured) signal function, generically denoted by
$S(E,\phi)$, representing in the following either the photoemission
probability into one hemisphere, $S_{+(-)}$ [Eq. (9)] or the photoelectron
asymmetry $A(E,\phi)$ [Eq. (8)], can be written in terms of a Fourier series
in the relative phase $\phi$. The emission signal takes the form [52]

$S(E,\phi)=c_{0}(E)+\sum_{i=1}^{\infty}c_{i}(E)\cos(i\phi-\delta_{i}(E)),$
(10)

where the leading term ($i=1$) provides the information of the relative
ionization phase $\delta_{1}(E)=\delta(E)$ in analogy to the RABBIT protocol
[52]. Higher-order Fourier components $c_{i}(E)$ for
$\left\\{i=2,3...\right\\}$ should provide an error estimate of the fit.
Extending the analogy to RABBIT, Zipp et al. [40] introduced an Eisenbud-
Wigner-Smith (EWS) -type time delay [53, 54] by mapping the phase delay
$\delta(E)$ onto a time delay as

$\tau(E)=\frac{\delta(E)}{2\omega}.$ (11)

Eq. (11) can be viewed as finite-difference approximation to the spectral
derivative $d\delta(E)/dE$ of the phase shift $\delta(E)$. We explore the
physical significance of $\delta(E)$ and $\tau(E)$ in more detail below.

##  III Energy dependence of phase delay in multi-photon ionization

As a representative example of $\omega-2\omega$ atomic ionization, we choose
the probe field with the fundamental frequency of a Ti:Sapphire laser of 800
nm wavelength in the near-infrared (NIR) region of the spectrum, and the pump
frequency as its second harmonic with a 400 nm wavelength in the visible (V)
region. In line with the experiment of Zipp et al. [40] we study atomic
ionization of argon by the two-color laser field in Eq. (1) with intensities
$I_{2\omega}=c/(8\pi)F_{2\omega}^{2}=8\times 10^{13}$ W/cm2 and
$I_{\omega}=c/(8\pi)F_{\omega}^{2}=4\times 10^{11}$ W/cm2. In Fig. 2 we
exhibit the results of our TDSE calculations in the SAE approximation [46,
47]. In Fig. 2a we show the variation of the total multiphoton spectrum
(integrated over all emission angles $\theta$) as a function of the relative
two-color phase $\phi$. We observe the typical multiphoton peak structure with
peak positions at energies predicted by Eq. (7), with the ionization potential
for argon of $I_{p}=15.78$ eV and ponderomotive energy $U_{p}=1.19$ eV. ATI
peaks at even multiples of $\omega$ result predominantly from absorption of
$N$ photons of frequency $2\omega$, while peaks at odd multiples of $\omega$
result from absorption or emission of at least one additional $\omega$ (probe)
photon. Following the convention of RABBIT [13, 27] we refer to the latter
group of peaks with energies near odd multiples of the NIR frequency $\omega$
as “sidebands” (SB). Unlike the total electron emission integrated over all
angles $\theta$ (Fig. 2a) whose $\phi$ dependence displays a $\pi$ periodicity
(emission near $\phi$ and $\phi+\pi$ are identical), the emission into the
forward hemisphere $S_{+}(E,\phi)$ given by Eq. (9) ($0\leq\theta\leq\pi/2$)
(Fig. 2b) and the asymmetry parameter $A(E,\phi)$ (Fig. 2c) display a $2\pi$
periodicity indicative of the parity-breaking contributions due to
$\omega-2\omega$ interferences. These structures are magnified in the close-up
figures 2d, 2e, and 2f.

Figure 2: (a) Total photoelectron spectrum (logarithmic scale), (b) forward
emission spectrum integrated in the +z hemisphere (logarithmic scale) and (c)
asymmetry parameter $S$ (linear scale) for argon as a function of the relative
phase $\phi$ and the electron energy (in eV) calculated within the TDSE. The
laser intensities are $I_{2\omega}=8\times 10^{13}$ W/cm2 and
$I_{\omega}=4\times 10^{11}$ W/cm2 for the respective frequencies $2\omega$
and $\omega=0.057$ a.u. with pulse duration $\tau=881.85$ a.u., corresponding
to eight full cycles of the latter. (d-f) Zoom in linear scale corresponding
to (a-c), respectively.

From the fit of the variation of the numerical data for $S_{+}(E,\phi)$ and
$A(E,\phi)$ to the Fourier expansion [Eq. (10)] at fixed $E$, the relative
ionization phase $\delta(E)$ can be extracted as the phase shift of the
$\cos\phi$ oscillation. Because of the broad Fourier width of the ultrashort
pulse [Eq. (1)] the multiphoton electron spectrum (Fig. 2a) is a continuous
function of $E$. Accordingly, also the phase shift $\delta(E)$ can be viewed
as a continuous function of $E$. Results for the energy dependence of
$\delta(E)$ predicted by the TDSE, the SFA, and the CVA calculations of
$S(E,\phi)$ are shown in Fig. 3. Most strikingly, the SFA jumps almost
discontinuously and periodically between $\pi$ near the ATI energies [even $n$
in Eq. (7)] and $0$ in the vicinity of side bands [odd $n$ in Eq. (7)]. The
CVA introduces modest variations to this SFA behavior which are a signature of
Coulomb scattering of the ionized electron. By contrast, the full TDSE
solution displays significant deviations from the SFA predictions indicating a
much more complex variation of the energy dependence interference phase
$\delta(E)$. Even at relatively large energies above the threshold ($\sim 20$
eV), no clear indication for the convergence towards the SFA limit as assumed
in previous analyses [40, 43] emerges. These strong variations of $\delta(E)$
and deviations from SFA appearing in the TDSE results are the signature of
simultaneous interaction of the escaping electron with both the atomic force
field and the strong laser fields, in particular intermediate off-shell bound-
bound and continuum-continuum (cc) transitions between field-dressed atomic
states [55]. Such contributions are absent in the SFA and the CVA.

Figure 3: Continuum phase shifts $\delta(E)$ extracted from the asymmetry
$A(E,\phi)$ as a function of the emission energy from the TDSE (thick black
solid line), SFA (dashed blue line), and CVA (thin red solid line) results.
Thick solid vertical gray lines denote ATI peak energies and dashed vertical
lines sideband energies according to Eq. (7). The horizontal dashed line
corresponds to the strong-field limit for ATI phase shifts [$\delta(E)=\pi$].

Interpretation of the phase shift $\delta(E)$ of the forward (or backward)
emission or asymmetry signal [Eq. (10)] requires a more detailed analysis of
the interfering quantum paths. Key point is that in the present
$\omega-2\omega$ multi-photon strong-field interference (MPSFI) scenario a
multitude of pathways contribute, a few of them shown for argon in Fig. 4,
well beyond the subset invoked in the analogy to the RABBITT protocol (Fig.
1a). This renders a quantitative analysis more challenging. For example, the
side-band energy $E_{n}=E_{15}$ can be reached not only by the path pair
$P_{1}$ (Fig. 4a), which resembles the RABBIT protocol, but also by other path
pairs with different sequences of absorption and emission events to the same
first order in the weak probe field (e.g. $P_{2}$, $P_{3}$,…), or to different
orders in the pump field (e.g. $P_{1}^{\prime},P_{2}^{\prime}$). The (virtual)
intermediate states reached by the probe photon may involve continuum (e.g.
$P_{2},P_{1}^{\prime}$) or bound states (e.g. $P_{3},P_{2}^{\prime}$). The
latter are expected to be more important when the path proceeds via a bound-
state resonance.

Figure 4: Examples of pairs of quantum paths reaching the sideband energy
$E_{n}=E_{15}$ (a) or the ATI (or main) peak energy $E_{n}=E_{14}$ (b) for
argon. In (a), the set $P_{i}(i=1,...)$ features pairs, each with absorption
(left L) or emission (right R) of one weak probe photon $\omega$ and the set
$P_{i}^{\prime}$ features pairs with one additional absorption and emission of
the strong pump photon (right R) compared to the direct path (left L) while
both absorbing one weak probe photon $\omega$. The absorption of the probe
photons may occur in the continuum ($P_{1}$, $P_{2}$ and $P_{1}^{\prime}$) or
in virtual intermediate bound states ($P_{3}$ and $P_{2}^{\prime}$). In (b),
the direct ATI process can, to lowest order, interfere with path pairs $P_{1}$
involving absorption or emission of two $\omega$ photons, $P_{2}$ are examples
of a process involving four $\omega$ photons. $P_{0}^{\prime}$ represents one
contribution to the dressing of the ATI electron by the probe field.

Multi-photon path interferences can be analyzed not only near the sidebands
(Fig. 4a) but also near ATI (or main) peaks (Fig. 4b). For example, at ATI
energy $E_{n}=E_{14}$ the direct path $P_{0}$ from the initial state to the
final state with $E_{14}$ via absorption of $7$ photons with frequency
$2\omega$ can interfere with a multitude of paths involving two probe photons
$P_{1}^{\prime\prime}$ (Fig. 4b), which are of the same order in the weak
field as the dressing of the ATI electron by the IR field $P_{0}^{\prime}$.
For a stronger probe field, even higher-order contributions may become
important; examples of which involve the absorption or emission of 2 or 4
$\omega$ photons are shown in Fig. 4b. It is important to realize that the set
of paths in Fig. 4 still do not fully reflect the complexity of the ensemble
of contributing interfering paths as the angular momentum degree of freedom is
omitted here for simplicity (see [56]). Each additional photon absorption or
emission process leads to a branching of paths to multiply degenerate states
of the same energy $E$ but different angular momenta
$\ell\rightarrow(\ell+1,\ell-1)$. Consequently, for an initial state with
angular momentum $\ell_{i}$ all partial waves $E$ within the interval
$\left[\max(0,\ell_{i}-N),\ell_{i}+N\right]$ can be coherently populated at
the final energy when the pulses are linearly polarized.

##  IV Analytical model for quantum path interferences in multi-photon
ionization

In order to provide an intuitive guide towards interpreting the ionization
phase shift $\delta(E)$ extracted from the quantum path interferences
contributing to MPSFI we present a simplified analysis based on a (lowest-
order) perturbative multi-photon description. Accordingly, the contribution of
the $N$-photon absorption path (e.g., $P_{0}$ in Fig. 4b) to electron emission
in the $\theta$ direction following the absorption of $N$ photons of frequency
$2\omega$ in the visible has the complex amplitude

$C\left(E_{2N},N\right)=\sum_{\ell}A_{N,\ell}\exp\left[i\left(N\phi-N\frac{\pi}{2}-\ell\frac{\pi}{2}+\eta_{\ell}\left(E_{2N},F\right)\right)\right]Y_{\ell}^{0}\left(\theta\right).$
(12)

In Eq. (12) $A_{N,\ell}$ is the modulus of the $N$-photon absorption amplitude
and $\eta_{\ell}\left(E_{2N},F\right)$ is the atomic ionization phase at
energy $E=E_{2N}$. In the weak field-limit this phase is expected to approach
the one-photon atomic ionization phase at the same energy and angular momentum
$\eta_{\ell}\left(E_{2N},F\rightarrow
0\right)=\eta_{\ell}\left(E_{2N}\right)$. However, in the present strong-field
setting, deviations from this limit are expected. The sum in Eq. (12) extends
over all orbital quantum numbers fulfilling the inequality
$\max\left[0,\ell_{i}-N\right]\leq\ell\leq\ell_{i}+N$. For estimating the
phases in Eq. (12) we have used that each photon absorption or emission event
contributes a phase $\pi/2$, each angular momentum change $\Delta\ell$ adds
another $\Delta\ell\pi/2$, and each absorption of a $2\omega$ pump photon
includes an additional relative phase $\phi$ of the pump field relative to the
probe field [see Eq. (1)]. Applying now Eq. (12) to the left (L) path of pair
$P_{1}$ (Fig. 4a) contributing near the sideband energy $E_{2N+1}$, the
combined amplitude for absorbing $N$ visible (V) $2\omega$ photons followed by
absorbing one NIR $\omega$ photon reads

$\displaystyle C_{P_{1},L}\left(E_{2N+1}\right)=$
$\displaystyle\sum_{\ell,\sigma=\pm
1}A_{N,\ell}^{\mathrm{V}}A_{1+,\sigma}^{\mathrm{NIR}}$
$\displaystyle\exp\left\\{i[N\phi-(N+1)\frac{\pi}{2}-(\ell+\sigma)\frac{\pi}{2}+\eta_{\ell}(E_{2N},F)+\varphi_{\ell+\sigma}^{cc,1+}(E_{2N},F)]\right\\}Y_{\ell+\sigma}^{0}(\theta)$

with $\sigma=\Delta\ell=\pm 1$ the change in angular momentum due to the
absorption of an additional NIR photon. $A_{1+,\sigma}^{\mathrm{NIR}}$ denotes
the modulus and $\varphi_{\ell+\sigma}^{cc,1+}(E_{n-1},F)$ the corresponding
additional phase of the absorption of one additional ($1+$) NIR photon. It
describes the continuum-continuum transition to the angular momentum sector
$\ell+\sigma$ in the sideband reached by the absorption of $N$ photons of
frequency $2\omega$ and one additional photon of frequency $\omega$, i.e.,
$n=2N+1$. In the perturbative limit, this phase is the analogue to the
corresponding phase in RABBIT which depends, in general, on $\ell$ [17].
However, for probe fields beyond the perturbative limit, the continuum-
continuum phase is expected to be dependent also on $F_{\omega}$. When both
pump and probe fields are simultaneously present [(Eq. (1)], the phases will
depend, in general, on the combined field $F$. The corresponding expression
for the right (R) of the path pair $P_{1}$ is accordingly given by

$\displaystyle C_{P_{1},R}\left(E_{2N+1}\right)=$
$\displaystyle\sum_{\ell,\sigma=\pm
1}A_{N+1,\ell}^{\mathrm{V}}A_{1-,\sigma}^{\mathrm{NIR}}Y_{\ell+\sigma}^{0}(\theta)$
$\displaystyle\exp\left\\{i[(N+1)\phi-(N+2)\frac{\pi}{2}-(\ell+\sigma)\frac{\pi}{2}+\eta_{\ell}(E_{2(N+1)},F)+\varphi_{\ell+\sigma}^{cc,1-}(E_{2(N+1)},F)]\right\\}$

where $A_{1-,\sigma}^{\mathrm{NIR}}$ denotes the modulus and
$\varphi_{\ell+\sigma}^{cc,1-}(E_{2(N+1)},F)$ the corresponding cc phase of
the emission amplitude of an IR photon . Note that the range of $\ell$
included in Eq. (IV) is different from that in Eq. (IV) and includes
$\max\left[0,\ell_{i}-(N+1)\right]\leq\ell\leq\ell_{i}+N+1$. When, e.g., only
the path pair $P_{1}$ in Fig. 4a is considered, the emission probability near
the sideband $E=E_{2N+1}$ [Eq. (9)] is now given by the coherent sum of Eq.
(IV) and (IV),

$S_{+(-)}(E_{n},\phi)=\int_{0(-1)}^{1(0)}\mathrm{d}\cos\theta\>\left|C_{P_{1},L}\left(E_{n}\right)+C_{P_{1},R}\left(E_{n}\right)\right|^{2}.$
(15)

The evaluation of Eq. (15) can be drastically simplified by including only the
dominant pathways along the so called “yrast line” well known from beam-foil
spectroscopy [57, 58] or, equivalently, assuming that only the pathways
preferred by the Fano propensity rule [59, 60] are realized. Accordingly, each
photoabsorption leads predominantly to an increase
$\left(1+\leftrightarrow\sigma=1\right)$ and photoemission to a decrease
$\left(1-\leftrightarrow\sigma=-1\right)$ by one unit of angular momentum.
Including only these dominant paths eliminates the summation over $\ell$ and
$\sigma$ in Eqs. (IV) and (IV). We note that this approximate selection rule
is only applicable to resonant bound-bound or continuum-continuum transitions
but not to tunneling or above-threshold ionization. For ATI peaks close to
threshold (th), the dominant $\ell$ values are delimited by [61]

$\ell\leq\ell_{\mathrm{th}}\leq\left(2Z_{T}\alpha\gamma\right)^{1/2}=\left(2\sqrt{2}Z_{T}\sqrt{\frac{N_{\mathrm{th}}}{2\omega}}\right)^{1/2}$
(16)

where $\alpha$ is the quiver amplitude, $\gamma$ the Keldysh parameter of the
laser field with frequency $2\omega$, and $N_{\mathrm{th}}$ the minimum number
photons of frequency $2\omega$ required to reach the continuum
($N_{\mathrm{th}}=6$ for argon). Accordingly, our TDSE calculations yield $f$
waves as dominant partial waves near threshold, which is very close to the
upper bound predicted by Eq. (16) $\ell_{\mathrm{th}}=4$ and well below the
prediction for the yrast line (or propensity rule [59, 62])
$\ell_{i}+N_{\mathrm{th}}=7$ as depicted in Fig. 5. The partial wave content
of the first ATI peak above threshold and starting point of the further spread
in angular momentum is thus centered at lower values of
$\ell\leq\ell_{\mathrm{th}}$. The evolution of the partial wave distribution
$p_{\ell}$ to higher partial waves with increasing ATI peak is discernible
(Fig. 5c). The first ATI peak exhibits a dominant angular momentum of
$\ell_{\mathrm{th}}=3$, whereas for the second ATI peak the dominant angular
momentum is $\ell=4$. The combined contribution of the $d$ and $g$ waves of
the second ATI peak produces a dominant $f$ wave ($\ell=3$) for the third ATI
peak but with an appreciable $\ell=5$ contribution, i.e., $p_{5}\simeq
0.5p_{3}$. Applying the approximate propensity rule to Eqs. (IV), (IV), and
(15) yields, e.g.,

$\displaystyle S_{+}(E_{2N+1},\phi)$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\mathrm{d}\cos\theta\left\\{(A_{N,\ell}^{\mathrm{V}})^{2}(A_{1+}^{\mathrm{NIR}})^{2}(Y_{\ell+1}^{0}(\theta))^{2}+(A_{N+1,\ell+1}^{\mathrm{V}})^{2}(A_{1-}^{\mathrm{NIR}})^{2}(Y_{\ell}^{0}(\theta))^{2}\right.$
$\displaystyle+$ $\displaystyle
2A_{N,\ell}^{\mathrm{V}}A_{N+1,\ell+1}^{\mathrm{V}}A_{1+}^{\mathrm{NIR}}A_{1-}^{\mathrm{NIR}}Y_{\ell+1}^{0}(\theta)Y_{\ell}^{0}(\theta)$
$\displaystyle\times$
$\displaystyle\left.\cos\left[\phi+\eta_{\ell+1}(E_{2(N+1)},F)-\eta_{\ell}(E_{2N},F)+\varphi^{cc,1-}_{\ell}(E_{2(N+1)},F)-\varphi^{cc,1+}_{\ell+1}(E_{2N},F)\right]\right\\}.$

with an analogous expression for $S_{-}\left(E,\phi\right)$. Consequently, the
asymmetry $A(E=E_{2N+1},\phi)$ given by Eq. (8) is proportional to

$\displaystyle A(E_{2N+1},\phi)$ $\displaystyle\sim$ $\displaystyle
S_{+}(E_{2N+1},\phi)-S_{-}(E_{2N+1},\phi)$ $\displaystyle\sim$ $\displaystyle
2A_{N,\ell}^{\mathrm{V}}A_{N+1,\ell+1}^{\mathrm{V}}A_{1+}^{\mathrm{NIR}}A_{1-}^{\mathrm{NIR}}\int_{0}^{1}\mathrm{d}\cos\theta
Y_{\ell+1}^{0}(\theta)Y_{\ell}^{0}(\theta)$ $\displaystyle\times$
$\displaystyle\cos\left[\phi+\eta_{\ell+1}(E_{2(N+1)},F)-\eta_{\ell}(E_{2N},F)+\varphi^{cc,1-}_{\ell}(E_{2(N+1)},F)-\varphi^{cc,1+}_{\ell+1}(E_{2N},F)\right].$

Comparison with Eq. (10) yields now an explicit analytic but approximate
expression of the phase delay between the two paths of the pair $P_{1}$ (Fig.
4a)

$\delta(E_{2N+1})\simeq\eta_{\ell}(E_{2N},F)-\eta_{\ell+1}(E_{2(N+1)},F)+\varphi_{\ell+1}^{cc,1+}(E_{2N},F)-\varphi_{\ell}^{cc,1-}(E_{2(N+1)},F).$
(19)

In the limit where all contributions to the phase of the wavepacket due to the
interplay with the atomic force field and the laser field can be neglected,
$\delta(E_{2N+1})\approx 0$, Eq. (IV) reduces to

$A(E_{2N+1},\phi)\propto S_{+}(E_{2N+1},\phi)-S_{-}(E_{2N+1},\phi)=C\cos\phi$
(20)

which agrees with the result in the SFA approximation first given by Zipp et
al. [40].

Figure 5:  Energy spectrum and angular momentum distribution after strong-
field ionization of argon by the one-color $2\omega$ field with the same
parameters as in Fig. 2. (a) Photoelectron spectrum, (b) electron distribution
as a function of the energy and angular momentum on a logarithmic scale
covering three orders of magnitude, and (c) normalized $p_{\ell}$ (integrated
over energy) for the first three ATI peaks from threshold.

A similar analysis for a pair of paths contributing to the asymmetry near the
ATI energy $E_{n}$, where now $n=2N$, taking into account only interference
between the direct ATI path $P_{0}$ and the path $P_{1}^{\prime\prime}(R)$
(Fig. 4b) involving absorption of two NIR photons yields

$\displaystyle A(E_{2N},\phi)$ $\displaystyle\sim$ $\displaystyle
S_{+}(E_{2N},\phi)-S_{-}(E_{2N},\phi)$ $\displaystyle\sim$ $\displaystyle
A_{N,\ell}^{\mathrm{V}}A_{N-1,\ell-1}^{\mathrm{V}}A_{2+}^{\mathrm{NIR}}\int_{0}^{1}\mathrm{d}\cos\theta
Y_{\ell+1}^{0}(\theta)Y_{\ell}^{0}(\theta)$ $\displaystyle\times$
$\displaystyle\cos\left[\phi+\pi+\eta_{\ell}(E_{2N},F)-\eta_{\ell-1}(E_{2(N-1)},F)-\varphi_{\ell+1}^{cc,2+}(E_{2(N-1)},F)\right].$

with $A_{2+}^{\mathrm{NIR}}$ $\left(\varphi_{\ell+1}^{cc,2+}\right)$ the
modulus (phase) of the two-photon transition amplitude from the ATI peak at
$E_{n-2}$ with $\ell-1$ to $\left(E_{n},\ell+1\right)$. Consequently, the
phase delay between these two paths $\delta(E)$ is given by

$\delta(E_{2N})\simeq-\pi-\eta_{\ell}(E_{2N},F)+\eta_{\ell-1}(E_{2(N-1)},F)+\varphi_{\ell+1}^{cc,2+}(E_{2(N-1)},F).$
(22)

In the limit that all atomic force field and laser field effects on the phase
delay can be neglected, the SFA limit would emerge as

$S_{+}(E_{2N},\phi)-S_{-}(E_{2N},\phi)=A\cos\left(\phi+\pi\right),$ (23)

which results, indeed, in a phase jump of $\pi$ between the sidebands [Eq.
(20)] and the ATI peaks [Eq. (23)] in agreement with our numerical results
(Fig. 3). Consequently, the deviations observed in the TDSE simulation and CVA
simulations from these SFA limit are an unambiguous signature of the interplay
between the atomic force field and laser fields in the atomic ionization
phases. It should be emphasized that the TDSE results include all paths
contributing to the multi-photon strong field interference for photoelectron
well beyond the simple “two-path double-slit” model [Eq. (IV) and (23)]
explicitly treated above.

The two-path model can provide guidance as to which information can be
extracted from MPSFI spectra. For example, the phase contributions $\eta$ and
$\phi^{cc}$ will be, in general, dependent on the field strengths
$F_{2\omega}$ and $F_{\omega}$ in a strong-field $\omega-2\omega$ scenario
fundamentally different from the standard RABBIT protocol. Moreover, while the
resulting phase delay $\delta(E)$ is a continuous function of $E$ (see Fig. 3)
the mapping of a phase delay onto a time delay according to Eq. (11) depends
on the specific position within the spectrum. Near sideband energies
$E_{2N+1}$, Eq. (19) has the appearance of a finite difference approximation
as implied by Eq. (11) and can thus be used to extract approximate time delays
$\tau=\delta(E_{2N+1})/2\omega$. Near ATI peaks [Eq. (22)], such
interpretation in terms of a finite-difference approximation fails as the
difference involves now different interfering zero- and two-IR photon paths.
Moreover, when all path pairs are included, a sum over many path pairs each of
which giving rise to terms of the form [Eq. (IV)] for sidebands and of the
form [Eq. (IV)] for ATI peaks will contribute to $A(E,\phi)$ rendering the
extraction of a spectral derivative for a specific phase difficult. Only in
cases where one path pair strongly dominates, in particular the pair $P_{1}$
for the sideband, approximate EWS time delays for a given partial wave can be
unambiguously assigned. With this caveat in place, we also give $\tau(E)$ in
Figs. 6, 7, and 8 for illustrative purposes.

Figure 6: (a) Electron spectra for a Yukawa potential [Eq. (3)] with $a=4$ and
$b=0.629$ calculated for one-color 2$\omega$ (black line) and two-color
$\omega-2\omega$ (red line) laser fields with $\phi=0$. (b) Phase delays
$\delta(E)$ in units of $\pi$ calculated from the asymmetry $A(E,\phi)$
integrated over hemispheres [see Eq. (9)]. For reference we also convert the
phase delay into a time delay [Eq. (11)] (right side axis). The laser
intensities are $I_{2\omega}=10^{11}$ W/cm2 and $I_{\omega}=5\times 10^{8}$
W/cm2. Other laser parameters are the same as in Fig. 2.

Before comparing simulations with experimental data, we illustrate the
partial-wave path-interference structure for a strongly simplified model
system in which the number of contributing paths and, thus, the complexity of
the ionizing process is drastically reduced. We consider an electron bound by
a Yukawa potential [Eq. (3)] with parameters ($a=4$, $b=0.629$) chosen such
that a single $2\omega$ photon is sufficient to reach the continuum and the
shallow potential supports only one 1s-like bound state with $E_{1s}=-0.08$.
Consequently, the energetic position of the first ATI peak coincides in this
case with the position of the standard photoionization peak. For later
reference we note that the screening length of this potential ($a=4$) is
sufficiently large as to include, despite being asymptotically short-ranged,
some Coulomb-laser coupling (CLC) or cc phase contributions [63]. Moreover, we
choose the intensities of the fields sufficiently low ($I_{2\omega}=10^{11}$
W/cm2, $I_{\omega}=5\times 10^{8}$ W/cm2) to be strictly in the perturbative
multi-photon regime. The photoelectron spectrum in both the presence and
absence of the weak probe field are displayed in Fig. 6a. Turning on the
$\omega$ field creates the side bands, as expected, while the ATI peaks remain
largely unaffected by the probe field. The absorption of a single V
($2\omega$) photon from the bound $1s$ initial state ionizes the model atom
creating a $p-$wave electron of energy corresponding to the first peak (ATI1).
The second peak (ATI2) results from the absorption of two-V($2\omega$)
photons, and is composed of the superposition of s and d waves due to the
selection rule of angular momentum $\Delta\ell=\pm 1$. We have determined the
angular momentum composition of ATI2 to contain $9.8\%$ of s character and
$90.1\%$ of d character consistent with the propensity rule invoked above. The
lowest sideband SB1 between the first $2\omega$ photoionization peak ATI1 and
the second peak ATI2 can be reached by either absorption of two photons [one
V:($2\omega$) and one NIR:($\omega$)] or absorption of two V ($2\omega$)
photons and emission of one NIR ($\omega$) photon. For the first sideband SB1
the angular momentum composition is given by $9.4\%$, $0.8\%$, and $89.8\%$
for the $s$, $p$, and $d$ states, respectively. The population of $s-$ and
$d-$partial waves in SB1 is close to that of ATI2 also in line with the
propensity for two-photon absorption irrespective of the different frequencies
involved. This distribution indicates the dominance of the
one-V($2\omega$)-one-NIR($1\omega$) absorption path to the SB1 over the
two-V($2\omega$) absorption and one-NIR($1\omega$) emission path in the
perturbative regime, which is expected since the latter path involves one more
photon from a weak field than the former and, consequently, is a higher-order
photoionization process. However, the latter path provides a small but crucial
contribution giving rise to a non-vanishing $\phi$ dependent contribution from
which the phase delay $\delta(E)$ can be extracted (Fig. 6b).

Remarkably, whereas $\delta(E)$ near the ATI peaks closely follows the SFA
predictions $\delta(E_{n})\simeq\pi$ [Eq. (23)], near the sideband peaks
strong deviations can be observed in Fig. 6b. For the first sidebands for
which this phase could be reliably extracted we find $\delta(E_{n})\simeq
0.3\pi$. For reference we also convert the phase delay $\delta(E)$ into an
EWS-type time delay following Eq. (11) and find for the sideband, within a
fairly small energy window ($3$ eV $\leq E\leq$ $10$ eV), an almost energy-
independent time delay of about $\tau\approx 200$ attoseconds. Using the
approximate expressions [Eqs. (20) and (23)] for a qualitative analysis of the
two-path interference these results suggest that the phase delay near the ATI
peaks is strongly dominated by the SFA contribution ($\sim\pi$) corresponding
to a time delay of $660$ as while atomic field corrections play only a minor
role. By contrast, near the sideband peaks the phase differences induced by
the atomic-field
$\eta_{\ell+1}(E_{2(N+1)})-\eta_{\ell}(E_{2N})+\varphi_{\ell}^{cc,1-}-\varphi_{\ell+1}^{cc,1+}$
are clearly visible. We note that the presence of a non-vanishing contribution
to the phase delay by the one-photon continuum-continuum transition
$\varphi^{cc,1\pm}$ for the Yukawa potential is consistent with the fact that
with increasing screening length ($a=4$ in the present case) an increasing
part of the full long-range Coulomb-laser coupling is restored [63].
Therefore, we can use Eq. (19) to estimate this contribution to the sideband
phase delay as

$\varphi_{\ell+1}^{cc,1+}(E_{2N})-\varphi_{\ell}^{cc,1-}(E_{2(N+1)})\simeq\delta(E_{2N+1})+\eta_{\ell+1}(E_{2(N+1)})-\eta_{\ell}(E_{2N}),$
(24)

where we have dropped the label $F$ because we consider the perturbative limit
($F\rightarrow 0$). The atomic ionization phases $\eta_{\ell}$ can be obtained
by the one-photon atomic ionization phase in a partial-wave expansion for the
Yukawa potential. By using Eq. (24), we estimate the cc phase contribution to
SB1 as $\varphi_{2}^{cc,1+}-\varphi_{1}^{cc,1-}\simeq 0.45$, for SB2 as
$\varphi_{3}^{cc,1+}-\varphi_{2}^{cc,1-}\simeq 0.8$, and for SB3 as
$\varphi_{4}^{cc,1+}-\varphi_{3}^{cc,1-}\simeq 0.92$, corresponding to time
delay contributions of approximately $11$, $19$, and $22$ as, respectively.
These phase contributions could shed some light on how the Yukawa potential
affects the cc contributions to the time delays. Besides, new studies on the
holographic angular streaking of electrons by corotating ($\omega-2\omega$)
fields suggest that nonadiabatic effects in the ionization process could be
responsible for such difference of the time delay with respect to the strong-
field approximation [64, 65]. The identification of non-adiabatic effects on
time delays (included in the TDSE calculations) are beyond the scope of this
paper. It is worth to mention that as the De Broglie’s wavelength of the
electron is longer than the screening length of the Yukawa short-range
potential, classical or semiclassical simulations are not valid for the energy
region shown in Fig. 6b.

##  V Comparison with experiment

Figure 7:  TDSE phase delays $\delta(E)$ calculated as a function of the
emission energy for (a) ATI peaks and (b) sidebands for the same pulse
parameters as in Fig. 2. Phase shifts extracted from data for forwards half
spheres $S_{+}(E,\phi)$ (squares), and asymmetry $A(E,\phi)$ (circles) with
integration over the energy window around each peak energy (full symbols) and
at the energy peak only (open symbols). Full green dots correspond to
experimental data by Zipp et al. [40] normalized to the TDSE result at the
highest sideband energy ($\sim 15$eV).

For a comparison with the experiment of Zipp et al. [40] we extract the multi-
photon ionization interference phase shifts $\delta(E)$ from the TDSE
simulation (Fig. 2). In view of the rapid variation with the energy $E$ (Figs.
2 and 3), we evaluate $\delta(E)$ not only at the ATI or sideband peaks
$E=E_{n}$ [Eq. (7)] but integrate the spectrum over an energy window of width
$\Delta E=0.3\omega$ centered around the peak. We show in Fig. 7 fits to
$\delta(E)$ for emission into forward hemisphere $S_{+}(E,\phi)$ [Eq. (9)] and
for the asymmetry $A(E,\phi)$ [Eq. (8)]. While minor differences of the order
of less than $0.05\pi$ between the different read-outs of $\delta(E)$ (via
$S_{+}$ or $A$) appear, the overall trends observed are independent of the
particular read-out protocol demonstrating that unambiguous information on the
phase delay can be extracted.

For further analysis and interpretation of the results of Fig. 7, two key
points should be taken into account. First, the experimental data for
$\delta(E)$ presented in [40] were relative and set to coincide with the SFA
value ($\delta=0$) at the highest energy measured ($E=15$eV) (a similar
renormalization was used in [43]). However, we observe significant deviations
in $\delta(E)$ from the SFA limit. Therefore, we instead renormalize the
experimental data to the full TDSE result at the highest experimental energy
in order to preserve this additional information on the absolute value of
$\delta(E)$. Accordingly, in Figs. 7a and b the experimental results are set
to coincide with the TDSE phase shifts calculated by integration over the
energy windows around the peaks and all angles in the forward hemisphere.
Overall, the trend in the experimental data is well reproduced by the
simulations. The sharp rise of the phase shift $\delta(E)$ for the first ATI
peak seen close to threshold in both the experiment and simulations was
recently interpreted in terms of transient trapping of the electron in Rydberg
states by the $\omega-2\omega$ field [43].

The second key feature is that the data in Fig. 7 were extracted at a
moderately strong NIR probe field with $I_{\omega}=4\times 10^{11}$ W/cm2. For
the standard RABBIT protocol or attosecond streaking field strengths
$F_{\omega}$ of that order of magnitude were found to be weak enough to
unambiguously extract atomic continuum-continuum or Coulomb-laser coupling
delays which are independent of the particular value of $I_{\omega}$ in line
with lowest-order perturbation theory [1]. However, in the present MPSFI
scenario the influence of the probe field $F_{\omega}$ beyond a lowest-order
perturbation theory must be considered.

Figure 8: Interference phase delay $\delta(E)$ as a function of the probe
laser intensity $I_{\omega}$ extracted from asymmetry parameter integrated
over hemispheres for three ATI peaks and three sidebands with energies as
indicated. All other laser parameters are the same as in Fig. 2. The
horizontal dashed line corresponds to the strong-field limit for ATI phase
shifts [$\delta(E)=\pi$], the SFA limit for the sidebands is $\delta(E)=0$
(not shown).

Indeed, exploring the variation of the extracted $\delta(E)$ at fixed pump
intensity $I_{2\omega}$ as a function of the probe intensity $I_{\omega}$
(Fig. 8) reveals a surprisingly strong dependence. The experimental value
$I_{\omega}=4\times 10^{11}$ W/cm2 is obviously well beyond the lowest-order
perturbative regime which precludes the direct applicability of a RABBIT-type
analysis. For sideband peaks, phase shifts $\delta(E)$ appears to converge to
the perturbative field-independent limit only for considerably lower fields
$I_{\omega}\lesssim 10^{10}$ W/cm2. These converged values differ, however,
significantly from the SFA limit even at the highest energy measured ($E=15.5$
eV). Near ATI peaks, variations are present even at such low intensities and
the approach to converged field-independent values is not yet obvious. It
appears that for the highest energies measured, e.g. $E=17.1$ eV and at the
lowest probe field $I_{\omega}\lesssim 10^{10}$ W/cm2 the phase near the ATI
peak may approach the SFA limit $\delta(E)\simeq\pi$. It should be noted,
however, that the interference contributions to ATI peaks, which are
responsible for the phase shift $\delta(E)$, result from (at least) a two-
photon absorption or emission event in the probe field ($P_{1}^{\prime\prime}$
as depicted in Fig. 4b), which becomes very weak at low $I_{\omega}$ rendering
the phase extraction uncertain. The non-negligible probe field dependence of
the extracted MPSFI phase delays $\delta(E)$, also indicated in Eqs. (19) and
(22) emerges as an important new feature, absent in standard RABBIT or
streaking measurements, that remains to be explored, experimentally as well as
theoretically.

##  VI Concluding remarks

We have presented simulations and the first detailed analysis of the phase
delays $\delta(E)$ in multi-photon ionization. They provide information on the
differences in ionization phases among different pathways open in a
$\omega-2\omega$ scenario for atomic ionization. We show that $\delta(E)$ is
determined by quantum path interferences between different sequences of photon
absorption and emission events. In the SFA limit these phases are given by
$\delta(E)=0$ at sideband energies and by $\delta(E)=\pi$ at the ATI peaks. We
find that the solutions of the time-dependent Schrödinger equation predict
phases strongly differing from these SFA limits even at relatively high
electron emission energies. We relate these phase shifts to the interplay
between the strong $\omega-2\omega$ field and the atomic force field not
accounted for by the SFA. We also point out the intrinsic difficulties to
relate the phase delays $\delta(E)$ to time delays in analogy to the standard
RABBIT protocol for one-photon ionization. A multitude of different
interfering pathways provides obstacles for a straightforward extraction of a
spectral derivative of the phase delay. We have found strong variation of
$\delta(E)$ with the intensities of the pump and probe fields. Our analysis
shows that further experimental insight into the multi-photon ionization phase
delay $\delta(E)$ can be gained by exploring its variation with both
$I_{2\omega}$ and $I_{\omega}$.

###### Acknowledgements.

This work was supported by CONICET PIP0386, PICT-2016-0296 PICT-2017-2945 and
PICT-2016-3029 of ANPCyT (Argentina), Austria-Argentina collaboration
AU/12/02, by the FWF special research programs SFB-041 (ViCoM), and doctoral
programme DK-W1243 (Solid4Fun), and by the European COST Action CA1822. The
computational results presented have been achieved using the Vienna Scientific
Cluster (VSC).

## References

  * [1] R. Pazourek, S. Nagele, and J. Burgdörfer, Rev. Mod. Phys. 87, 765 (2015).
  * [2] F. Krausz and M. Ivanov, Rev. Mod. Phys. 81, 163-234 (2009).
  * [3] V. Véniard, R. Taïeb, and A Maquet, Phys. Rev. Lett. 74, 4161 (1995).
  * [4] J. M. Schins, P. Breger, P. Agostini, R. C. Constantinescu, H. G. Muller, A. Bouhal, G. Grillon, A. Antonetti, and A Mysyrowicz, J. Opt. Soc. Am. B 13, 197 (1996).
  * [5] T. E. Glover, R. W. Schoenlein, A. H. Chin, and C. V. Shank, Phys. Rev. Lett. 76 2468 (1996).
  * [6] J. Hummert, M. Kubin, S. D. López, J. I. Fuks, F. Morales, M. J. J. Vrakking, O. Kornilov, and D. G. Arbó, J. Phys. B: At. Mol. Opt. Phys. 53, 154003 (2020).
  * [7] J. Itatani, F. Quéré, G. L. Yudin, M. Yu. Ivanov, F. Krausz and P. B. Corkum, Phys. Rev. Lett., 88, 173903 (2002).
  * [8] E. Goulielmakis et al., Science 305, 1267 (2004).
  * [9] E. Goulielmakis et al., Science 320, 1614 (2008).
  * [10] V. Véniard, R. Taïeb, and A. Maquet, Phys. Rev. A54, 721 (1996).
  * [11] P. M. Paul, E. S. Toma, P. Breger, G. Mullot, F. Augé, Ph. Balcou, H. G. Muller, and P. Agostini, Science 292, 1689 (2001).
  * [12] M. Schultze, M. Fieß, N. Karpowicz, J. Gagnon, M. Korbman, M. Hofstetter, S. Neppl, A. L. Cavalieri, Y. Komninos, Th. Mercouris, C. A. Nicolaides, R. Pazourek, S. Nagele, J. Feist, J. Burgdörfer, A. M. Azzeer, R. Ernstorfer, R. Kienberger, U. Kleineberg, E. Goulielmakis, F. Krausz, and V. S. Yakovlev, Science 328, 1658 (2010).
  * [13] K. Klünder, J. M. Dahlström, M. Gisselbrecht, T. Fordell, M. Swoboda, D. Guénot, P. Johnsson, J. Caillat, J. Mauritsson, A. Maquet, R. Taïeb, and A. L’Huillier, Phys. Rev. Lett. 106, 143002 (2011).
  * [14] D. Guénot, K. Klünder, C. L. Arnold, D. Kroon, J. M. Dahlström , M. Miranda, T. Fordell, M. Gisselbrecht, P. Johnsson, J. Mauritsson, E. Lindroth, A. Maquet, R. Taïeb, A. L’Huillier, and A. S. Kheifets, Phys. Rev. A, 85, 053424, (2012).
  * [15] D. Guénot, D. Kroon, E. Balogh, E. W. Larsen, M. Kotur, M. Miranda, T. Fordell, P. Johnsson, J. Mauritsson, M. Gisselbrecht, K. Varjú, C. L. Arnold, T. Carette, A. S. Kheifets, E. Lindroth, A. L’Huillier, and J. M. Dahlström, Journal of Physics B: Atomic, Molecular, Optical Physics 47, 245602, (2014).
  * [16] J.M. Dahlström, D. Guénot, K. Klunder, M. Gisselbrecht, J. Mauritsson, A. L’Huillier, A. Maquet, and R. Taïeb, Chemical Physics 414, 53-64 (2013).
  * [17] Jaco Fuchs, Nicolas Douguet, Stefan Donsa, Fernando Martin, Joachim Burgdörfer, Luca Argenti, Laura Cattaneo, and Ursula Keller, Optica 2, 154 (2020).
  * [18] M. Huppert, I. Jordan, D. Baykusheva, A. von Conta, and H. J. Wörner et al., Phys. Rev. Lett 117, 093001 (2016).
  * [19] Beaulieu et al., Science 358, 1288 (2017).
  * [20] A. Cavalieri et al., Nature 449, 1029 (2007).
  * [21] C. Lemell, S. Neppl, G. Wachter, K. Tokesi, R. Ernstorfer, P. Feulner, R. Kienberger, and J. Burgdörfer, Phys. Rev. B 91, 241101(R) (2015).
  * [22] S. Haessler, T. Balciunas, G. Fan, T. Witting, R. Squibb, L. Chipperfield, A. Zaïr, G. Andriukaitis, A. Pugzlys, J. W. G. Tisch, J. P. Marangos, and A. Baltuska, Ultrafast Phenomena XIX, Springer Proceedings in Physics, Volume 162, p. 72. Springer International Publishing Switzerland, 2015.
  * [23] S. Nagele, R. Pazourek, J. Feist, and J. Burgdörfer, Phys. Rev. A 85, 033401 (2012).
  * [24] R. Pazourek, J. Feist, S. Nagele, and J. Burgdörfer, Phys. Rev. Lett. 108, 163001 (2012).
  * [25] R. Della Picca, A. A. Gramajo, S. D. López, and D. G. Arbó, Journal of Physics: Conference Series 1412, 042002 (2020).
  * [26] R. Della Picca, M. F. Ciappina, M. Lewenstein, and D. G. Arbó, Phys. Rev. A102, 043106 (2020).
  * [27] J. M. Dahlström, A. L’Huillier, and A. Maquet, J. of Phys. B: Atomic, Molecular and Optical Physics 45, 183001, (2012).
  * [28] A. S. Kheifets, Phys. Rev. A 87, 063404 (2013).
  * [29] J. Feist, O. Zatsarinny, S. Nagele, R. Pazourek, J. Burgdörfer, X. Guan, K. Bartschat and B. I. Schneider, Phys. Rev. A 89, 033417, (2014).
  * [30] Jing Su, Hongcheng Ni, Andreas Becker, and Agnieszka Jaron-Becker, Phys. Rev. A87, 033420 (2013).
  * [31] D. W. Schumacher, F. Weihe, H. G. Muller, and P. H. Bucksbaum, Phys. Rev. Lett. 73, 1344 (1994).
  * [32] D. G. Arbó, C. Lemell, S. Nagele, N. Camus, L. Fechner, A. Krupp, T. Pfeifer, S. D. López, R. Moshammer, and J. Burgdörfer Phys. Rev. A 92, 023402 (2015).
  * [33] F. Ehlotzky, Phys. Rep. 345, 175 (2001).
  * [34] X. Xie, S. Roither, D. Kartashov, E. Persson, D. G. Arbó, L. Zhang, S. Gräfe, M. S. Schöffler, J. Burgdörfer, A. Baltuška, and M. Kitzler, Phys. Rev. Lett., 108, 193004 (2012).
  * [35] D. G. Arbó, S. Nagele, X.-M. Tong, X. Xie, M. Kitzler, and J. Burgdörfer, Phys. Rev. A, 89, 043414 (2014).
  * [36] D. You, K. Ueda, E. V. Gryzlova, A. N. Grum-Grzhimailo, M. M. Popova, E. I. Staroselskaya, O. Tugs, Y. Orimo, T. Sato, K. L. Ishikawa, et al., Phys. Rev. X 10, 031070 (2020).
  * [37] J. Fuchs et al., Phys. Rev. Lett. (2020, submitted) (https://arxiv.org/abs/2012.07426).
  * [38] S. Donsa, N. Douguet, J. Burgdörfer, I. Brezinová, and L. Argenti, Phys. Rev. Lett. 123, 133203 (2019).
  * [39] G. Laurent, W. Cao, H. Li, Z. Wang, I. Ben-Itzhak, and C. L. Cocke, Phys. Rev. Lett. 109, 083001 (2012).
  * [40] L. J. Zipp, A. Natan, and P. H. Bucksbaum, Optica 1, 361-364 (2014).
  * [41] Yudi Feng, Min Li, Siqiang Luo, Kun Liu, Baojie Du, Yueming Zhou, and Peixiang Lu, Phys. Rev. A 100, 063411 (2019).
  * [42] M. Bertolino and J.M. Dahlström. Phys. Rev. Res 3, 013270 (2021).
  * [43] Xiaohong Song, Guangluo Shi, Guojun Zhang, Jingwen Xu, Cheng Lin, Jing Chen, and Weifeng Yang, Phys. Rev. Lett. 121, 103201 (2018).
  * [44] Anatoli S. Kheifets and Alexander W. Bray, Phys. Rev. A103, L011101 (2021).
  * [45] H. G. Muller, Phys. Rev. A 60, 1341 (1999).
  * [46] X.-M. Tong and Shih.-I. Chu, Chemical Physics, 217, 119 (1997).
  * [47] X.-M. Tong and Shih.-I. Chu, Phys. Rev. A61, 031401(R) (2000).
  * [48] V. Keldysh, Zh. Eksp. Theo. Fiz. 47, 1945 (1964); Sov. Phys. JETP 20, 1307 (1965).
  * [49] F. H. M. Faisal, J. Phys. B 6, L89 (1973).
  * [50] H. R. Reiss, Phys. Rev. A 22, 1786 (1980).
  * [51] M. Jain and N. Tzoar, Phys. Rev. A 18, 538 (1978).
  * [52] Stefan Donsa, Manuel Ederer, Renate Pazourek, Joachim Burgdörfer, and Iva Brezinová, Phys. Rev. A102, 033112 (2020).
  * [53] E. P. Wigner, Phys. Rev. 98, 145 (1955).
  * [54] F. T. Smith, Phys. Rev. 118, 349 (1960); erratum Phys. Rev. 119, 2098 (1960).
  * [55] C.-H. Zhang and U. Thumm, Phys. Rev. A82, 043405 (2010).
  * [56] Divya Bharti, David Atri-Schuller, Gavin Menning, Kathryn R. Hamilton, Robert Moshammer, Thomas Pfeifer, Nicolas Douguet, Klaus Bartschat, and Anne Harth, Phys. Rev. A103, 022834 (2021).
  * [57] F. Bell, G. Trollmann, H. Böckl, H.-D. Betz, Nuclear Instruments and Methods in Physics Research, 194, 423-427 (1982).
  * [58] F. Bell, G. Trollmann, H. D. Betz, Phys. Lett. A 88, 37-39 (1982).
  * [59] U. Fano, Phys. Rev A 32, 617 (1985).
  * [60] D. Busto, J. Vinbladh, S. Zhong, M. Isinger, S. Nandi, S. Maclot, et al., Phys. Rev. Lett. 123, 133201 (2019).
  * [61] D. G. Arbó, C. Lemell, J., Burgdörfer, J. Phys. Conf. Ser. 635, 012003 (2015).
  * [62] Mattias Bertolino, David Busto, Felipe Zapata, and Jan Marcus Dahlström, J. Phys. B 53, 144002 (2020).
  * [63] S. Nagele, R. Pazourek, M. Wais, G. Wachter, J., Burgdörfer, J. Phys. Conf. Ser. 488, 012004 (2014).
  * [64] S. Eckart, Phys. Rev. Res. 2, 033248 (2020).
  * [65] D. Trabert, S. Brennecke, K. Fehre, N. Anders, A. Geyer, S. Grundmann, M. S. Schöffler, L. Ph. H. Schmidt, T. Jahnke, R. Dr̈ner, M. Kunitski, and S. Eckart, Nat. Commun. 12, 1697 (2021).


# On a Fast Solution Strategy for a Surface-Wire Integral Formulation of the
Anisotropic Forward Problem in Electroencephalography

Carlo Baronio(1) Giulio Cosentino(1) Paolo Ricci  Clément Henry  Maxime Y.
Monin  Adrien Merlini  Francesco P. Andriulli

###### Abstract

This work focuses on a quasi-linear-in-complexity strategy for a hybrid
surface-wire integral equation solver for the electroencephalography forward
problem. The scheme exploits a block diagonally dominant structure of the wire
self block—that models the neuronal fibers self interactions—and of the
surface self block—modeling interface potentials. This structure leads to two
Neumann iteration schemes further accelerated with adaptive integral methods.
The resulting algorithm is linear up to logarithmic factors. Numerical results
confirm the performance of the method in biomedically relevant scenarios.

On a Fast Solution Strategy for a Surface-Wire Integral Formulation of the
Anisotropic Forward Problem in Electroencephalography

${}^{(\ref{})}$

(1) Galileo Ferraris High School, Turin, Italy

(2) Politecnico di Torino, Turin, Italy

(3) IMT Atlantique, Brest, France

##  I Introduction

Several neuro-pathologies require precise functional brain imaging as part of
their diagnostic or therapeutic protocols (see [1] and references therein).
Among non-invasive strategies, high resolution electroencephalography (HR-
EEG), that images the electric activity of the brain from scalp potentials, is
widely used. In HR-EEG the volume currents are retrieved from the measurements
of the electric potentials on the scalp by solving the EEG inverse problem.
Solving this inverse problem requires multiple solutions of the EEG forward
problem (FP) in which the surface potential generated by a known current
configuration is computed. Boundary element methods (BEMs) are very popular in
the biomedical community to model the FP and a recent hybrid formulation [1]
has introduced the possibility of modeling white matter anisotropies by
coupling surface BEM with an integral equation for partially conducting wires.
In this work we present a fast matrix-vector multiplication algorithm for this
hybrid formulation which, by exploiting the block diagonal dominance structure
(induced by the presence of neuronal fibers in the model) and coupling this
matrix structure with adaptive integral methods, obtains a scheme with
$O(N\log N)$ complexity in the $N$ degrees of freedom. Theoretical and
algorithmic considerations will be complemented by numerical experiments
showing the impact of the formulation on medical scenarios.

##  II Background and Notation

Consider a sequence of nested compartments $\Omega_{i}$, $i=1,\ldots,C$
modeling the different layers of the head medium, characterized by homogeneous
and isotropic conductivities $\sigma_{i}$. The boundary of each compartment is
denoted by $\Gamma_{i}$. Following the strategy in [1], the inhomogeneity and
anisotropy of the head medium is modeled by populating the white matter with
wires of finite anisotropic conductivity contrast
$\bar{\bar{\chi}}(\boldsymbol{r})=\left(\sigma_{i_{w}}\boldsymbol{\mathsf{I}}-\bar{\bar{\sigma}}(\boldsymbol{r})\right)\bar{\bar{\sigma}}^{-1}(\boldsymbol{r})$
with respect to the background conductivity $\sigma_{i_{w}}$ of the white
matter’s compartment. In this setting, the EEG FP consists in finding the
electric potential $\phi(\boldsymbol{r})$ on the scalp surface $\Gamma_{C}$
generated by a primary current $\boldsymbol{J}\\!_{p}(\boldsymbol{r})$. To do
so, the surface $\xi$ and wire $\boldsymbol{J}\\!_{\mathrm{eq}}$ unknowns (see
[1] for their physical definition) are expanded with discrete basis functions,
i.e. $\xi\approx\sum_{i=1}^{N_{s}}(\boldsymbol{\alpha}_{s})_{i}p_{i}$ and
$\boldsymbol{J}\\!_{\mathrm{eq}}\approx\sum_{i=1}^{N_{w}}(\boldsymbol{\alpha}_{w})_{i}\boldsymbol{h}_{i}$
where $p_{i}$ and $\boldsymbol{h}_{i}$ are the 2D and 1D linear Lagrange
interpolants, respectively. Following a Galerkin approach leads to a linear
system of $N=N_{s}+N_{w}$ unknowns

$\boldsymbol{\mathsf{Z}}\boldsymbol{\alpha}=\begin{bmatrix}\boldsymbol{v}_{s}\\\
\boldsymbol{v}_{w}\end{bmatrix}\text{ with
}\boldsymbol{\mathsf{Z}}=\begin{bmatrix}-\boldsymbol{\mathsf{G}}_{ss}+\boldsymbol{\mathsf{D}}^{*}_{ss}&-\boldsymbol{\mathsf{D}}^{*}_{sw}\\\
-\boldsymbol{\mathsf{S}}_{ws}&\boldsymbol{\mathsf{G}}_{ww}+\boldsymbol{\mathsf{S}}_{ww}\end{bmatrix}\,,$
(1)

and where
$\boldsymbol{\alpha}=[\boldsymbol{\alpha}_{s}\,\boldsymbol{\alpha}_{w}]^{\mathsf{T}}$,
$(\boldsymbol{\mathsf{G}}_{ss})_{ij}=\tfrac{\sigma_{i_{n}}+\sigma_{i_{n}+1}}{2(\sigma_{i_{n}+1}-\sigma_{i_{n}})}\langle
p_{i},p_{j}\rangle_{\Gamma}$, $(\boldsymbol{\mathsf{D}}^{*}_{ss})_{ij}=\langle
p_{i},\mathcal{D}^{*}p_{j}\rangle_{\Gamma}$,
$(\boldsymbol{\mathsf{D}}^{*}_{sw})_{ij}=\tfrac{1}{\sigma_{i_{n}}}\langle
p_{i},\mathcal{D}^{*}_{v}\bar{\bar{\chi}}\boldsymbol{h}_{j}\rangle_{\Gamma}$,
$(\boldsymbol{\mathsf{S}}_{ws})_{ij}=\langle\boldsymbol{h}_{i},\nabla\mathcal{S}p_{j}\rangle_{\Omega}$,
$(\boldsymbol{\mathsf{G}}_{ww})_{ij}=\langle\boldsymbol{h}_{i},\bar{\bar{\chi}}^{-1}\boldsymbol{h}_{j}\rangle_{\Omega}$,
$(\boldsymbol{\mathsf{S}}_{ww})_{ij}=\tfrac{1}{\sigma_{i_{n}}}\langle\boldsymbol{h}_{i},\nabla\mathcal{S}_{v}\boldsymbol{h}_{j}\rangle_{\Omega}$,
$(\boldsymbol{v}_{s})_{i}=-\tfrac{1}{\sigma_{s}}\langle
p_{i},\mathcal{D}^{*}_{v}\boldsymbol{J}\\!_{p}\rangle_{\Gamma}$, and
$(\boldsymbol{v}_{w})_{i}=-\tfrac{1}{\sigma_{p}}\langle\boldsymbol{h}_{i},\nabla\mathcal{S}_{v}\boldsymbol{J}\\!_{p}\rangle_{\Omega}$,
with
$(\mathcal{S}f)(\boldsymbol{r})=\int_{S}G(\boldsymbol{r},\boldsymbol{r}^{\prime})f(\boldsymbol{r}^{\prime})\,\mathrm{d}S^{\prime}$,
$(\mathcal{D}^{*}f)(\boldsymbol{r})=\int_{S}\hat{\boldsymbol{n}}\cdot\nabla
G(\boldsymbol{r},\boldsymbol{r}^{\prime})f(\boldsymbol{r}^{\prime})\,\mathrm{d}S^{\prime}$,
$(\mathcal{S}_{v}\boldsymbol{f})(\boldsymbol{r})=\int_{V}G(\boldsymbol{r},\boldsymbol{r}^{\prime})\nabla^{\prime}\cdot\boldsymbol{f}(\boldsymbol{r}^{\prime})\,\mathrm{d}V^{\prime}$,
and
$(\mathcal{D}^{*}_{v}\boldsymbol{f})(\boldsymbol{r})=\int_{V}\hat{\boldsymbol{n}}\cdot\nabla
G(\boldsymbol{r},\boldsymbol{r}^{\prime})\nabla^{\prime}\cdot\boldsymbol{f}(\boldsymbol{r}^{\prime})\,\mathrm{d}V^{\prime}$.
Above $\hat{\boldsymbol{n}}$ denotes the unit normal vector pointing outwards
$\Gamma_{i}$ and
$G(\boldsymbol{r},\boldsymbol{r}^{\prime})=\frac{1}{4\pi\left\|\boldsymbol{r}-\boldsymbol{r}^{\prime}\right\|}$
is the static Green function. Once (1) is solved, $\mathcal{S}$ and
$\mathcal{S}_{v}$ can be applied to $\boldsymbol{\alpha}$ to get the potential
$\phi(\boldsymbol{r})$ on $\Gamma_{C}$.

##  III A Fast Solution Strategy

With respect to a standard integral formulation for isotropic media,
corresponding to the left diagonal block in (1), the inclusion of the white
matter anisotropy adds a new wire-wire diagonal block and two coupling blocks
in the system. First, the new scheme aims at decoupling the surface and wire
solution via block diagonal inversion and Neumann series solution of the
remainder: after separating diagonal and off-diagonal blocks:
$\boldsymbol{\mathsf{Z}}=\boldsymbol{\mathsf{Z}}_{\mathrm{self}}+\boldsymbol{\mathsf{Z}}_{\mathrm{coupl}}$
with
$\boldsymbol{\mathsf{Z}}_{\mathrm{self}}=\left[\boldsymbol{\mathsf{Z}}_{ss},\boldsymbol{\mathsf{0}};\boldsymbol{\mathsf{0}},\boldsymbol{\mathsf{Z}}_{ww}\right]$
and
$\boldsymbol{\mathsf{Z}}_{\mathrm{coupl}}=\left[\boldsymbol{\mathsf{0}},\boldsymbol{\mathsf{Z}}_{sw};\boldsymbol{\mathsf{Z}}_{ws},\boldsymbol{\mathsf{0}}\right]$,
we solve (1) as
$\big{(}\boldsymbol{\mathsf{I}}+\boldsymbol{\mathsf{Z}}_{\mathrm{self}}^{-1}\boldsymbol{\mathsf{Z}}_{\mathrm{coupl}}\big{)}\boldsymbol{\alpha}=\boldsymbol{\mathsf{Z}}_{\mathrm{self}}^{-1}\boldsymbol{v}$
via a Neumann series approach enabled by the block diagonal dominance, in
cases of practical relevance, of the original matrix (i.e. for the spectral
radius
$\rho_{Z}=\rho\\!\left(\boldsymbol{\mathsf{Z}}_{\mathrm{self}}^{-1}\boldsymbol{\mathsf{Z}}_{\mathrm{coupl}}\right)<1$).
Thus we have
$\boldsymbol{\alpha}_{s}^{(k+1)}=\boldsymbol{\mathsf{Z}}_{ss}^{-1}(\boldsymbol{v}_{s}-\boldsymbol{\mathsf{Z}}_{sw}\boldsymbol{\alpha}_{w}^{(k)})$
and
$\boldsymbol{\alpha}_{w}^{(k+1)}=\boldsymbol{\mathsf{Z}}_{ww}^{-1}(\boldsymbol{v}_{w}-\boldsymbol{\mathsf{Z}}_{ws}\boldsymbol{\alpha}_{s}^{(k)})$
whose complexity reduces to the one of the two inversions and of the
multiplication of the coupling terms. The multiplication of the coupling terms
can be done efficiently if a fast matrix vector product algorithm is
available. We have opted for an adaptive integral method (AIM) [2]. In other
words, all kernel interactions in $\boldsymbol{\mathsf{D}}^{*}_{ss}$,
$\boldsymbol{\mathsf{S}}_{ww}$, $\boldsymbol{\mathsf{S}}_{ws}$ and
$\boldsymbol{\mathsf{D}}^{*}_{sw}$ between all Gaussian quadrature points are
interpolated on the same Cartesian grid with a number of nodes proportional to
the number of unknowns $N$ and handled via FFT in $\mathcal{O}\\!\left(N\log
N\right)$ complexity. As is standard in AIM [2], a near field precorrection is
required for all kernels: a generic $\boldsymbol{\mathsf{D}}^{*}$ and
$\boldsymbol{\mathsf{S}}$ (for surface, wire, and off diagonal couplings) is
written as
$\boldsymbol{\mathsf{D}}^{*}=\boldsymbol{\mathsf{D}}^{*}_{\mathrm{near}}-\widetilde{\boldsymbol{\mathsf{D}}}^{*}_{\mathrm{near}}+\boldsymbol{\mathsf{\Phi}}_{p}\boldsymbol{\mathsf{\Lambda}}^{\mathsf{T}}\boldsymbol{\mathsf{g}}_{D}\boldsymbol{\mathsf{\Lambda}}\boldsymbol{\mathsf{\Phi}}_{f}$
and
$\boldsymbol{\mathsf{S}}=\boldsymbol{\mathsf{S}}_{\mathrm{near}}-\widetilde{\boldsymbol{\mathsf{S}}}_{\mathrm{near}}+\boldsymbol{\mathsf{\Phi}}_{p}\boldsymbol{\mathsf{\Lambda}}^{\mathsf{T}}\boldsymbol{\mathsf{g}}_{S}\boldsymbol{\mathsf{\Lambda}}\boldsymbol{\mathsf{\Phi}}_{f}$,
where $\boldsymbol{\mathsf{D}}^{*}_{\mathrm{near}}$ and
$\boldsymbol{\mathsf{S}}_{\mathrm{near}}$ are the uncompressed near fields,
$\widetilde{\boldsymbol{\mathsf{D}}}^{*}_{\mathrm{near}}$ and
$\widetilde{\boldsymbol{\mathsf{S}}}_{\mathrm{near}}$ are the FFT
precorrections, $\boldsymbol{\mathsf{\Lambda}}$ is unique for every product
and interpolates the quadrature points, and $\boldsymbol{\mathsf{\Phi}}_{p}$
and $\boldsymbol{\mathsf{\Phi}}_{f}$ map quadrature points to basis functions;
all these matrices are sparse. The FFT is applied to the Toeplitz matrices
$\boldsymbol{\mathsf{g}}_{S}$ and $\boldsymbol{\mathsf{g}}_{D}$ that, because
of the translation invariance of all Green functions involved, require
$\mathcal{O}\\!\left(N\right)$ memory storage. Since the double layer kernel
is $\hat{\boldsymbol{n}}\cdot\nabla
G(\boldsymbol{r},\boldsymbol{r}^{\prime})=\hat{\boldsymbol{n}}\cdot\frac{\boldsymbol{r}^{\prime}-\boldsymbol{r}}{4\pi\left\|\boldsymbol{r}^{\prime}-\boldsymbol{r}\right\|^{3}}$,
the product of $\boldsymbol{\mathsf{g}}_{D}$ with a vector is split into three
scalar components.

Since the $\boldsymbol{\mathsf{Z}}_{ss}$ block corresponds to the classical
homogeneous multilayer BEM formulation, once a fast matrix vector product
algorithm is available, it can be inverted iteratively with standard
techniques (see [1] and references therein). Regarding
$\boldsymbol{\mathsf{Z}}_{ww}$, the near field kernel interactions are
extracted with an octree and the resulting sparse matrix
$\boldsymbol{\mathsf{N}}$ is used as a preconditioner of the linear system
$\boldsymbol{\mathsf{Z}}_{ww}\boldsymbol{x}=\boldsymbol{b}$. The near field
dominance of $\boldsymbol{\mathsf{Z}}_{ww}$—due to the electric current
flowing along the fibers, i.e.
$\rho_{w}=\rho\\!\left(\boldsymbol{\mathsf{N}}^{-1}(\boldsymbol{\mathsf{Z}}_{ww}-\boldsymbol{\mathsf{N}})\right)<1$—enables
a second usage of a Neumann series from which
$\boldsymbol{x}^{(k+1)}=\boldsymbol{x}^{(k)}+\boldsymbol{\mathsf{N}}^{-1}\big{(}\boldsymbol{b}-\boldsymbol{\mathsf{Z}}_{ww}\boldsymbol{x}^{(k)}\big{)}$.
A sparse solver is used to invert $\boldsymbol{\mathsf{N}}$ in
$\mathcal{O}\\!\left(N\right)$ time complexity and the multiplication of
$\boldsymbol{\mathsf{Z}}_{ww}$ is done in $\mathcal{O}\\!\left(N\log N\right)$
with the AIM.

TABLE I: Performance Comparison | Setup time | Storage | Time per RHS
---|---|---|---
Standard Iterative Solution | $24\,657\text{\,}\mathrm{s}$ | $31.2\text{\,}\mathrm{G}\mathrm{B}$ | $89.15\text{\,}\mathrm{s}$
This work | $2498\text{\,}\mathrm{s}$ | $0.4\text{\,}\mathrm{G}\mathrm{B}$ | $26.42\text{\,}\mathrm{s}$

Figure 1: Solution time vs number of unknowns.

##  IV Numerical Results

The favorable complexity scaling of the proposed scheme has been verified on a
set of canonical geometries composed of spherical surfaces and orthogonal
brain fibers. The total timings are reported in Fig. 1 and clearly confirm
that the scheme we propose is, up to logarithmic factors, linear in
complexity. The relevance of our fast solution strategy for real case
scenarios has been studied on a realistic head model obtained with magnetic
resonance imaging (MRI) data that includes white matter neuronal fibers with a
tangential anisotropic conductivity of
$1.3\text{\,}\mathrm{S}\text{\,}{\mathrm{m}}^{-1}$ and four layers (gray
matter, cerebrospinal fluid, skull, scalp) with conductivities
$0.13\text{\,}\mathrm{S}\text{\,}{\mathrm{m}}^{-1}$,
$1.79\text{\,}\mathrm{S}\text{\,}{\mathrm{m}}^{-1}$,
$0.01\text{\,}\mathrm{S}\text{\,}{\mathrm{m}}^{-1}$, and
$0.43\text{\,}\mathrm{S}\text{\,}{\mathrm{m}}^{-1}$ respectively. The obtained
current on the neuronal fibers is shown in Fig. 2. For this problem the radius
of the fibers is chosen to match a total volume of
$450\text{\,}{\mathrm{mm}}^{3}$. The total number of unknowns is $63\,922$.
The two spectral radii are $\rho_{Z}=0.439$ and $\rho_{w}=0.799$, both less
than one, thus allowing the Neumann strategy. For this experiment we have
compared in Table I the method proposed in this work with the uncompressed
solution. In both cases the tolerance iterative schemes has been set to
$10^{-3}$ and the results show the advantage of the new scheme.

Figure 2: Brain fibers current density on the MRI head model.

## Acknowledgment and Contributions

The work of this paper and the associated Early-Research Program for talented
high school students has received funding from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement No 724846, project 321) and from the EU H2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement n° 955476 (project COMPETE). The first two authors are listed in
alphabetic order and have contributed equally to this work.

## References

  * [1] M. Y. Monin, L. Rahmouni, A. Merlini and F. P. Andriulli, “A hybrid volume-surface-wire integral equation for the anisotropic forward problem in electroencephalography,” IEEE J. Electromagn. RF Microw. Med. Biol., vol. 4, no. 4, pp. 286–293, Dec. 2020.
  * [2] E. Bleszynski, M. Bleszynski and T. Jaroszewicz, “AIM: Adaptive integral method for solving large-scale electromagnetic scattering and radiation problems,” in Radio Science, vol. 31, no. 5, pp. 1225–1251, Sept.-Oct. 1996.