Solving for the low-rank tensor components of a scattering wave function
Jacob Snoeijer and Wim Vanroose

Atomic and molecular breakup reactions, such as multiple-ionisation,
are described by a driven Schrödinger equation.  This equation is
equivalent to a high-dimensional Helmholtz equation and it has
solutions that are outgoing waves, emerging from the target.  We
show that these waves can be described by a low-rank
approximation.  For 2D problems this it a matrix product of two
low-rank matrices, for 3D problems it is a low-rank tensor decomposition.
We propose an iterative method that solves, in an alternating way, for
these low-rank components of the scattered wave. We illustrate the
method with examples in 2D and 3D.

§ INTRODUCTION

An in-coincidence experiment measures simultaneously the outgoing
momenta of multiple products of a microscopic reaction
[21]. It is an instrument that can study the
correlations in reactions involving multiple particles.  In double
ionization, for example, a single photon ionizes, simultaneously, two
electrons and the outgoing momenta of both particles are captured
[1]. The reaction probes the correlation between
two electrons in, for example, a chemical bound at the moment of
photon impact.  The outgoing wave of the two electrons is described by
a 6D correlated wave and results in a cross section that depends on
four angles, the directions of the first and the second electron.

Free-electron lasers, and similar experiments around the world, are
expected to generate a wealth of this high-dimensional scattering
data.  This will result in high-dimensional forward and inverse wave
problems that need to be solved to interpret the data.

The experimental cross section are often smooth functions as a
function of the angles. Similarly, some parts of the scattering
solution, such as single-ionization, only probes a limited subspace
of the possible full solution space. The scattering solution can then
be described by a low-rank wave function, a product of one-particle
bound states with scattering waves in the other coordinates.

This paper introduces a low-rank representation for the
scattering solutions, not only for the single ionization but also for
double and triple ionization waves that appear in breakup reaction.

We also propose and analyze an alternating direction algorithm that
directly solves for the low-rank components that describe the
solution. This reduces a large-scale linear system to smaller,
low-dimensional, scattering problems that are solved in a iterative
sequence.  The proposed method can be generalized to high-dimensional
scattering problems where a low-rank tensor decomposition is used to
represent the full scattering wave function.

Efficient low-rank tensor representations are used in quantum physics
for quite some time already [11, 14].  They are also used in the applied mathematics
literature to approximate high-dimension problems, for a review see
[8, 9, 12].
Methods such ALS [10], DMRG
[16], and AMEn [6] use in
alternating directions, a small linear system to determine the
low-rank components of a tensor decomposition. These innovations have
not found their application in computational scattering theory.

To calculate cross sections, from first principles, we start
from a multi-particle Schrödinger equation. The equation is
reformulated into a driven Schrödinger equation with an unknown scattering
wave function and a right hand side that describes the
excitation, for example, a dipole operator working on the initial

Since the asymptotic behaviour of a scattering function for multiple
charged particles is in many cases unknown, absorbing boundary
conditions [22, 2] are used.  Here, an artificial layer is
added to the numerical domain that dampens outgoing waves. The
outgoing wave boundary conditions are then replaced with
homogeneous Dirichlet boundary conditions at the end of the artificial
layer. These boundary do not require any knowledge about the
asymptotic behaviour, which becomes very complicated for these
multiple charged particles.

The resulting equation is discretized on a grid and results in a
large, sparse indefinite linear system.  It is typically solved by a
preconditioned Krylov subspace method [5].  However,
the preconditioning techniques for indefinite systems are not as
efficient as preconditioners for symmetric and positive definite
system.  And solving the resulting equation is still a
computationally expensive task, often requiring a distributed
calculation on a supercomputer.

To compare the resulting theoretical cross sections with experimental
data, a further postprocessing step is necessary. The cross section is
the farfield map and this is calculated through integrals of the
scattering wave function, which is the solution of the linear system,
and a Greens function [13].

The main result of the paper is that we show that scattering waves
that describe multiple ionization can be represented by a low-rank
tensor. We first show this for a 2D wave and then generalize the
results to 3D waves.  The methodology can be generalized to higher

The outline of the paper is a follows. In section <ref> we
review the methodology that solves a forward scattering problem. It
results in a driven Schrödinger equation with absorbing boundary
conditions.  From the solution we can extract the cross section using
an integral.  In section <ref> we illustrate, in 2D, that a
solution can be approximated by a truncated low-rank approximation. We
also show that these low-rank components can be calculated directly
with an iterative method.  In section <ref> we show that
this methodology generalizes to 3D and higher dimensional problems.
We use a truncated tensor decomposition and determine the components
with a similar iterative method.
A discussion of some numerical results and a comparison of the different
presented versions of the method is given in section <ref>.
In the final section, Sec. <ref>, we summarize some
conclusions and discuss some possible extensions of the presented method.

§ STATE OF THE ART

This section summarize the methodology that solves a forward
break-up problems with charged particles. The methodology is developed
in a series of papers [20, 13]
and applied to solve the impact-ionization problem
[19] and double ionization of molecules
[27, 26].  These methods are
being extended to treat, for example, water

The helium atom, He, is the simplest system with double ionization
[3]. It has two electrons with coordinates
$\mathbf{r}_1 \in \mathbb{R}^3$ and $\mathbf{r}_2 \in \mathbb{R}^3$
relative to the nucleus positioned at $\mathbf{0}$.  The driven
Schrödinger equation for $u(\mathbf{r}_1,\mathbf{r}_2) \in C^2$ then

\begin{equation}\label{DrivenSchrodinger}
\left( {- \frac{1}{2} \Delta_{\mathbf{r}_1} - \frac{1}{2} \Delta_{\mathbf{r}_2} }  { - \frac{1}{\|\mathbf{r}_1\|}- \frac{1}{\|\mathbf{r}_2\|}} + {\frac{1}{\|\mathbf{r}_1-\mathbf{r}_2\|} }- { E} \right) u(\mathbf{r}_1,\mathbf{r}_2) =   \mathcal{\mu}  \phi_0(\mathbf{r}_1,\mathbf{r}_2) \quad  \forall \mathbf{r}_1, \mathbf{r}_2 \in \mathbb{R}^3,
\end{equation}

where the right hand side is
the dipole operator $\mathcal{\mu}$ working on the ground state $\phi_0$, the
eigenstate with the lowest energy $\lambda_0$.  The operators
$-\frac{1}{2} \Delta_{\mathbf{r}_1}$ and
$-\frac{1}{2}\Delta_{\mathbf{r}_2}$ are the Laplacian operators for
the first and second electron and model the kinetic energy.  The
nuclear attraction is $-1/\|\mathbf{r}_1\|$ and $-1/\|\mathbf{r}_2\|$
and the electron-electron repulsion is

The total energy $E = h\nu + \lambda_0$ is the energy
deposited in the system by the photon, $h\nu$, and the energy $\lambda_0$ of
the ground state. If the $E>0$, both electrons can
escape simultaneously from the system. The solution
$u(\mathbf{r}_1,\mathbf{r}_2)$ then represents a 6D wave emerging from
the nucleus.

The equation can be interpreted as a Helmholtz equation with a
space-dependent wave number, $k^2(\mathbf{r}_1,\mathbf{r}_2)$,

\begin{equation}\label{eq:helmholtz}
\left( - \Delta_{6D} - {  k^2(\mathbf{r}_1,\mathbf{r}_2)} \right)
u(\mathbf{r}_1,\mathbf{r}_2) = {  f(\mathbf{r}_1,\mathbf{r}_2)} \quad \forall (\mathbf{r}_1,
\mathbf{r}_2) \in \mathbb{R}^6.
\end{equation}

In this paper we prefer to write this Helmholtz equation as

\begin{equation}
\left( - \Delta_{6D}   -   k_0^2  (1 + \chi(\mathbf{r}_1,\mathbf{r}_2) ) \right)  u(\mathbf{r}_1,\mathbf{r}_2)  =  f(\mathbf{r}_1,\mathbf{r}_2) \quad  \forall (\mathbf{r}_1, \mathbf{r}_2) \in \mathbb{R}^6,
\end{equation}

where $k_0^2$ is a constant wave number, in this case related to the
total energy $E$, and a space-dependent function $\chi: \mathbb{R}^6
\rightarrow \mathbb{R}$, that goes to zero if $\|\mathbf{r}_1\|
\rightarrow \infty $ or $\|\mathbf{r}_1\| \rightarrow \infty $ that
represents all the potentials.

§.§ Expansion in spherical waves and absorbing boundary conditions

For small atomic and molecular system, where spherical symmetry is relevant, the
system is typically written in spherical coordinates and expanded in
spherical harmonics.  With $\mathbf{r}_1(\rho_1, \theta_1,
\varphi_1)$ and $\mathbf{r}_2(\rho_2, \theta_2, \varphi_2)$ we can

\begin{equation}\label{eq:partialwaves}
u(\mathbf{r}_1,\mathbf{r}_2) = \sum_{l_1}^\infty \sum_{m_1 = -l_1}^{l_1} \sum_{l_2=0}^{\infty} \sum_{m_2=-l_2}^{l_2}  u_{l_1m_1,l_2m_2}(\rho_1,\rho_2)  Y_{l_1m_1}(\theta_1,\varphi_1) Y_{l_2m_2}(\theta_2,\varphi_2),
\end{equation}

where $Y_{l_1m_1}(\theta_1, \varphi_1)$ and $Y_{l_2m_2}(\theta_2,
\varphi_2)$ are spherical Harmonics, the eigenfunctions of the angular
part of a 3D Laplacian in spherical coordinates. In practice the sum in equation (<ref>) is
truncated. The expansion is then a low-rank, truncated, tensor decomposition
of a 6D tensor describing the solution.

For each $l_1$, $m_1$, $l_2$ and $m_2$ combination, the radial function
$u_{l_1m_1,l_2m_2}(\rho_1,\rho_2)$ describes an outgoing wave that
depends on the distances $\rho_1$ and $\rho_2$ of the two electrons to
the nucleus.

A coupled equation that simultaneously solves for all the
$u_{l_1m_1l_2m_2}(\rho_1,\rho_2)$'s is found by inserting the
truncated sum in (<ref>), multiplying with $Y_{l_1m_1}^*(\theta_1,
\varphi_1)$ and $Y_{l_2m_2}^*(\theta_2,\varphi_2)$ and integrating
over all the angular coordinates,

\begin{equation}\label{eq:coupledpartialwaves}
\begin{aligned}\left(-\frac{1}{2}\frac{d^2}{d \rho_1^2} + \frac{l_1(l_1+1)}{2\rho_1^2 }
-\frac{1}{2}\frac{d^2}{d \rho_2^2} + \frac{l_2(l_2+1)}{2\rho_2^2 } + V_{l_1m_1l_2m_2}
(\rho_1,\rho_2) -E \right) u_{l_1m_1l_2m_2}(\rho_1,\rho_2)   \\ + \sum_{l^\prime_1m^\prime_1l^\prime_2m^\prime_2} V_{l_1m_1l_2m_2, l^\prime_1m^\prime_1l^\prime_2m^\prime_2}(\rho_1,\rho_2) u_{l^\prime_1m^\prime_1l^\prime_2m^\prime_2}(\rho_1,\rho_2) = f_{l_1m_1,l_2m_2}(\rho_1,\rho_2) \quad \forall l_1m_1 l_1m_2. \quad \forall \rho_1,\rho_2 \in [0,\infty[
\end{aligned}
\end{equation}

with boundary conditions $u(\rho_1=0,\rho_2)=0$ for all $\rho_2 \ge 0$ and $u(\rho_1,\rho_2=0)=0$ for all $\rho_1 \ge 0$.

The equation (<ref>) is typically discretized
on a spectral elements quadrature grid [20].

To reflect the physics, where electrons are emitted from the system,
outgoing wave boundary conditions need to be applied at the outer
boundaries. There are many ways to implement outgoing wave boundary
conditions. Exterior complex scaling (ECS) [22] for
example, is frequently used in the computational atomic and molecular
physics literature. In the computational electromagnetic scattering a
perfectly matched layers (PML) [2] is used, which can also be
interpreted as a complex scaled grid. [4]

§.§ Calculation of the amplitudes

To correctly predict the probabilities of the arriving particles at
the detector, we need the amplitudes of the solution far away from the
molecule. These are related to the asymptotic amplitudes of the
wave functions.

Let us go back to the formulation with Helmholtz equation, as given in
(<ref>). Suppose that we have solved the following
Helmholtz equation with absorbing boundary conditions, in any

\begin{equation} \label{Helmholtz}
\left(- \Delta -k_0^2 \left( 1+ \chi(\mathbf{x}) \right) \right) u_{\text{sc}}(\mathbf{x}) = f(\mathbf{x}),  \quad \forall \mathbf{x} \in [-L,L]^d,
\end{equation}

where $f$ is only non-zero on the real part of the grid $[-L,L]^d
\subset \mathbb{R}^d$. Similarly, $\chi(\mathbf{x})$ is only non-zero
on the box $[-L,L]^d$.

The calculation of the asymptotic amplitudes requires the solution
$u_{\text{sc}}(\mathbf{x})$ for an $\mathbf{x}$ outside of the box
$[-L, L]^d$. To that end, we reorganize equation (<ref>),
after we have solved it, as follows

\begin{equation} \label{eq:reorganized}
\left( - \Delta -k_0^2 \right) u_{\text{sc}} = f + k_0^2  \, \chi \, u_{\text{sc}}.
\end{equation}

The right hand side of (<ref>) is now only non-zero on
$[-L,L]^d$, since both $f$ and $\chi$ are only non-zero
there. Furthermore, since we have solved (<ref>) we also
know $u_\text{sc}$ on $[-L,L]^d$.  So the full right hand side of
(<ref>) is known. The remaining left hand side of
(<ref>) is now a Helmholtz equation with a constant
wave number $k_0^2$.  For this equation the Greens function is known

\begin{equation}
\begin{aligned}
u_{\text{sc}}(\mathbf{x}) &= \int  G(\mathbf{x},\mathbf{y}) \left( f(\mathbf{y}) +   k_0^2 \,\chi(\mathbf{y}) u_\text{sc}(\mathbf{y})\right) d\mathbf{y} \\
&=   \int\limits_{[-L,L]^d} G(\mathbf{x},\mathbf{y}) \left( f(\mathbf{y}) +   k_0^2\, \chi(\mathbf{y}) u_\text{sc} (\mathbf{y}) \right) d\mathbf{y} \quad \forall x \in \mathbb{R}^d,
\end{aligned}
\end{equation}

where $f$ and $\chi$ are limited to $[-L,L]^d$ thus we can truncate the
integral to the box $[-L,L]^d$.

This methodology was successfully applied to calculate challenging
break up problems, for example [19].

§.§ Single ionization versus double ionization

Let us discuss the qualitative behaviour of the solution for single
and double ionization. To illustrate the behaviour, we truncate the
partial wave expansion, (<ref>), to the first
term. This is known as a $s$-wave expansion. The 6D wave function is
then approximated as

\begin{equation}
u(\mathbf{r}_1,\mathbf{r}_2) \approx  u(\rho_1,\rho_2)  Y_{00}(\theta_1,\varphi_1) Y_{00}(\theta_2,\varphi_2).
\end{equation}

The radial wave, $u(\rho_1,\rho_2)$, then fits a 2D Helmholtz equation

\begin{equation}\label{eq:radialequation}
\left(- \frac{1}{2} \frac{d^2}{d\rho_1^2} - \frac{1}{2} \frac{d^2}{d\rho_2^2} + V_1(\rho_1) + V_2(\rho_2) + V_{12}(\rho_1,\rho_2) - E \right) u(\rho_1,\rho_2) = f(\rho_1,\rho_2), \quad \forall \rho_1,\rho_2 \in [0,\infty[,
\end{equation}

where $V_1(\rho_1)$ and $V_2(\rho_2)$ represents the one-particle
potentials and $V_{12}(\rho_1,\rho_2)$ the two-particle repulsion.
This model is known as a $s$-wave or Temkin-Poet model
[24, 18].

Before the photo-ionization, the atom is in a two-particle ground
state. In this $s$-wave model, it is the eigenstate of

\begin{equation}
\left(- \frac{1}{2} \frac{d^2}{d\rho_1^2} - \frac{1}{2} \frac{d^2}{d\rho_2^2} +
V_1(\rho_1) + V_2(\rho_2) + V_{12}(\rho_1,\rho_2)
\right)\phi_0(\rho_1,\rho_2) = \lambda_0 \phi_0(\rho_1,\rho_2).
\end{equation}

with the lowest energy. Simultaniously,
there are one-particle states that are eigenstates of

\begin{equation}\label{eq:eigenstate_H1}
\left(-\frac{1}{2} \frac{d^2}{d\rho_1^2} + V_1(\rho_1) \right) \phi_i(\rho_1)  =  \mu_i \phi_i(\rho_1)
\end{equation}

\begin{equation}
\label{eq:eigenstate_H2}
\left(- \frac{1}{2} \frac{d^2}{d\rho_2^2} + V_2(\rho_2) \right) \varphi_i(\rho_2)  =  \nu_i \varphi_i(\rho_2).
\end{equation}

When (<ref>) is solved with the energy $E = \hbar
\nu + \lambda<0$, there is only single ionization. Only one of the two
coordinates $\rho_1$ or $\rho_2$ can become large and the solution,
as can be seen in Figure <ref>, is
localized along both axis. The solution is a product of an outgoing
wave in one coordinate and a bound state in the other
coordinate. For example, along the $\rho_2$-axis, the solution is
described by $A_i(\rho_2)\phi_i(\rho_1)$, where $A_i(\rho_2)$ is a
one-dimensional outgoing wave, with an energy $E-\mu_i$ and
$\phi_i(\rho_1)$ is a bound state of (<ref>) in the
first coordinate with energy $\mu_i$.  Similary, there is a wave,
along the $\rho_1$ axis, that is an outgoing wave of the form
$B_l(\rho_1) \varphi_l(\rho_2)$, with a scattering wave in the first coordinate, $\rho_1$,
and a bound state in the coordinate $\rho_2$, solution of

When (<ref>) is solved with energy $E = \hbar \nu +
\lambda \ge 0$ there is also double ionization and both coordinates
$\rho_1$ and $\rho_2$ can become large. We see, in
Figure <ref>, a (spherical) wave in the middle of
the domain, where both coordinates can be become large. To describe
this solution the full coordinate space is necessary.  Note that
these solutions still show single ionization along the axes. Even for
$E> 0$, one particle can take away all the energy and leave the other
particle as a bound state.

Single ionization

Double ionization

Left: When the energy $E = \hbar \nu + \lambda_0 <0$, there is
only single ionization. The solution is then localized along the
edges, where the solution is a combination of an outgoing wave in
the $\rho_1$ and a bound state in $\rho_2$, or vice-versa. Right: For
energy $E>0$, there is,  in addition to single ionization with
solution localized along the edges, a double ionization wave where
both coordinates can become large.

Rank-1 approximation

Rank-2 approximation

Rank-3 approximation

Rank-4 approximation

Rank-5 approximation

Full rank wave

Contour plots of the double ionization wave function (bottom, right) and low-rank approximations for increasing rank.

§.§ Coupled channel model for single ionization waves

In this section, we write the single ionization solution as a low-rank
decomposition and derive the equations for the low-rank
components. When there is only single ionization, the total wave can
be written as

\begin{equation}\label{expansion}
u(\rho_1,\rho_2) = \sum_{m=1}^{M} \phi_m(\rho_1) A_m(\rho_2)  +  \sum_{l=1}^L  B_l(\rho_1)\varphi_l(\rho_2),
\end{equation}

where $\phi_m(\rho_1)$ and $\varphi_l(\rho_2)$ are the  bound state eigenstates,
defined in (<ref>) and (<ref>).  The
first term is localized along the $\rho_2$-axis, the second term is
localized along the $\rho_1$-axis with $\mu_i < 0$ and $\nu_i<0$.

As discussed in [5], this expansion is not unique. We
can add multiples of $\gamma_i \varphi_i(\rho_2)$ to $A_i(\rho_2)$ and
simultaneously subtract $\gamma_l \phi_l(\rho_1)$ from $B_l(\rho_1)$
without contaminating the result.  Indeed, for any choice of $\gamma_i
\in \mathbb{C}$ and $L=M$ holds that

\begin{equation}
u(\rho_1,\rho_2) = \sum_{m=1}^{M} \phi_m(\rho_1) \left(A_m(\rho_2) + \gamma_m \varphi_m(\rho_2) \right)  +  \sum_{l=1}^L  \left(B_l(\rho_1)-\gamma_i \phi_i(\rho_1))\varphi_l(\rho_2) \right) = u(\rho_1,\rho_2).
\end{equation}

To make the expansion unique, [5] chooses to select $A_i
\perp \varphi_j$ when $j \ge i$ and $B_i\perp \phi_m$ when $i \ge j$.

In this paper, we choose to make the functions in the set $\{\phi_{i
\in\{1,\ldots,m\}}, B_{l \in \{1,\ldots, L \}}\}$ orthogonal.  We
also assume that $V_{12}(\rho_1,\rho_2) \approx \sum_{i=1}^L\sum_{l}^M
\phi_i(\rho_1) \varphi(\rho_2) \int
\phi_i^*(\rho_1)\varphi^*_l(\rho_2) V_{12}(\rho_1,\rho_2) d\rho_1

Given a function $f(\rho_1,\rho_2)$, a right hand side, we can now
derive the equations for $A_m$ and $B_l$. When we insert the low-rank
decomposition of the expansion (<ref>), in the 2D Hamiltonian
(<ref>), multiply with $\phi_i$ and integrate over
$\rho_1$ we find

\begin{equation}
(H_2 + \mu_j -E) A_j(\rho_2) + \sum_{l=1}^M V_{ij}(\rho_2) A_l(\rho_2) =  \int_0^\infty \phi^*_j(\rho_1) f(\rho_1,\rho_2) d\rho_1, \quad \text{for} \quad j=1, \ldots, M \quad \text{and}\quad  \forall \rho_2 \in [0,L],
\end{equation}

\begin{equation}
V_{ij}(\rho_2)=  \int_0^\infty \phi^*_j(\rho_1) V_{12}(\rho_1,\rho_2) \phi_l(\rho_1) d\rho_1 .
\end{equation}

We have used that $\phi_i \perp B_l$ to eliminate the second term in
the expansion (<ref>).

Similarly, for $B_l$, we find

\begin{equation}
\begin{aligned}
(H_1 + \nu_l -E) B_l(\rho_1) + \sum_{m=1}^L W_{l}(\rho_1) B_k(\rho_1) =  \int_0^\infty\!\!\!\! \varphi_l(\rho_2) \left(f(\rho_1,\rho_2)-\sum_i^M \phi_i \phi^* f(\rho_1,\rho_2) \right) d\rho_2 \quad\\
\text{for} \quad i=1,\ldots, M, \quad \forall \rho_1 \in [0,\infty[
\end{aligned}
\end{equation}

\begin{equation}
W_{lk}(\rho_1):=\left(\int_0^\infty \!\!\!\!  \varphi^*_l(\rho_2) V_{12}(\rho_1,\rho_2) \varphi_k(\rho_2) d\rho_2 \right).
\end{equation}

§ LOW-RANK MATRIX REPRESENTATION OF A 2D WAVE FUNCTION THAT INCLUDES BOTH SINGLE AND DOUBLE IONIZATION

§.§ Low rank of the double ionization solution

We now discuss the main result of the paper. We will derive a
coupled channel equation that gives a low-rank approximation for the
double ionization wave function, as shown in Figure <ref>.

In section <ref> we have shown
that the single ionization wave can be represented by a low-rank
decomposition. In this section, we show that also the double
ionization wave can be written as a similar low-rank decomposition.

We first illustrate that the solution of a 2D driven Schrödinger equation
that contains both single and double ionization, it is a solution of
(<ref>) with $E>0$, can be represented by a similar low-rank

In Figure <ref> we solve the Helmholtz equation with a
space-dependent wave number, $k(x,y)$, in the first quadrant where $x
\ge 0$ and $y \ge 0$. The equation is

\begin{equation}\label{driven}
\left(  - \Delta_2 - k^2(x,y) \right) u_{sc}(x,y) = f(x,y),
\end{equation}

where $\Delta_2$ is the 2D Laplacian and the solution $u_{sc}$
satisfies homogeneous boundary conditions $u_{sc}(x,0)=0$ for all
$x\ge 0$ and $u_{sc}(0,y)=0$ for all $y \ge 0$. On the other
boundaries we have outgoing boundary conditions.

The right hand side $f(x,y)$ has a support that is limited to $[0,b]^2
\subset [0,L]^2 \subset \mathbb{R}_+^2$, i.e. $f(x,y)=0$, for all
$x\ge b$ or $y\ge b$.

The wave number $k(x,y)$ can be split in a constant part, $k_0^2$ and
variable part $\chi(x,y)$. The variable part is also only non-zero on

\begin{equation}
k^2(x,y) = \begin{cases}
k_0^2 \left(1  + \chi(x,y) \right)  \quad &\text{if} \quad  x<b \quad \text{and} \quad y<b,\\
k_0^2   \quad &\text{if} \quad  x\ge b \quad \text{or} \quad y\ge b.\\
\end{cases}
\end{equation}

We extend the domain with exterior complex scaling (ECS) absorbing
boundary condition [13].

The wave function $u_{sc}$ is discretized on the two dimensional mesh and can be
represented as a matrix $\mv{A}\in \mathbb{C}^{n \times n}$. We can
compute the singular decomposition of this matrix, $\mv{A} = \mv{U}
\mv{\Sigma} \mv{V}^\H$ where $U,Σ, V ∈ℂ^n ×n$ where $U^U = I$, $V^V = I$ and $Σ$ is a diagonal matrix with the
singular value $σ_i$ on the diagonal.

\begin{figure}
\begin{tabular}{ll}
\includegraphics[width=0.47\textwidth]{wave-plot-M=1000-E=2.png} &     \includegraphics[width=0.47\textwidth]{wave-plot-M=1000-E=16.png}\\
\input{wave-svd-M=1000-E=2.tikz} & \input{wave-svd-M=1000-E=16.tikz}
\end{tabular}
\caption{The top figures show the wave function for two different
energies (i.e. $E=2$ and $E=16$). The singular values of the
matrix-representation of the discretized functions are shown in
the bottom figures.  we show the solution of a 2D Helmholtz
equation with a space-dependent wave number $k^2(x,y)$ given
by $k^2(x,y) = E - e^{-|x-y|}$. The right hand side $f(x,y)$ on
the finite domain $[0, b[^2$ is given by $f(x,y) = -
e^{-x^2-y^2}$.  Finite difference discretization is done on a
uniform mesh with $M=1000$ interior mesh points per
dimension. At the boundaries $x=b$ and $y=b$ the domain is
extended with exterior complex scaling under an angle
$\frac{\pi}{6}$ where $33\%$ additional discretization points
are added, so $n = 1333$. We have $b=10$.}
\label{fig:wavefunction}
\end{figure}

The results are shown in Figure~\ref{fig:wavefunction} and show that
the singular values rapidly decrease. Thus the wave function can
efficiently be approximated by a truncated representation,
\begin{equation}\label{eq:truncated}
\mv{A} \approx \sum_{i=1}^r \mv{u}_i \sigma_i \mv{v}^\H_i,
\end{equation}
where $u_i ∈ℂ^n$ are columns of $U$ and $v_i ∈ℂ^n$ are rows from $V$ and $σ_i$ are largest $r$ singular values. Thus $A$ is approximated by it low-rank representation with rank $r$.
This truncated decomposition drops all contributions with $σ_i<τ$ below a threshold $τ$, for example the expected discretisation error.

Figure \ref{fig:crosssection} illustrates that a low-rank approximation to the wave function is sufficient
calculate an accurate approximation of the cross section.

\begin{figure}
\begin{center}
\input{sdcs.tikz}
\end{center}
\caption{The cross section based on different low-rank approximations of the wave function.}
\label{fig:crosssection}
\end{figure}

\subsection{Determining the low-rank components directly}
\label{sec:lowrankcomponents}
In the example of the previous section, we have first calculated a
matrix representation of the solution, $A ∈ℂ^n ×n$,
and then approximated it by low-rank components.  The aim is now to
develop a method that calculates, directly, these components without
first calculating the full solution $A$. This approach avoids expensive

We start from the 2D Helmholtz equation as given in \eqref{driven}.
In matrix form this is given by
\begin{equation} \label{eq:matrixform}
- \mv{D}_{xx} \mv{A} - \mv{A} \mv{D}^\T_{yy} - \mv{K} \circ \mv{A} = \mv{F},
\end{equation}
where $D_xx ∈ℂ^n ×n$ and $D_yy ∈ℂ^n ×n$ are sparse matrices that represent the discretization
of the second derivatives, $K$ is the matrix that represents the
space-dependent wave number, $k^2(x,y)$, on the grid and $A ∈ℂ^n×n$ is the matrix that describes the unknown
partial wave. The right hand side $F ∈ℂ^n ×n$ is given. The Hadamard product, $∘$, multiplies the
matrices point wise, element by element.

We now make the approximation  $A ≈UV^$̋, with low-rank matrices
$\mv{U} \in \mathbb{C}^{n \times r}$ and $\mv{V} \in \mathbb{C}^{n
\times r}$ where $r \ll n$ and write

\begin{equation} \label{eq:matrixequation}
- \mv{D}_{xx} \mv{U} \mv{V}^\H - \mv{U} \mv{V}^\H \mv{D}_{yy} - \mv{K} \circ \mv{U} \mv{V}^\H = \mv{F}.
\end{equation}

We start with a guess for $\mv{V} \in \mathbb{C}^{n \times r}$ with
orthogonal columns such that $\mv{V}^\H \mv{V}=\mv{I}_r$. We can now
multiply, (<ref>), from the right,  by $\mv{V}$ and obtain

\begin{equation}\label{eq:2dvarwavenr}
- \mv{D}_{xx} \mv{U} \mv{V}^\H \mv{V} - \mv{U} \mv{V}^\H \mv{D}_{yy} \mv{V}  -  \left(\mv{K} \circ \mv{U}\mv{V}^\H\right) \mv{V} = \mv{F}\mv{V}.
\end{equation}

where $\mv{U} \in \mathbb{C}^{n \times r}$ is now the remaining unknown.  We use
the vectorizing identities $\mvec{\mv{A} \circ \mv{B} } = \mvec{\mv{A}
}\circ \mvec{\mv{B}}$ for $\mv{A}, \mv{B}~\in~\mathbb{C}^{l \times p}$
and $\mvec{\mv{A}\mv{B}} = \left(\mv{I}_m \otimes \mv{A} \right) \mvec{\mv{B}} =
\left(\mv{B}^\T \otimes \mv{I}_k \right) \mvec{\mv{A}}$, for $\mv{A} \in
\mathbb{C}^{k \times l}$ and $\mv{B} \in \mathbb{C}^{l \times m}$ and

\begin{equation}
- ( \mv{I}_r \otimes \mv{D}_{xx}) \mvec{\mv{U}}  - \left( (\mv{V}^\H \mv{D}_{yy} \mv{V})^\T \otimes \mv{I}_r  \right)\mvec{\mv{U}} - \mvec{\left(\mv{K} \circ \mv{U}\mv{V}^\H \right) \mv{V} }  = \mvec{\mv{F}\mv{V}}.
\end{equation}

We can simplify the last term of the left hand side and write

\begin{equation}\label{eq:vecThirdTerm}
\begin{aligned}
\mvec{  \left(\mv{K} \circ \mv{U}\mv{V}^\H \right) \mv{V} } &= \left(\mv{V}^\T \otimes \mv{I} \right) \mvec{\mv{K} \circ \mv{U} \mv{V}^\H},\\
&= \left(\mv{V}^\T \otimes \mv{I} \right) \left( \mvec{\mv{K}} \circ \mvec{\mv{U} \mv{V}^\H} \right),\\
&= \left(\mv{V}^\T \otimes \mv{I} \right) \diag{\mvec{\mv{K}}}  \mvec{\mv{U} \mv{V}^\H},\\
&= \left(\mv{V}^\T \otimes \mv{I} \right) \diag{\mvec{\mv{K}}}  \left(\left(\mv{V}^\H\right)^\T \otimes \mv{I} \right) \mvec{\mv{U}}.\\
\end{aligned}
\end{equation}

This results in

\begin{equation}\label{eq:first}
\left[ - ( \mv{I} \otimes \mv{D}_{xx})  - \left( \left(\mv{V}^\H \mv{D}_{yy} \mv{V}\right)^\T \otimes \mv{I}  \right) -  \left(\mv{V}^\T \otimes \mv{I} \right) \diag{\mvec{\mv{K}}}  \left(\left(\mv{V}^\H\right)^\T \otimes \mv{I} \right)\right] \mvec{\mv{U}} = \mvec{\mv{F}\mv{V}}.
\end{equation}

This is a linear system for the remaining unknown columns of the matrix $\mv{U} \in \mathbb{C}^{n \times r}$.

In (<ref>), we have approximated $\mv{A}$ as
$\mv{U}\mv{\Sigma}\mv{V}^\H$ where $U,V ∈ℂ^n
×r$and $Σ ∈ℂ^r ×r$ are
truncated matrices. With an orthogonal guess for $V$, we solve
for a $U$ in \eqref{eq:first}.  Since we approximate $A$
now by the product $UV^$̋, hence $\mv{U}$, solution of
(<ref>), includes the diagonal matrix with singular values.

We can now do a QR decomposition of $\mv{U}$ to arrive at a guess for
the orthogonal matrix $\mv{U}$.

The next step is to improve the guess for $\mv{V}$ in a similar way. The
equation (<ref>) becomes, when we multiply
from the left by $\mv{U}^\H$,
\begin{equation}
- \mv{U}^\H \mv{D}_{xx} \mv{U}\mv{V}^\H - \mv{U}^\H\mv{U}\mv{V}^\H \mv{D}_{yy} - \mv{U}^\H \left(\mv{K} \circ \mv{U}\mv{V}^\H\right) = \mv{U}^\H \mv{F}.
\end{equation}
Using the vectorizing identities this results in
\begin{equation}\label{eq:second}
\left[ - \left(\mv{I} \otimes  \mv{U}^\H \mv{D}_{xx}\mv{U}\right) - \left(\mv{D}_{yy}^\T \otimes \mv{I} \right) -  \left( \mv{I} \otimes \mv{U}^\H\right) \diag{\mvec{\mv{K}}}  \left( \mv{I} \otimes \mv{U} \right) \right] \mvec{\mv{V}^\H}   = \mvec{\mv{U}^\H \mv{F}},
\end{equation}
where we use that
\begin{equation}
\begin{aligned}
\mv{U}^\H \left(  \mv{K} \circ \mv{U}\mv{V}^\H\right)  & = \left( \mv{I} \otimes \mv{U}^\H \right) \mvec{\mv{K}} \circ \mvec{\mv{U}\mv{V}^\H} \\
& = \left( \mv{I} \otimes \mv{U}^\H \right) \diag{\mvec{\mv{K}}}  \left( \mv{I} \otimes \mv{U} \right) \mvec{\mv{V}^\H}. \\
\end{aligned}
\end{equation}

Combining the equations \eqref{eq:first} and \eqref{eq:second} we can
now propose an algorithm that updates $U$ and $V$ in an
alternating way. The steps are described in the

\begin{algorithm}[t]
\SetAlgoLined
Choose $\mv{V} \in \mathbb{C}^{n \times r}$ as initial guess\;
$[\mv{V},\mv{R}] = \qr{\mv{V},0}$\;
\While{not converged}{
Solve $\left[ - \left( \mv{I} \otimes \mv{D}_{xx}\right)  - \left( \left(\mv{V}^\H \mv{D}_{yy} \mv{V}\right)^\T \otimes \mv{I}  \right) -  \left(\mv{V}^\T \otimes \mv{I} \right) \diag{\mvec{\mv{K}}}  \left(\left(\mv{V}^\H\right)^\T \otimes \mv{I} \right)\right] \mvec{\mv{U}} = \mvec{\mv{F}\mv{V}}$\;
$[\mv{U},\mv{R}] = \qr{\mv{U},0}$\;
Solve $\left[ - \left(\mv{I} \otimes  \mv{U}^\H \mv{D}_{xx} \mv{U} \right) - \left(\mv{D}_{yy}^\T \otimes \mv{I}\right)  -  \left( \mv{I} \otimes \mv{U}^\H\right) \diag{\mvec{\mv{K}}}  \left( \mv{I} \otimes \mv{U} \right) \right] \mvec{\mv{V}^\H}   = \mvec{\mv{U}^\H \mv{F}}$\;
$[\mv{V},\mv{R}] = \qr{\mv{V},0}$\;
$\mv{A} = \mv{U} \mv{R}^\H \mv{V}^\H$\;
\caption{Solve for the low-rank matrix decomposition of the solution $\mv{A} \approx \mv{U}\mv{V}^\H$ of a 2D Helmholtz problem with space-dependent wave number.}
\label{alg:2d}
\end{algorithm}

\subsection{Comparison between coupled channel and a low-rank decomposition}
We now compare the coupled channel approach from section
\ref{sec:singleionization} with the low-rank approach from the
previous section, Sec.~\ref{sec:lowrankcomponents}.  In the coupled channel
calculation we use the eigenfunctions of one-particle subsystems as a
basis for the expansion. After integration over one of the coordinates
this results in a coupled set of one-dimensional equations. Let us
illustrate that equation \eqref{eq:first} and \eqref{eq:second} reduce
to a the coupled channel equations when we choose the eigenfunctions,
from \eqref{eq:eigenstate_H1} and \eqref{eq:eigenstate_H2}, as columns
for $\mv{U}$.

Let use take a look at the equation \eqref{eq:second} and assume that
$\mv{K} = E\mv{I} -V_1(\mv{x})-V_2(\mv{y})-V_{12}(\mv{x},\mv{y})$ and
that the columns of $\mv{U}$ are the eigenstates of
$-\mv{D}_{xx}+V_1(\mv{x})$.  We can then write
\begin{equation*}
- \left( \mv{I} \otimes \mv{U}^\H\right) \diag{\mvec{\mv{K}}}  \left( \mv{I} \otimes \mv{U} \right)
=  - E\,  \mv{I} \otimes \mv{I} + V_2(\mv{y}) \otimes \mv{I}  + \mv{I} \otimes \mv{U}^\H V_1(\mv{x}) \mv{U} + \left( \mv{I} \otimes \mv{U}^\H\right) \diag{ \mvec{V_{12}(\mv{x},\mv{y})}}  \left( \mv{I} \otimes \mv{U} \right).
\end{equation*}
Then equation \eqref{eq:second} becomes
\begin{equation*}
\left[  \mv{I} \otimes  \mv{U}^\H \left(-\mv{D}_{xx} + V_1(\mv{x})\right) \mv{U}  +  \left(-\mv{D}_{yy}^\T + V_2(\mv{y}) -E \mv{I}\right) \otimes \mv{I}  -   \left( \mv{I} \otimes \mv{U}^\H\right) \diag{ \mvec{V_{12}(\mv{x},\mv{y})}} \left( \mv{I} \otimes \mv{U} \right)  \right] \mvec{\mv{V}^\H}   = \mvec{\mv{U}^\H \mv{F}}.
\end{equation*}
When we use that the columns of $\mv{U}$ are eigenfunctions of
$-\mv{D}_{xx}+V_1(\mv{x})$ with eigenvalues $\mu_i$, equation
\eqref{eq:second} becomes
\begin{equation}
\left[  \mv{I} \otimes \diag{\mv{\mu}} + \left(-\mv{D}_{yy}^\T + V_2(\mv{y}) -E \, \mv{I}\right) \otimes \mv{I} -  \left( \mv{I} \otimes \mv{U}^\H\right) \diag{\mvec{V_{12}(\mv{x}, \mv{y})}}  \left( \mv{I} \otimes \mv{U} \right)  \right] \mvec{\mv{V}^\H}   = \mvec{\mv{U}^\H \mv{F}}.
\end{equation}
The term $\left( \mv{I} \otimes \mv{U}^\H\right)
\diag{\mvec{V_{12}(\mv{x}, \mv{y})}} \left( \mv{I} \otimes \mv{U}
\right)$ couples the columns of $\mv{V}$. It should be interpreted as
a discretized version of $\int \phi^*_i(x)V(x,y)\phi_j(x) dx$, where
we integrate over one of the coordinates. The columns of $\mv{U}$ are
the $\phi_i$ represented on a integration grid.

However, in general, the columns are $\mv{U}$ are not eigenfunctions of the
operator. The term $\mv{I} \otimes \diag{\mv{\mu}}$ then becomes a matrix that
also couples the different components of $\mv{V}$.

In short, each iteration we are solving a generalized coupled channel
\subsection{Convergence with projection operators}
\label{sec:projection2d}
We will now write both linear systems, \eqref{eq:first} for
$\mvec{\mv{U}}$, and, \eqref{eq:second} for $\mvec{\mv{V}}$, as
projection operators applied to the residual of the matrix equation,
\eqref{eq:matrixform}.

We denote by $\mv{L}$ the discretized 2D Helmholtz operator on the full grid
\begin{equation} \label{eq:linearoperator}
\mv{L} = \left( \mv{I} \otimes (-\mv{D}_{xx}) \right) + \left( (-\mv{D}_{yy}) \otimes \mv{I} \right) - \diag{\mvec{\mv{K}}} \left(\mv{I} \otimes \mv{I} \right).
\end{equation}
We can now explicitly write equation \eqref{eq:first} in terms of projections and this linear operator:
\begin{equation} \label{eq:solutionU}
\left(\mv{V}^\T \otimes \mv{I} \right) \mv{L} \left(\conj{\mv{V}} \otimes \mv{I} \right) \mvec{\mv{U}} = \left(\mv{V}^\T \otimes \mv{I} \right)\mvec{\mv{F}}.
\end{equation}

The residual matrix, $\mv{R}$, is given by
\begin{equation}
\begin{split}
\mv{R} &= \mv{F} - (-\mv{D}_{xx}) \mv{U} \mv{V}^\H - \mv{U} \mv{V}^\H (-\mv{D}_{yy}) + \mv{K} \circ \left(\mv{U}\mv{V}^\H\right).
\end{split}
\end{equation}
In vector form this reads, using that $\mv{U}$ is a solution of equation \eqref{eq:solutionU},
\begin{equation}
\begin{split}
\mvec{\mv{R}} &= \mvec{\mv{F}} - \left( \conj{\mv{V}} \otimes (- \mv{D}_{xx}) + (-\mv{D}_{yy}^\T) \conj{\mv{V}} \otimes \mv{I} - \diag{\mvec{\mv{K}}} \left( \conj{\mv{V}} \otimes \mv{I} \right)  \right) \mvec{\mv{U}}, \\
&= \left(\mv{I} - \left( \conj{\mv{V}} \otimes (-\mv{D}_{xx}) + (-\mv{D}_{yy}^\T) \conj{\mv{V}} \otimes \mv{I} - \diag{\mvec{\mv{K}}} \left( \conj{\mv{V}} \otimes \mv{I} \right)  \right) \left[ \left(\mv{V}^\T \otimes \mv{I} \right) \mv{L} \left(\conj{\mv{V}} \otimes \mv{I} \right) \right]^{-1} \left(\mv{V}^\T \otimes \mv{I} \right) \right) \mvec{\mv{F}} \\
&= \left(\mv{I} - \mv{L} \left( \conj{\mv{V}} \otimes \mv{I} \right) \left[ \left(\mv{V}^\T \otimes \mv{I} \right) \mv{L} \left(\conj{\mv{V}} \otimes \mv{I} \right) \right]^{-1} \left(\mv{V}^\T \otimes \mv{I} \right) \right) \mvec{\mv{F}} \\
&= P_{\mv{V}}\mvec{\mv{F}},
\end{split}
\end{equation}
where $P_{\mv{V}}$ is given by
\begin{equation}\label{eq:projectorVarOnV}
\begin{split}
P_{\mv{V}} &:= \mv{I} - \mv{L} \left( \conj{\mv{V}} \otimes \mv{I} \right) \left[ \left(\mv{V}^\T \otimes \mv{I} \right) \mv{L} \left(\conj{\mv{V}} \otimes \mv{I} \right) \right]^{-1} \left(\mv{V}^\T \otimes \mv{I} \right) \\
&:= \mv{I} - \mv{X}.
\end{split}
\end{equation}
The operator $P_{\mv{V}}$ is a projection operator. Indeed, observe
that the terms between the two inverses cancel, in the next equation,
against one of the inverse factors:
\begin{equation}
\begin{split}
\mv{X}^2 &= \mv{L} \left( \conj{\mv{V}} \otimes \mv{I} \right) \left[ \left(\mv{V}^\T \otimes \mv{I} \right) \mv{L} \left(\conj{\mv{V}} \otimes \mv{I} \right) \right]^{-1} \left(\mv{V}^\T \otimes \mv{I} \right) \mv{L} \left( \conj{\mv{V}} \otimes \mv{I} \right) \left[ \left(\mv{V}^\T \otimes \mv{I} \right) \mv{L} \left(\conj{\mv{V}} \otimes \mv{I} \right) \right]^{-1} \left(\mv{V}^\T \otimes \mv{I} \right) \\
&= \mv{L} \left( \conj{\mv{V}} \otimes \mv{I} \right) \left[ \left(\mv{V}^\T \otimes \mv{I} \right) \mv{L} \left(\conj{\mv{V}} \otimes \mv{I} \right) \right]^{-1} \left(\mv{V}^\T \otimes \mv{I} \right)\\
&= \mv{X}.
\end{split}
\end{equation}
We then have that $P_{\mv{V}}^2=(1-\mv{X})(1-\mv{X}) =
1-2\mv{X}+\mv{X}^2 = 1- \mv{X} = P_{\mv{V}}$.

This projection operator removes all components from the residual matrix
that can be corrected by the subspace spanned by $\mv{V}$. It is
similar as a deflation operator, often used in preconditioning

A similar derivation results in a projection operator $P_{\mv{U}}$ for
the update of $\mv{V}$:
\begin{equation}\label{eq:projectorVarOnU}
\begin{split}
P_{\mv{U}} &= \mv{I} - \mv{L} \left( \mv{I} \otimes \mv{U} \right) \left[\left(\mv{I} \otimes \mv{U}^\H \right) \mv{L} \left(\mv{I} \otimes \mv{U} \right)\right]^{-1} \left(\mv{I} \otimes \mv{U}^\H \right) \\
&= \mv{I} - \mv{Y}.
\end{split}
\end{equation}
This is again a projection operator.

So, Algorithm~\ref{alg:2d} repeatedly projects the residual matrix, $\mv{R}$, on a
subspace. Alternating between a subspace that is orthogonal to the
subspace spanned by the columns of $\mv{V}^{(k)}$ at iteration $k$ and
a subspace that is orthogonal to the columns of $\mv{U}^{(k)}$ at
iteration $k$. The residual matrix $\mv{R}^{(k)}$ after $k$ iterations
is the result of a series of projections
\begin{equation}
\mv{R}^{(k)} = P_{\mv{U}^{(k)}}P_{\mv{V}^{(k)}} \, P_{\mv{U}^{(k-1)}}P_{\mv{V}^{(k)}} \ldots \,P_{\mv{U}^{(0)}}P_{\mv{V}^{(0)}} \,\mv{R}^{(0)}.
\end{equation}
It is similar to the \textit{method of Alternating Projections}
[28] that goes back to Neumann
[15], where a solution is projected on two
alternating subspaces resulting in a solution that lies in the
intersection between the two spaces.  However, here the columns for
$\mv{U}^{(k)}$ and $\mv{V}^{(k)}$ are changing each iteration.

However, when the rank of $\mv{U}^{(k)}$ and $\mv{V}^{(k)}$ is
sufficiently large, the only intersection between the changing
subspaces  is the $\mv{0}$ matrix.
So the residual converges to $0$.

\section{Solving for the low-rank tensor approximation of 3D Helmholtz equations}
\label{sec:3dhelmholtz}
Algorithm~\ref{alg:2d}, as introduced in section~\ref{sec:lowrankcomponents}
for two dimensional problems, can be extended to higher dimensions.
We illustrate now this extension for 3D Helmholtz problems. First, we will
discuss the problem with a constant wave number and then extend the results
to space-dependent wave numbers.

We solve the Helmholtz equation with a constant or space-dependent wave number,
$k(x,y,z)$, in the first quadrant, where $x\ge 0, y \ge 0$ and $z \ge 0$. The equation is
\begin{equation}\label{eq:driven3d}
\left(  -\Delta_3 - k^2(x,y,z) \right) u_{sc}(x,y,z) = f(x,y,z),
\end{equation}
where $\Delta_3$ is the 3D Laplacian and the solution $u_{sc}$
satisfies homogeneous boundary conditions $u_{sc}(0,y,z)=0$ for all
$y,z\ge 0$, $u_{sc}(x,0,z)=0$ for all $x,z \ge 0$ and $u_{sc}(x,y,0)=0$
for all $x,y \ge 0$. On the other boundaries we have outgoing boundary conditions.

The right hand side $f(x,y,z)$ has a support that is limited to $[0,b]^3
\subset [0,L]^3 \subset \mathbb{R}_+^3$, i.e. $f(x,y,z)=0$, for all
$x \ge b$, $y \ge b$ or $z \ge b$.

The wave number $k(x,y,z)$ can be split in a constant part, $k_0^2$, and
variable part $\chi(x,y,z)$. The variable part is also only non-zero on
\begin{equation}
k^2(x,y,z) = \begin{cases}
k_0^2 \left(1  + \chi(x,y,z) \right)  \quad &\text{if} \quad  x<b \quad \text{and} \quad y<b \quad \text{and} \quad z<b,\\
k_0^2   \quad &\text{if} \quad  x\ge b \quad \text{or} \quad y\ge b \quad \text{or} \quad z\ge b.\\
\end{cases}
\end{equation}
We extend the domain with exterior complex scaling (ECS) absorbing
boundary condition [13].

The wave function is discretized on a three
dimensional mesh with $n_1 \times n_2 \times n_3$ unknowns
and can be represented by a Tucker tensor decomposition [25]
$\mt{M} \in \C^{n_1 \times n_2 \times n_3}$ with multi-linear rank
$\mv{r} = (r_1, r_2, r_3)$:
\begin{equation}
\mt{M} = \mt{G} \times_1 \mv{U}_1 \times_2 \mv{U}_2 \times_3 \mv{U}_3 \in \C^{n_1 \times n_2 \times n_3}.
\end{equation}
Here the tensor $\mt{G} \in \C^{r_1 \times r_2 \times r_3}$ is called
the core tensor and the factor matrices $\mv{U}_i \in \C^{n_i \times r_i}$
have orthonormal columns for $i = 1, 2, 3$.  Here
$r_i$ refers to the rank for each direction and $n_i$ to the number of
mesh points in each direction. So, to store this tensor only one core
tensor and $d$ factor matrices need to be stored, so the
storage costs\footnote{Here, we used the notation $r^d := \prod_{i=1}^{d} r_i$ and $r := \sqrt[d]{r^d}$,
to deal with possible unequal number of discretization points $n_i$ or ranks $r_i$ in different directions for a Tucker tensor.
%But, in the examples of this paper we will always use $r = r_1 = r_2 = r_3$.
scales $\Oh{r^d + dnr}$.

Let $\mathcal{L}$ be the discretization of
the three dimensional Helmholtz operator as given in
\eqref{eq:driven3d}. Observe that the operator $\mathcal{L}$ can be
written as a sum of Kronecker-products, where the matrix representation
$\mv{L}$ is of the following form
\begin{equation}\label{eq:operatorLasSumKroneckerProducts}
\mv{L} = -\mv{I} \otimes \mv{I} \otimes \mv{D}_{xx} - \mv{I} \otimes \mv{D}_{yy} \otimes \mv{I}  - \mv{D}_{zz} \otimes \mv{I} \otimes \mv{I} - \diag{\mvec{\mt{K}}} \mv{I} \otimes \mv{I} \otimes \mv{I}.
\end{equation}
Here $\mv{D}_{xx} \in \mathbb{C}^{n_1 \times n_1}, \mv{D}_{yy} \in \mathbb{C}^{n_2 \times n_2}$ and $\mv{D}_{zz} \in \mathbb{C}^{n_3 \times n_3}$ are sparse matrices that represent the discretization
of the second derivatives and $\mt{K}$ is the tensor that represents the
constant or space-dependent wave number, $k^2(x,y, z)$, discretized on the grid.

\subsection{Helmholtz equation with constant wave number}
First, consider the 3D Helmholtz problem with a constant wave number, so $k^2(x,y,z) \equiv k^2$.
The application of the Helmholtz operator $\mathcal{L}$ on
tensor $\mt{M}$ in Tucker tensor format is given by
\begin{equation}\label{eq:start3d}
\begin{split}
\mathcal{L} \mt{M} &= \mt{F} \\
\mathcal{L} \mt{M} &= -\mt{G} \times_1 \mv{D}_{xx} \mv{U}_1 \times_2 \mv{U}_2 \times_3 \mv{U}_3\\
&- \mt{G} \times_1 \mv{U}_1 \times_2 \mv{D}_{yy} \mv{U}_2 \times_3 \mv{U}_3\\
&- \mt{G} \times_1 \mv{U}_1 \times_2 \mv{U}_2 \times_3 \mv{D}_{zz} \mv{U}_3\\
&- k^2 \mt{G} \times_1 \mv{U}_1 \times_2 \mv{U}_2 \times_3 \mv{U}_3\\
&= \mt{F}
\end{split}
\end{equation}
where $\mv{U}_i^\H \mv{U}_i = \mv{I}$ for $i = 1,2,3$ and $\mt{F}$ is a tensor representation of the right hand side function $f$ discretized on the grid.

\subsubsection{Solving for product of basis functions and core terms (version 1)}\label{sec:version1}
Similar to the two dimensional case, we can derive equations to
iteratively solve for the factors $\mv{U}_1$, $\mv{U}_2$ and
$\mv{U}_3$. To derive the equations for $\mv{U}_1$ we start from
\eqref{eq:start3d} and multiply with $\mv{U}_2$ and $\mv{U}_3$ in the
second and third direction, respectively:
\begin{equation*}
\mathcal{L}\mt{M} \times_2 \mv{U}_2^\H \times_3 \mv{U}_3^\H = \mt{F} \times_2 \mv{U}_2^\H \times_3 \mv{U}_3^\H.
\end{equation*}
For a review of tensors, tensor decompositions and tensor
operations like this tensor-times-matrix product denoted by
the symbol $\times_i$ we refer to [12].

Using explicitly the Tucker tensor representation for $\mt{M}$ and that the columns
of $\mv{U}_i$ are orthonormal, the following expression is derived for
the Helmholtz operator applied on a tensor in Tucker format
\begin{equation*}
\mathcal{L}\mt{M} \times_2 \mv{U}_2^\H \times_3 \mv{U}_3^\H
= \mt{G} \times_1 (-\mv{D}_{xx} - k^2 \mv{I}) \mv{U}_1
- \mt{G} \times_1 \mv{U}_1 \times_2 \mv{U}_2^\H \mv{D}_{yy} \mv{U}_2
- \mt{G} \times_1 \mv{U}_1 \times_3 \mv{U}_3^\H \mv{D}_{zz} \mv{U}_3.
\end{equation*}
Writing this tensor equation in the first unfolding leads to a matrix
equation, recall the first unfolding is given by $\mv{M}_{(1)} =
\conj{\mv{U}_1} \mv{G}_{(1)} \left( \mv{U}_3 \otimes \mv{U}_2
\right)^\H$, see also [12]:
\begin{equation}\label{eq:firstUnfolding}
\conj{\left(-\mv{D}_{xx} - k^2 \mv{I} \right) \mv{U}_1} \mv{G}_{(1)}
- \conj{\mv{U}_1} \mv{G}_{(1)} \left(I \otimes \mv{U}_2^\H \mv{D}_{yy} \mv{U}_2\right)^\H
- \conj{\mv{U}_1} \mv{G}_{(1)} \left(\mv{U}_3^\H \mv{D}_{zz} \mv{U}_3  \otimes I\right)^\H
=  \mv{F}_{(1)} \left(\mv{U}_3 \otimes \mv{U}_2\right).
\end{equation}
To solve this equation for $\mv{U}_1$, it is written in vectorized form as
\begin{equation}\label{eq:eqnU1}
\left\{
\mv{I} \otimes \conj{\left(-\mv{D}_{xx} - k^2 \mv{I} \right)} +
\left[
-\conj{\left(\mv{I} \otimes \mv{U}_2^\H \mv{D}_{yy} \mv{U}_2\right)} - \conj{\left(\mv{U}_3^\H \mv{D}_{zz} \mv{U}_3 \otimes \mv{I} \right)}
\right] \otimes \mv{I}
\right\} \mvec{\underbrace{\conj{\mv{U}_1}\mv{G}_{(1)}}_{\mv{X}_1}} = \mvec{\mv{F}_{(1)} \left(\mv{U}_3 \otimes \mv{U}_2\right)}.
\end{equation}
Observe that this is a square system with $n_1 \times r_2 r_3$ unknowns,
where the solution in matrix form $\mv{X}_1$ could have, in general, a
rank $r > r_1$.
In a similar way, equations for $\mv{U}_2$ and $\mv{U}_3$ are
derived by multiplying \eqref{eq:start3d} with the other factor
matrices in the appropriate directions:
\begin{equation}\label{eq:eqnU2}
\left\{
\mv{I} \otimes \conj{\left(-\mv{D}_{yy} - k^2 \mv{I} \right)} +
\left[
\conj{-\left(\mv{I} \otimes \mv{U}_1^\H \mv{D}_{xx} \mv{U}_1\right)} - \conj{\left(\mv{U}_3^\H \mv{D}_{zz} \mv{U}_3 \otimes \mv{I} \right)}
\right]  \otimes \mv{I}
\right\} \mvec{\underbrace{\conj{\mv{U}_2}\mv{G}_{(2)}}_{\mv{X}_2}}
= \mvec{\mv{F}_{(2)} \left(\mv{U}_3 \otimes \mv{U}_1\right)},
\end{equation}
\begin{equation}\label{eq:eqnU3}
\left\{
\mv{I} \otimes \conj{\left(-\mv{D}_{zz} - k^2 \mv{I} \right)} +
\left[-
\conj{\left(\mv{I} \otimes \mv{U}_1^\H \mv{D}_{xx} \mv{U}_1\right)} - \conj{\left(\mv{U}_2^\H \mv{D}_{yy} \mv{U}_2 \otimes \mv{I} \right)}
\right] \otimes \mv{I}
\right\} \mvec{\underbrace{\conj{\mv{U}_3}\mv{G}_{(3)}}_{\mv{X}_3}}
= \mvec{\mv{F}_{(3)} \left(\mv{U}_2 \otimes \mv{U}_1\right)}.
\end{equation}
Alternating between solving for $\mv{U}_1$, $\mv{U}_2$ and $\mv{U}_3$
using \eqref{eq:eqnU1}, \eqref{eq:eqnU2} or \eqref{eq:eqnU3} results
in algorithm that approximates the low-rank solutions for three
dimensional problems as given in \eqref{eq:driven3d}. This algorithm
is summarized in Algorithm~\ref{alg:const3dv1}.
Also in the three dimensional case the orthogonality of the columns
of $\mv{U}_1$, $\mv{U}_2$ and $\mv{U}_3$ are maintained by
additional QR factorizations.
\begin{algorithm}[t]
\SetAlgoLined
[$\mt{G}, \mv{U}_1, \mv{U}_2, \mv{U}_3$] = hosvd(initial guess)\;
\While{not converged}{
\For{i = 1, 2, 3}{
Solve for $\mv{X}_i = \conj{\mv{U}_i}\mv{G}_{(i)} \in \C^{n \times r^{d-1}}$ using \eqref{eq:eqnU1}, \eqref{eq:eqnU2} or \eqref{eq:eqnU3}\;
$\conj{\mv{U}_{i}} \mv{G}_{(i)} = \qr{\mv{X}_i(:,~ 1:r_i), 0}$;
$\mt{G} = \texttt{reconstruct} \left[\mv{G}_{(i)}, ~i\right]$\;
$\mt{M} = \mt{G} \times_1 \mv{U}_1  \times_2 \mv{U}_2  \times_3 \mv{U}_3$\;
\caption{Solve for the low-rank tensor decomposition of the solution $\mt{M}$ of a 3D Helmholtz with constant wave number (version 1).}
\label{alg:const3dv1}
\end{algorithm}
Observe that we solve for a large matrix $\mv{X}_i \in \C^{n_i \times
r_1r_2r_3/r_i}$. So, in general the rank of this matrix could be
$\min\left(n_i,~ r_1r_2r_3/r_i\right)$. But it is also known that
$\mv{X}_i = \conj{\mv{U}_i}\mv{G}_{(i)}$ which leads to the fact that
the rank of $\mv{X}_i$ should be at most $r_i$. Selecting the
first $r_i$ columns of $\mv{X}_i$ and computing its QR decomposition
is sufficient to derive a new orthonormal basis as factor matrix

Finally, observe that solving for $\mv{X}_i$ using \eqref{eq:eqnU1},
\eqref{eq:eqnU2} or \eqref{eq:eqnU3} is computationally not
efficient. In all iterations, we solve for a total of $d nr^{d-1}$
unknowns, while there are only $r^d + dnr$ unknowns in the Tucker
tensor factorization. Furthermore, solving equations \eqref{eq:eqnU1},
\eqref{eq:eqnU2} and \eqref{eq:eqnU3} is also expensive. Indeed,
computing a symmetric reverse Cuthill-McKee permutation of the system
matrix one observes a matrix with a bandwidth $\Oh{r^{d-1}}$. For
example when $d = 3, n_i = n = 168, r_i = r = 18$ one obtains the
sparsity pattern on the diagonal of the matrix as shown in Figure
\ref{fig:spytop_symrcm_UGsystem_alg1}. So solving a system as given in
\eqref{eq:eqnU1}, \eqref{eq:eqnU2} or \eqref{eq:eqnU3} has a
computational cost of $\Oh{nr^{2(d-1)}}$.

\begin{figure}[t]
\centering
\begin{subfigure}{0.4\textwidth}
\includegraphics[width=\textwidth]{spytop_symrcm_UGsystem_alg1.png}
\caption{Sparsity pattern of the top of the symmetric reverse Cuthill-McKee permutation of the system matrix to solve for $\mv{X}_i$ using \eqref{eq:eqnU1}. Note: only the first 4.5 of the 168 blocks are shown.}
\label{fig:spytop_symrcm_UGsystem_alg1}
\end{subfigure}
\quad
\begin{subfigure}{0.4\textwidth}
\includegraphics[width=\textwidth]{spytop_symrcm_Gsystem_alg2.png}
\caption{Sparsity pattern of the symmetric reverse Cuthill-McKee permutation of the system matrix to solve for $\mv{G}_{(1)}$ using \eqref{eq:eqnGv2}.}
\label{fig:spytop_symrcm_Gsystem_alg2}
\end{subfigure}
\caption{Sparsity patterns of the symmetric reverse Cuthill-McKee permutation of certain system matrices ($d = 3, n = 168, r = 18$).}
\end{figure}

\subsubsection{Solving for the basis functions and the core tensor separately (version 2)}\label{sec:version2}
To circumvent solving the large systems in \eqref{eq:eqnU1},
\eqref{eq:eqnU2} and \eqref{eq:eqnU3}, we can pre-compute the
QR factorization of the unfolding of the core tensor, $\mv{G}_{(i)}$,
and project the equations onto the obtained $\mv{Q}_i$.  Indeed, this
will further reduce the number of unknowns in these linear systems to exactly
the number of unknowns that are needed for the factor matrices
$\conj{\mv{U}}_i$, for $i = 1,2,3$.

Let us discuss the details. We start again from equation \eqref{eq:firstUnfolding} and use the QR
factorization of $\mv{G}_{(1)}^\H$, $\mv{Q}_1 \mv{R}_1^\H = \qr{\mv{G}^\H_{(1)}}$. This yields
\begin{equation*}
\conj{\left(-\mv{D}_{xx} - k^2 \mv{I} \right) \mv{U}_1} \mv{R}_1\mv{Q}_1^\H
- \conj{\mv{U}_1} \mv{R}_1\mv{Q}_1^\H \left(I \otimes \mv{U}_2^\H \mv{D}_{yy} \mv{U}_2\right)^\H
- \conj{\mv{U}_1} \mv{R}_1\mv{Q}_1^\H \left(\mv{U}_3^\H \mv{D}_{zz} \mv{U}_3  \otimes I\right)^\H
=  \mv{F}_{(1)} \left(\mv{U}_3 \otimes \mv{U}_2\right).
\end{equation*}
Post multiplication of this equation by $\mv{Q}_1$ yields
\begin{equation*}
\conj{\left(-\mv{D}_{xx} - k^2 \mv{I} \right) \mv{U}_1} \mv{R}_1
- \conj{\mv{U}_1} \mv{R}_1\mv{Q}_1^\H \left(I \otimes \mv{U}_2^\H \mv{D}_{yy} \mv{U}_2\right)^\H \mv{Q}_1
- \conj{\mv{U}_1} \mv{R}_1\mv{Q}_1^\H \left(\mv{U}_3^\H \mv{D}_{zz} \mv{U}_3  \otimes I\right)^\H \mv{Q}_1
=  \mv{F}_{(1)} \left(\mv{U}_3 \otimes \mv{U}_2\right)\mv{Q}_1.
\end{equation*}
To solve this equation for $\mv{U}_1$, it is written is vectorized form
\begin{equation}\label{eq:eqnU1v2}
\left\{
\mv{I} \otimes \conj{\left(-\mv{D}_{xx} - k^2 \mv{I} \right)} +
\mv{Q}_1^\T\left[ -\conj{\left(\mv{I} \otimes \mv{U}_2^\H \mv{D}_{yy} \mv{U}_2\right)} - \conj{\left(\mv{U}_3^\H \mv{D}_{zz} \mv{U}_3 \otimes \mv{I} \right)}
\right]\conj{\mv{Q}_1} \otimes \mv{I}
\right\} \mvec{\underbrace{\conj{\mv{U}_1}\mv{R}_1}_{\mv{X}_1}} = \mvec{\mv{F}_{(1)} \left(\mv{U}_3 \otimes \mv{U}_2\right) \mv{Q}_1}.
\end{equation}
In a similar way, the update equations for $\mv{U}_2$ and $\mv{U}_3$ are
derived by multiplying \eqref{eq:start3d} with the other factor
matrices in the appropriate dimensions and using the QR factorizations
of $\mv{G}_{(i)}^\H$:
\begin{equation}\label{eq:eqnU2v2}
\left\{
\mv{I} \otimes \conj{\left(-\mv{D}_{yy} - k^2 \mv{I} \right)} +
\mv{Q}_2^\T\left[-
\conj{\left(\mv{I} \otimes \mv{U}_1^\H \mv{D}_{xx} \mv{U}_1\right)} - \conj{\left(\mv{U}_3^\H \mv{D}_{zz} \mv{U}_3 \otimes \mv{I} \right)}
\right]\conj{\mv{Q}_2}  \otimes \mv{I}
\right\} \mvec{\underbrace{{\conj{\mv{U}_2}\mv{R}_2}}_{\mv{X}_2}}
= \mvec{\mv{F}_{(2)} \left(\mv{U}_3 \otimes \mv{U}_1\right) \mv{Q}_2}.
\end{equation}
\begin{equation}\label{eq:eqnU3v2}
\left\{
\mv{I} \otimes \conj{\left(-\mv{D}_{zz} - k^2 \mv{I} \right)} +
\mv{Q}_3^\T\left[ -
\conj{\left(\mv{I} \otimes \mv{U}_1^\H \mv{D}_{xx} \mv{U}_1\right)} - \conj{\left(\mv{U}_2^\H \mv{D}_{yy} \mv{U}_2 \otimes \mv{I} \right)}
\right]\conj{\mv{Q}_3} \otimes \mv{I}
\right\} \mvec{\underbrace{\conj{\mv{U}_3}\mv{R}_3}_{\mv{X}_3}}
= \mvec{\mv{F}_{(3)} \left(\mv{U}_2 \otimes \mv{U}_1\right) \mv{Q}_3}.
\end{equation}
All these equations are cheap to solve. Indeed,
$\mvec{\conj{\mv{U}_i}\mv{R}_i}$ has length $n_i r_i$. Computing a
symmetric reverse Cuthill-McKee permutation of these system matrices one
observes a matrix with a bandwidth $\Oh{r}$, so solving these
equations has a computational cost $\Oh{nr^2}$.

Of course, this only updates the factor matrices as basis vectors in each
direction. As a single final step, we still have to compute the core
tensor $\mt{G}$. This will be the computationally most expensive

Core tensor $\mt{G}$ can be obtained by multiplying \eqref{eq:start3d}
with all the $d$ factor matrices in the matching directions.
Unfolding this equation in a certain direction (eg. the first folding)
leads again to a matrix equation. In vectorized form, it is given by
\begin{equation}\label{eq:eqnGv2}
\left\{
\mv{I} \otimes \conj{\mv{U}_1^\H\left(-\mv{D}_{xx} - k^2 \mv{I} \right)\mv{U}_1} +
\left[
-\conj{\left(\mv{I} \otimes \mv{U}_2^\H \mv{D}_{yy} \mv{U}_2\right)} - \conj{\left(\mv{U}_3^\H \mv{D}_{zz} \mv{U}_3 \otimes \mv{I} \right)}
\right] \otimes \mv{I}
\right\} \mvec{\mv{G}_{(1)}} = \mvec{\mv{U}_1^\T \mv{F}_{(1)} \left(\mv{U}_3 \otimes \mv{U}_2\right)}.
\end{equation}
Indeed, considering again an example where $d = 3, n_i = n = 168, r_i
= r = 18$ one obtains a matrix with a sparsity pattern that is shown
in  Figure~\ref{fig:spytop_symrcm_Gsystem_alg2}. Hence, this matrix has
not a limited bandwidth anymore. It coupled all functions to all other
functions. Although this equation has to be solved only once in the
algorithm, when the rank increases, it will rapidly dominate the computational
cost of this algorithm.

\begin{algorithm}[t]
[$\mt{G}, \mv{U}_1, \mv{U}_2, \mv{U}_3$] = hosvd(initial guess)\;
\While{not converged}{
\For{i = 1, 2, 3}{
$\mv{Q}_i \widetilde{\mv{R}} = \qr{\mv{G}_{(i)}^\H, 0}$\;
Solve for $\mv{X}_i = \conj{\mv{U}_i}\mv{R}_i \in \C^{n_i \times r_i}$ using \eqref{eq:eqnU1v2}, \eqref{eq:eqnU2v2} or \eqref{eq:eqnU3v2}\;
$\conj{\mv{U}_{i}} \mv{R}_i = \qr{\mv{X}_i, 0}$\;
$\mt{G} = \texttt{reconstruct} \left[\mv{R}_i\mv{Q}_i^\H, ~i\right]$;
Solve for $\mv{G}_{(1)} \in \C^{r_1 \times r_2r_3}$ using \eqref{eq:eqnGv2}\;
$\mt{G} = \texttt{reconstruct} \left[\mv{G}_{(1)}, ~1\right]$\;
$\mt{M} = \mt{G} \times_1 \mv{U}_1  \times_2 \mv{U}_2  \times_3 \mv{U}_3$\;
\caption{Solve for the low-rank tensor decomposition of the solution $\mt{M}$ of a 3D Helmholtz problem with constant wave number (version 2).}
\label{alg:const3dv2}
\end{algorithm}

\subsubsection{Efficient combination of version 1 and version 2 into new algorithm (version 3)}\label{sec:version3}
In the first version of the algorithm, (see section~\ref{sec:version1}),
an update for $\mv{G}_{(i)}$ is computed for each direction in each
iteration. This leads to a too expensive algorithm.
Then we changed the algorithm such that the costs for the updates in each
direction is reduced, (see section~\ref{sec:version2}). But, in that version
almost all information for a full update of core tensor $\mt{G}$ is
lost. Therefore a final, but potential too expensive, equation needs to be solved.

Observe that the expensive computation for the full core tensor, in
version 2, can now be replaced by a single solve per iteration as done in
version 1. This leads to a third version of the algorithm.
It avoids repeatedly solving the large systems (like version 1) and it
does not solve too expensive systems (like version 2). The computational
complexity of this algorithm is equal to the complexity of version 1,
so $\Oh{nr^{2(d-1)}}$.
Furthermore, the systems that need to be solved, each iteration, have
exactly the same number of unknowns as the representation of the
tensor in low-rank Tucker tensor format.
In summary, this final version of the algorithm is given by Algorithm~\ref{alg:const3dv3}.

\begin{algorithm}[t]
[$\mt{G}, \mv{U}_1, \mv{U}_2, \mv{U}_3$] = hosvd(initial guess)\;
\While{not converged}{
\For{i = 1, 2}{
$\mv{Q}_i \widetilde{\mv{R}} = \qr{\mv{G}_{(i)}^\H, 0}$\;
Solve for $\mv{X}_i = \conj{\mv{U}_i}\mv{R}_i \in \C^{n_i \times r_i}$ using \eqref{eq:eqnU1v2} or \eqref{eq:eqnU2v2}\;
$\conj{\mv{U}_{i}} \mv{R}_i = \qr{\mv{X}_i, 0}$\;
$\mt{G} = \texttt{reconstruct} \left[\mv{R}_i\mv{Q}_i^\H, ~i\right]$;
Solve for $\mv{X}_3 = \conj{\mv{U}_3}\mv{G}_{(3)} \in \C^{n_3 \times r^{d-1}}$ using \eqref{eq:eqnU3}\;
$\conj{\mv{U}_{3}} \mv{G}_{(3)} = \qr{\mv{X}_3, 0}$\;
$\mt{G} = \texttt{reconstruct} \left[\mv{G}_{(3)},~3\right]$\;
$\mt{M} = \mt{G} \times_1 \mv{U}_1  \times_2 \mv{U}_2  \times_3 \mv{U}_3$\;
\caption{Solve for the low-rank tensor decomposition of the solution $\mt{M}$ of a 3D Helmholtz problem with constant wave number (version 3).}
\label{alg:const3dv3}
\end{algorithm}

\subsubsection{Numerical comparison of three versions for 3D Helmholtz equation}
Consider a three dimensional domain $\Omega = [-10, ~10]^3$ that is
discretized with $M=100$ equidistant mesh points per direction in the
interior of the domain. The domain is extended with exterior complex
scaling to implement the absorbing boundary conditions. Hence, in total
there are $n = n_x = n_y = n_z = 168$ unknowns per direction. As
constant wave number we use $\omega = 2$ and a right hand side $f(x,y,z)
= -e^{-x^2-y^2-z^2}$.
By symmetry, we expect a low rank factorization with a equal low-rank
in each direction, so we fix $r = r_x = r_y = r_x$.

The convergence of the residuals of the three versions are given in
the left column of Figure~\ref{fig:residual+runtime-const3d}.
It is clear that all three versions converge to a good low-rank
approximation of the full solution. By increasing the maximal
attainable rank $r$, a better low-rank solution is obtained, as
expected. Remarkably,  for $r=30$, in version 2, the final residual
is larger then the residuals obtained by both other algorithms while
the compute-time for version 2 is larger than the other algorithms.

\begin{figure}
\centering
\input{residual-v1-M=100.tikz} ~ \input{runtime-v1-3d-M=100.tikz} \par
\input{residual-v2-M=100.tikz} ~ \input{runtime-v2-3d-M=100.tikz} \par
\input{residual-v3-M=100.tikz} ~ \input{runtime-v3-3d-M=100.tikz}
\caption{Left: Plot of residual per iteration for constant wave number in 3D Helmholtz problem.
Right: Plot of runtime of most time consuming parts for constant wave number in 3D Helmholtz problem. Both problems have $M=100$.
Top: Algorithm~\ref{alg:const3dv1} (version 1),
middle: Algorithm~\ref{alg:const3dv2} (version 2),
bottom: Algorithm~\ref{alg:const3dv3} (version 3).}
\label{fig:residual+runtime-const3d}
\end{figure}

The compute-time for the most time-consuming parts in the different
versions of the algorithm can be measured as a function of the maximal
attainable rank $r$. For the three versions of the algorithm the
runtimes are shown in the right column of
Figure~\ref{fig:residual+runtime-const3d}. For all parts the expected
and measured dependence on the rank $r$ are given. For all versions of
the algorithm 10 iterations are applied.

\begin{figure}
\centering
\input{runtimes-v123-M=100-iter=4.tikz}
\caption{Plot of runtime for 4 iteration with constant wavenumber in 3D using the three different algorithms ($M=100$)
\label{fig:runtime-const3d}}
\end{figure}

Comparing the total runtime for the three different versions one obtains
results as shown in Figure~\ref{fig:runtime-const3d}. Indeed, as
expected version 3 is approximately 3 times faster than version 1 and
the runtime scales similar in rank $r$. Further, for small rank $r$
version 2 is faster than both other versions. But when the rank
increases the expensive solve for the core tensor $\mt{G}$ starts to
dominate the runtime. The total runtime will increase dramatically.

\subsection{Projection operator for constant wave number}
Also in three dimensions we can write the linear systems
\eqref{eq:eqnU1} for $\mv{U}_1$, \eqref{eq:eqnU2} for
$\mv{U}_2$ and \eqref{eq:eqnU3} for $\mv{U}_3$ as
projection operators applied to the residual of the tensor equation,
\eqref{eq:start3d}.

Consider a tensor $\mt{M}$ in Tucker format and factorized as $\mt{M}
= \mt{G} \times_1 \mv{U}_1 \times_2 \mv{U}_2 \times_3 \mv{U}_3$, with
unknowns $\mt{G}, \mv{U}_1, \mv{U}_2$ and $\mv{U}_3$. Discretization
of \eqref{eq:driven3d} leads to a linear operator $\mathcal{L}$
applied on tensors. Its matrix representation $\mv{L}$ has a sum
of Kronecker products structure, as given in
\eqref{eq:operatorLasSumKroneckerProducts}.

Solving for an unknown factors $\mv{U}_1$, $\mv{U}_2$ or $\mv{U}_3$
(and the core-tensor $\mt{G}$) using \eqref{eq:eqnU1},
\eqref{eq:eqnU2} or \eqref{eq:eqnU3} can be interpreted as a
projection operator applied on the residual.  For example,
\eqref{eq:eqnU1} can be interpreted as
\begin{equation}
\conj{\left(\mv{U}_3^\H \otimes \mv{U}_2^\H \otimes \mv{I} \right) \mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)} \mvec{\conj{\mv{U}_1}\mv{G}_{(1)}} = \left( \mv{U}_3^\T \otimes \mv{U}_2^\T \otimes \mv{I} \right) \mvec{\mv{F}_{(1)}}.
\end{equation}
The residual, in tensor format, is given by
\begin{equation}
\begin{split}
\mt{R} &= \mt{F} - \mathcal{L}\mt{M}, \\
&= \mt{F} - \mt{G} \times_1 \left(-\mv{D}_{xx} - k^2 \mv{I}\right) \mv{U}_1 \times_2 \mv{U}_2 \times_3 \mv{U}_3 + \mt{G} \times_1 \mv{U}_1 \times_2 \mv{D}_{yy}\mv{U}_2 \times_3 \mv{U}_3 + \mt{G} \times_1 \mv{U}_1 \times_2 \mv{U}_2 \times_3 \mv{D}_{zz}\mv{U}_3.
\end{split}
\end{equation}
Writing this tensor equation in the first unfolding leads to the following matrix equation
\begin{equation}
\mv{R}_{(1)} = \mv{F}_{(1)} - \conj{\left(-\mv{D}_{xx} - k^2 \mv{I}\right)\mv{U}_1}\mv{G}_{(1)} \left(\mv{U}_3 \otimes \mv{U}_2 \right)^\H + \conj{\mv{U}_1}\mv{G}_{(1)} \left(\mv{U}_3 \otimes \mv{D}_{yy}\mv{U}_2 \right)^\H + \conj{\mv{U}_1}\mv{G}_{(1)} \left(\mv{D}_{zz}\mv{U}_3 \otimes \mv{U}_2 \right)^\H,
\end{equation}
which can be vectorized as
\begin{equation}
\begin{split}
\mvec{\mv{R}_{(1)}} &= \mvec{\mv{F}_{(1)}} - \left(\conj{\left(\mv{U}_3 \otimes \mv{U}_2 \otimes \left(-\mv{D}_{xx}-k^2\mv{I}\right)\right)} - \conj{\left(\mv{U}_3 \otimes \mv{D}_{yy}\mv{U}_2 \otimes \mv{I}\right)} - \conj{\left(\mv{D}_{zz}\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)}\right)\mvec{\conj{\mv{U}_1}\mv{G}_{(1)}} \\
&= \mvec{\mv{F}_{(1)}} - \conj{\mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)}\mvec{\conj{\mv{U}_1}\mv{G}_{(1)}} \\
&= \mvec{\mv{F}_{(1)}} - \conj{\mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)} \left[\conj{\left(\mv{U}_3^\H \otimes \mv{U}_2^\H \otimes \mv{I} \right) \mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)}\right]^{-1}\left( \mv{U}_3^\T \otimes \mv{U}_2^\T \otimes \mv{I} \right) \mvec{\mv{F}_{(1)}} \\
&= P_{23}\mvec{\mv{F}_{(1)}},
\end{split}
\end{equation}
where operator $P_{23}$ is given by
\begin{equation}\label{eq:projector3dconst}
\begin{split}
P_{23} &= \mv{I} - \conj{\mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)} \left[\conj{\left(\mv{U}_3^\H \otimes \mv{U}_2^\H \otimes \mv{I} \right) \mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)}\right]^{-1}\left( \mv{U}_3^\T \otimes \mv{U}_2^\T \otimes \mv{I} \right) \\
&= \mv{I} - \mv{X}.
\end{split}
\end{equation}
This operator $P_{23}$ is indeed a projection operator. Observe that the terms between the two inverses cancel against one of the inverse factors:
\begin{equation*}
\begin{split}
\mv{X}^2 &= \conj{\mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)} \left[\conj{\left(\mv{U}_3^\H \otimes \mv{U}_2^\H \otimes \mv{I} \right) \mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)}\right]^{-1}\left( \mv{U}_3^\T \otimes \mv{U}_2^\T \otimes \mv{I} \right)\conj{\mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)} \left[\conj{\left(\mv{U}_3^\H \otimes \mv{U}_2^\H \otimes \mv{I} \right) \mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)}\right]^{-1}\left( \mv{U}_3^\T \otimes \mv{U}_2^\T \otimes \mv{I} \right) \\
&= \conj{\mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)} \left[\conj{\left(\mv{U}_3^\H \otimes \mv{U}_2^\H \otimes \mv{I} \right) \mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)}\right]^{-1}\left( \mv{U}_3^\T \otimes \mv{U}_2^\T \otimes \mv{I} \right) \\
&= \mv{X}.
\end{split}
\end{equation*}
This operator is a natural extension to higher dimensions of the two
dimensional operators as derived in section \ref{sec:projection2d}.

A similar derivation results in projection operators $P_{13}$ and $P_{12}$ for the updates in $\mv{U}_2$ and $\mv{U}_3$, respectively.
\begin{equation}\label{eq:projectors3dconst}
\begin{split}
P_{23} &= \mv{I} - \conj{\mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)} \left[\conj{\left(\mv{U}_3^\H \otimes \mv{U}_2^\H \otimes \mv{I} \right) \mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)}\right]^{-1}\left( \mv{U}_3^\T \otimes \mv{U}_2^\T \otimes \mv{I} \right), \\
P_{13} &= \mv{I} - \conj{\mv{L} \left(\mv{U}_3 \otimes \mv{I} \otimes \mv{U}_1\right)} \left[\conj{\left(\mv{U}_3^\H \otimes \mv{I} \otimes \mv{U}_1^\H \right) \mv{L} \left(\mv{U}_3 \otimes \mv{I} \otimes \mv{U}_1\right)}\right]^{-1}\left( \mv{U}_3^\T \otimes \mv{I} \otimes \mv{U}_1^\T \right), \\
P_{12} &= \mv{I} - \conj{\mv{L} \left(\mv{I} \otimes \mv{U}_2 \otimes \mv{U}_1\right)} \left[\conj{\left(\mv{I} \otimes \mv{U}_2^\H \otimes \mv{U}_1^\H \right) \mv{L} \left(\mv{I} \otimes \mv{U}_2 \otimes \mv{U}_1\right)}\right]^{-1}\left( \mv{I} \otimes \mv{U}_2^\T \otimes \mv{U}_1^\T \right).
\end{split}
\end{equation}
The successive application of these projection operators on the
residual results in an updated residual that lies in the intersection
of all subspaces.

\subsection{Helmholtz equation with space-dependent wave number}
The presented algorithms  with constant wave number can be
extended to space-dependent wave numbers. So, lets
consider a 3D Helmholtz problem  where $\mt{K} = k^2(x,y,z)$ represents the space-dependent
wave number on the discretized mesh.

Further, we assume that a Canonical Polyadic decomposition of the
space-dependent wave number tensor
$\mt{K}$ is known, i.e.
\begin{equation}\label{eq:CP4K}
\mt{K} = \sum_{i=1}^s \sigma_i \left(\mv{v}_i^{(1)} \circ \mv{v}_i^{(2)} \circ \cdots \circ	\mv{v}_i^{(d)}\right),
\end{equation}
where $s \in \N_+$ is the CP-rank of $\mt{K}$ and $\mv{v}_i^{(j)} \in \C^{n_j}$ for $i = 1,2,\ldots, s; j = 1,2,\ldots, d$ are vectors. Further, $\sigma_i$ is a tensor generalization of a singular value and $\circ$ denotes the vector outer product.

The application of the space-dependent Helmholtz operator $\mathcal{L}$ on tensor $\mt{M}$ is given by
\begin{equation}\label{eq:start3dvar}
\begin{split}
\mathcal{L} \mt{M} &= \mt{F} \\
\mathcal{L} \mt{M} &= -\mt{G} \times_1 \mv{D}_{xx} \mv{U}_1 \times_2 \mv{U}_2 \times_3 \mv{U}_3\\
&- \mt{G} \times_1 \mv{U}_1 \times_2 \mv{D}_{yy} \mv{U}_2 \times_3 \mv{U}_3\\
&- \mt{G} \times_1 \mv{U}_1 \times_2 \mv{U}_2 \times_3 \mv{D}_{zz} \mv{U}_3\\
&- \mt{K} \circ \left(\mt{G} \times_1 \mv{U}_1 \times_2 \mv{U}_2 \times_3 \mv{U}_3 \right)\\
&= \mt{F},
\end{split}
\end{equation}
where $\mv{U}_i^\H \mv{U}_i = \mv{I}$ for $i = 1,2,3$ and $\mt{F}$
is a tensor representation of the right hand side function $f$
discretized on the used grid.
Here $\circ$ denotes the Hadamard product for tensors.

In a similar way as in the three dimensional constant wave number
case, we can derive equations to iteratively solve for the factors
$\mv{U}_1$, $\mv{U}_2$ and $\mv{U}_3$. We start from \eqref{eq:start3dvar} and multiply with
$\mv{U}_2$ and $\mv{U}_3$ in the second and third direction,
respectively. Using that the columns of $\mv{U}_i$ are orthonormal,
the following expression is derived:
\begin{equation*}
\mathcal{L} \mt{M}  \times_2 \mv{U}_2^\H \times_3 \mv{U}_3^\H = -\mt{G} \times_1 \mv{D}_{xx}
- \mt{G} \times_1 \mv{U}_1 \times_2 \mv{U}_2^\H\mv{D}_{yy} \mv{U}_2
- \mt{G} \times_1 \mv{U}_1 \times_3 \mv{U}_3^\H\mv{D}_{zz} \mv{U}_3
- \left[\mt{K} \circ \left(\mt{G} \times_1 \mv{U}_1 \times_2 \mv{U}_2 \times_3 \mv{U}_3 \right) \right] \times_2 \mv{U}_2^\H \times_3 \mv{U}_3^\H
\end{equation*}
Written in the first unfolding, the multiplication with $\mv{U}_2^\H$ and $\mv{U}_3^\H$ in, respectively, the second and third direction is equivalent to post-multiplication with the matrix
$\left(\mv{I} \otimes \mv{U}_2^\H \right)^\H \left(\mv{U}_3^\H \otimes \mv{I} \right)^\H = \left(\mv{U}_3^\H \otimes \mv{U}_2^\H \right)^\H = \left(\mv{U}_3 \otimes \mv{U}_2 \right)$.

Most of the terms are equal to the case where we had a constant wave
number, see also \eqref{eq:start3d}. Let us focus on the last term
that contains the Hadamard product with the space-dependent wave number, i.e.:
\begin{equation}\label{eq:hadamardProduct}
\begin{split}
\mt{K} &\circ \left(\mt{G} \times_1 \mv{U}_1 \times_2 \mv{U}_2 \times_3 \mv{U}_3 \right).
\end{split}
\end{equation}
For Hadamard products of tensors, $\mt{Z} = \mt{X} \circ \mt{Y}$, the following property for the $k$-th unfolding holds $\mv{Z}_{(k)} = \mv{X}_{(k)} \circ \mv{Y}_{(k)}$.
Thus, written in the first unfolding \eqref{eq:hadamardProduct} is given by
\begin{equation}\label{eq:hadamardProduct1stUnfolding}
\begin{split}
\mv{K}_{(1)} &\circ \mv{M}_{(1)}\\
\mv{K}_{(1)} &\circ \left( \conj{\mv{U}_1} \mv{G}_{(1)} \left( \mv{U}_3 \otimes \mv{U}_2 \right)^\H \right).
\end{split}
\end{equation}

As the Hadamard product-term \eqref{eq:hadamardProduct1stUnfolding} is written in the first unfolding and multiplication with $\mv{U}_2^\H$ and $\mv{U}_3^\H$ in respectively the second and third dimension results in
\begin{equation}\label{eq:lastTerm}
\left[\underbrace{\mv{K}_{(1)}}_{\mv{K}} \circ \underbrace{\conj{\mv{U}_1}\mv{G}_{(1)}}_{\mv{U}} \underbrace{\left( \mv{U}_3 \otimes \mv{U}_2 \right)^\H}_{\mv{V}^\H} \right] \underbrace{\left(\mv{U}_3 \otimes \mv{U}_2 \right)}_{\mv{V}}.
\end{equation}

The derivation of the other terms of \eqref{eq:start3dvar} are equal to the constant wave number case.
The equation in the first unfolding leads to a matrix equation:
\begin{equation}\label{eq:firstUnfoldingVar}
-\conj{\mv{D}_{xx}\mv{U}_1} \mv{G}_{(1)}
- \conj{\mv{U}_1} \mv{G}_{(1)} \left(I \otimes \mv{U}_2^\H \mv{D}_{yy} \mv{U}_2\right)^\H
- \conj{\mv{U}_1} \mv{G}_{(1)} \left(\mv{U}_3^\H \mv{D}_{zz} \mv{U}_3  \otimes I\right)^\H
- \left[\mv{K}_{(1)} \circ \conj{\mv{U}_1}\mv{G}_{(1)} \left( \mv{U}_3 \otimes \mv{U}_2 \right)^\H \right] \left(\mv{U}_3 \otimes \mv{U}_2 \right)
=  \mv{F}_{(1)} \left(\mv{U}_3 \otimes \mv{U}_2\right).
\end{equation}

Vectorization of the last term, i.e. \eqref{eq:lastTerm}, results again in an expression for the space-dependent wave number of the form $\left(\mv{K} \circ \mv{U}\mv{V}^\H\right)\mv{V}$, similar to the two dimensional case which was given in \eqref{eq:2dvarwavenr}.
Using again \eqref{eq:vecThirdTerm}, the vectorization of this expression is given by
\begin{equation}\label{eq:tmp1}
\left(\mv{U}_3^\T \otimes \mv{U}_2^\T \otimes \mv{I} \right) \diag{\mvec{\mv{K}_{(1)}}} \conj{\left(\mv{U}_3 \otimes 	\mv{U}_2 \otimes \mv{I} \right)} \mvec{\conj{\mv{U}_1}\mv{G}_{(1)}}.
\end{equation}

Because $\mt{K}$ is known in a Canonical Polyadic tensor (CP tensor) decomposition\footnote{Otherwise a Canonical Polyadic tensor decomposition can be computed using for example an CP-ALS algorithm [12].}, as given in \eqref{eq:CP4K}, we have
\begin{equation}
\diag{\mvec{\mv{K}_{(1)}}} = \sum_{i=1}^{s} \sigma_i \diag{\mv{v}_{i}^{(3)}} \otimes \diag{\mv{v}_{i}^{(2)}} \otimes \diag{\mv{v}_{i}^{(1)}}.
\end{equation}

So, using the CP tensor representation of the space-dependent wave number the vectorization in \eqref{eq:tmp1} simplifies even further:
\begin{equation*}
\sum_{i=1}^{s} \sigma_i \left(\mv{U}_3^\T \otimes \mv{U}_2^\T \otimes \mv{I} \right) \left[ \diag{\mv{v}_{i}^{(3)}} \otimes \diag{\mv{v}_{i}^{(2)}} \otimes \diag{\mv{v}_{i}^{(1)}} \right] \conj{\left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)} \mvec{\conj{\mv{U}_1}\mv{G}_{(1)}}
\end{equation*}
which reduced to
\begin{equation*}
\underbrace{\sum_{i=1}^{s} \sigma_i
\left(\mv{U}_3^\T \diag{\mv{v}_{i}^{(3)}} \conj{\mv{U}_3} \right) \otimes
\left(\mv{U}_2^\T \diag{\mv{v}_{i}^{(2)}} \conj{\mv{U}_2} \right) \otimes
\left(\diag{\mv{v}_{i_1}^{(1)}}\right)}_{K_1} \mvec{\conj{\mv{U}_1}\mv{G}_{(1)}}.
\end{equation*}
In this way the $K_1$ operator is defined and can be applied to $\mvec{\conj{\mv{U}_1}\mv{G}_{(1)}}$. Observe that this expansion is only advantageous if the space-dependent wave number has low rank, which is typical the case for our applications.

In a similar way, the $K_2$ and $K_3$ operators can be derived:
\begin{equation}\label{eq:Koperators}
\begin{split}
K_1 &= \sum_{i=1}^{s}\sigma_i \left(\mv{U}_3^\T \diag{\mv{v}_{i}^{(3)}} \conj{\mv{U}_3} \right) \otimes \left(\mv{U}_2^\T \diag{\mv{v}_{i}^{(2)}} \conj{\mv{U}_2} \right) \otimes \left(\diag{\mv{v}_{i}^{(1)}}\right) \\
K_2 &= \sum_{i=1}^{s}\sigma_i \left(\mv{U}_3^\T \diag{\mv{v}_{i}^{(3)}} \conj{\mv{U}_3} \right) \otimes \left(\mv{U}_1^\T \diag{\mv{v}_{i}^{(1)}} \conj{\mv{U}_1} \right) \otimes \left(\diag{\mv{v}_{i}^{(2)}}\right) \\
K_3 &= \sum_{i=1}^{s}\sigma_i \left(\mv{U}_2^\T \diag{\mv{v}_{i}^{(2)}} \conj{\mv{U}_2} \right) \otimes \left(\mv{U}_1^\T \diag{\mv{v}_{i}^{(1)}} \conj{\mv{U}_1} \right) \otimes \left(\diag{\mv{v}_{i}^{(3)}}\right).
\end{split}
\end{equation}

So, we find the following linear system to solve for $\mvec{\conj{\mv{U}_1}\mv{G}_{(1)}}$:
\begin{equation}\label{eq:eqnU1var}
\left\{
-\mv{I} \otimes \conj{\mv{D}_{xx}} -
\left[
\conj{\left(\mv{I} \otimes \mv{U}_2^\H \mv{D}_{yy} \mv{U}_2\right)} - \conj{\left(\mv{U}_3^\H \mv{D}_{zz} \mv{U}_3 \otimes \mv{I}\right)}
\right] \otimes \mv{I}
- K_1
\right\}
\mvec{\underbrace{\conj{\mv{U}_1}\mv{G}_{(1)}}_{\mv{X}_1}} = \mvec{\mv{F}_{(1)} \left(\mv{U}_3 \otimes \mv{U}_2\right)}.
\end{equation}
Observe this is a square system with $n_1 \times r_2 r_3$ unknowns
(where the solution in matrix form $\mv{X}_1$ is typical for rank $r >
r_1$).  In a similar way, update equations for $\mv{U}_2$ and
$\mv{U}_3$ are derived by multiplying \eqref{eq:start3dvar} with the
other factor matrices in the appropriate directions:
\begin{equation}\label{eq:eqnU2var}
\left\{
-\mv{I} \otimes \conj{\mv{D}_{yy}} +
\left[ -
\conj{\left(I \otimes \mv{U}_1^\H \mv{D}_{xx} \mv{U}_1\right)} - \conj{\left(\mv{U}_3^\H \mv{D}_{zz} \mv{U}_3 \otimes \mv{I}\right)}
\right]  \otimes \mv{I}
- K_2
\right\}
\mvec{\underbrace{\conj{\mv{U}_2}\mv{G}_{(2)}}_{\mv{X}_2}} = \mvec{\mv{F}_{(2)} \left(\mv{U}_3 \otimes \mv{U}_1\right)}.
\end{equation}
\begin{equation}\label{eq:eqnU3var}
\left\{-\mv{I} \otimes \conj{\mv{D}_{zz}} +
\left[-
\conj{\left(\mv{I} \otimes \mv{U}_1^\H \mv{D}_{xx} \mv{U}_1\right)} - \conj{\left(\mv{U}_2^\H \mv{D}_{yy} \mv{U}_2 \otimes I\right)}
\right] \otimes I
- K_3
\right\} \mvec{\underbrace{\conj{\mv{U}_3}\mv{G}_{(3)}}_{\mv{X}_3}} = \mvec{\mv{F}_{(3)} \left(\mv{U}_2 \otimes \mv{U}_1\right)}.
\end{equation}
Alternating between solving for $\mv{U}_1$, $\mv{U}_2$ and $\mv{U}_3$ using \eqref{eq:eqnU1var}, \eqref{eq:eqnU2var} or \eqref{eq:eqnU3var} results in an algorithm to approximate low-rank tensor solutions for three dimensional problems as given in \eqref{eq:driven3d}. Also in this case the orthogonality of the columns of $\mv{U}_1$, $\mv{U}_2$ and $\mv{U}_3$ are maintained by additional QR factorizations.
So, we derive the algorithm as formulated in Algorithm \ref{alg:var3dv1}. The generalization for dimensions $d > 3$ is straight forward.

\begin{algorithm}[t]
\SetAlgoLined
[$\mt{G}, \mv{U}_1, \mv{U}_2, \mv{U}_3$] = hosvd(initial guess)\;
[$\mv{\Sigma}, \mv{V}_1, \mv{V}_2, \mv{V}_3$] = cp\_als($\mt{K}$)\;
\While{not converged}{
\For{i = 1, 2, 3}{
Compute $K_i$ using \eqref{eq:Koperators}\;
Solve for $\mv{X}_i = \conj{\mv{U}_i}\mv{G}_{(i)} \in \C^{n_i \times r^{d-1}}$ using \eqref{eq:eqnU1var}, \eqref{eq:eqnU2var} or \eqref{eq:eqnU3var}\;
$\conj{\mv{U}_{i}} \mv{G}_{(i)} = \qr{\mv{X}_i(:,~ 1:r_i), 0}$;
$\mt{G} = \texttt{reconstruct} \left[\mv{G}_{(i)}, ~i\right]$\;
$\mt{M} = \mt{G} \times_1 \mv{U}_1  \times_2 \mv{U}_2  \times_3 \mv{U}_3$\;
\caption{Solve for the low-rank tensor decomposition of the solution $\mt{M}$ of a 3D Helmholtz problem with space-dependent wave number (version 1).}
\label{alg:var3dv1}
\end{algorithm}

Similar to the discussion for the constant wave number algorithms,
observe that we solve again for a large matrix $\mv{X}_i \in \C^{n_i \times	r_1r_2r_3/r_i}$.
So, in general the rank of this matrix could be
$\min\left(n_i,~ r_1r_2r_3/r_i\right)$. But it is also known that
$\mv{X}_i = \conj{\mv{U}_i}\mv{G}_{(i)}$ leads to the fact that
the rank of $\mv{X}_i$ should be at most $r_i$.  So selecting the
first $r_i$ columns of $\mv{X}_i$ and computing its QR decomposition
is sufficient to derive a new orthonormal basis as factor matrix

Algorithm \ref{alg:var3dv1} is exactly the space-dependent wave number
equivalent of Algorithm \ref{alg:const3dv1}. The same ideas can be applied
to derive space-dependent wave number alternatives of the algorithms
corresponding to version 2 and 3.
Again, to circumvent solving large systems, we can pre-compute the
QR factorization of $\mv{G}_{(i)}$ and project these equations onto
the obtained $\mv{Q}_i$. Indeed, this will reduce the number of
unknowns in these linear systems to exactly the number of unknowns
as needed for the factor matrices $\conj{\mv{U}}_1$ and

Let us discuss the details. We start again from equation \eqref{eq:firstUnfoldingVar} and use the QR
factorization of $\mv{G}_{(1)}^\H$, $\mv{Q}_1 \mv{R}_1^\H = \qr{\mv{G}^\H_{(1)}}$. This yields
\begin{equation*}
-\conj{\mv{D}_{xx}\mv{U}_1} \mv{R}_1\mv{Q}_1^\H
- \conj{\mv{U}_1} \mv{R}_1\mv{Q}_1^\H \left(I \otimes \mv{U}_2^\H \mv{D}_{yy} \mv{U}_2\right)^\H
- \conj{\mv{U}_1} \mv{R}_1\mv{Q}_1^\H \left(\mv{U}_3^\H \mv{D}_{zz} \mv{U}_3  \otimes I\right)^\H
- \left[\mv{K}_{(1)} \circ \conj{\mv{U}_1}\mv{R}_1\mv{Q}_1^\H \left( \mv{U}_3 \otimes \mv{U}_2 \right)^\H \right] \left(\mv{U}_3 \otimes \mv{U}_2 \right)
=  \mv{F}_{(1)} \left(\mv{U}_3 \otimes \mv{U}_2\right).
\end{equation*}
Post multiplication of the left hand side of this equation by $\mv{Q}_1$ yields
\begin{equation*}
-\conj{\mv{D}_{xx}\mv{U}_1} \mv{R}_1
- \conj{\mv{U}_1} \mv{R}_1\mv{Q}_1^\H \left(I \otimes \mv{U}_2^\H \mv{D}_{yy} \mv{U}_2\right)^\H\mv{Q}_1
- \conj{\mv{U}_1} \mv{R}_1\mv{Q}_1^\H \left(\mv{U}_3^\H \mv{D}_{zz} \mv{U}_3  \otimes I\right)^\H\mv{Q}_1
- \left[\mv{K}_{(1)} \circ \conj{\mv{U}_1}\mv{R}_1\mv{Q}_1^\H \left( \mv{U}_3 \otimes \mv{U}_2 \right)^\H \right] \left(\mv{U}_3 \otimes \mv{U}_2 \right)\mv{Q}_1.
%=  \mv{F}_{(1)} \left(\mv{U}_3 \otimes \mv{U}_2\right)\mv{Q}_1.
\end{equation*}
To solve this equation for $\mv{U}_1$, it is written in vectorized form as
\begin{equation}\label{eq:eqnU1v2forvar}
\left\{
-\mv{I} \otimes \conj{\mv{D}_{xx}} +
\mv{Q}_1^\T\left[-
\conj{\left(\mv{I} \otimes \mv{U}_2^\H \mv{D}_{yy} \mv{U}_2\right)} - \conj{\left(\mv{U}_3^\H \mv{D}_{zz} \mv{U}_3 \otimes \mv{I} \right)}
\right]\conj{\mv{Q}_1} \otimes \mv{I}
-\mv{Q}_1^\T K_1 \conj{\mv{Q}_1}
\right\} \mvec{\underbrace{\conj{\mv{U}_1}\mv{R}_1}_{\mv{X}_1}} = \mvec{\mv{F}_{(1)} \left(\mv{U}_3 \otimes \mv{U}_2\right) \mv{Q}_1}.
\end{equation}
In a similar way, the update equations for $\mv{U}_2$ and $\mv{U}_3$ are
derived by multiplying \eqref{eq:start3dvar} with the other factor
matrices in the appropriate dimensions and using the QR factorizations
of $\mv{G}_{(i)}^\H$:
\begin{equation}\label{eq:eqnU2v2forvar}
\left\{
-\mv{I} \otimes \conj{\mv{D}_{yy}} +
\mv{Q}_2^\T\left[
-\conj{\left(\mv{I} \otimes \mv{U}_1^\H \mv{D}_{xx} \mv{U}_1\right)} - \conj{\left(\mv{U}_3^\H \mv{D}_{zz} \mv{U}_3 \otimes \mv{I} \right)}
\right]\conj{\mv{Q}_2}  \otimes \mv{I}
-\mv{Q}_2^\T K_2 \conj{\mv{Q}_2}
\right\} \mvec{\underbrace{{\conj{\mv{U}_2}\mv{R}_2}}_{\mv{X}_2}}
= \mvec{\mv{F}_{(2)} \left(\mv{U}_3 \otimes \mv{U}_1\right) \mv{Q}_2}.
\end{equation}
\begin{equation}\label{eq:eqnU3v2forvar}
\left\{
-\mv{I} \otimes \conj{\mv{D}_{zz}} +
\mv{Q}_3^\T\left[
-\conj{\left(\mv{I} \otimes \mv{U}_1^\H \mv{D}_{xx} \mv{U}_1\right)} - \conj{\left(\mv{U}_2^\H \mv{D}_{yy} \mv{U}_2 \otimes \mv{I} \right)}
\right]\conj{\mv{Q}_3} \otimes \mv{I}
\right\} \mvec{\underbrace{\conj{\mv{U}_3}\mv{R}_3}_{\mv{X}_3}}
-\mv{Q}_3^\T K_3 \conj{\mv{Q}_3}
= \mvec{\mv{F}_{(3)} \left(\mv{U}_2 \otimes \mv{U}_1\right) \mv{Q}_3}.
\end{equation}

All these equations are cheap to solve. Indeed,
$\mvec{\conj{\mv{U}_i}\mv{R}_i}$ has length $n_i r_i$. Computing a
symmetric reverse Cuthill-McKee permutation of the system matrix one
observes a matrix with a bandwidth $\Oh{r}$, so solving these
equations has a computational cost $\Oh{nr^2}$.

Alternating between solving for $\mv{U}_1$, $\mv{U}_2$ and $\mv{U}_3$ using \eqref{eq:eqnU1v2forvar}, \eqref{eq:eqnU2v2forvar} or \eqref{eq:eqnU3var} results again in an algorithm to approximate low-rank solutions for three dimensional space-dependent Helmholtz problems. Also in this case the orthogonality of the columns of $\mv{U}_1$, $\mv{U}_2$ and $\mv{U}_3$ are maintained by additional QR factorizations.
So, we derive the algorithm as formulated in Algorithm \ref{alg:var3dv3}.
Algorithm \ref{alg:var3dv3} is exactly the space-dependent wave number equivalent of Algorithm \ref{alg:const3dv3}.

\begin{algorithm}[t]
\SetAlgoLined
[$\mt{G}, \mv{U}_1, \mv{U}_2, \mv{U}_3$] = hosvd(initial guess)\;
[$\mv{\Sigma}, \mv{V}_1, \mv{V}_2, \mv{V}_3$] = cp\_als($\mt{K}$)\;
\While{not converged}{
\For{i = 1, 2}{
Compute $K_i$ using \eqref{eq:Koperators}\;
$\mv{Q}_i \widetilde{\mv{R}} = \qr{\mv{G}_{(i)}^\H, 0}$\;
Solve for $\mv{X}_i = \conj{\mv{U}_i}\mv{R}_i \in \C^{n_i \times r_i}$ using \eqref{eq:eqnU1v2forvar} or \eqref{eq:eqnU2v2forvar}\;
$\conj{\mv{U}_{i}} \mv{R}_i = \qr{\mv{X}_i, 0}$\;
$\mt{G} = \texttt{reconstruct} \left[\mv{R}_i\mv{Q}_i^\H, ~i\right]$;
Compute $K_3$ using \eqref{eq:Koperators}\;
Solve for $\mv{X}_3 = \conj{\mv{U}_3}\mv{G}_{(3)} \in \C^{n_3 \times r^{d-1}}$ using \eqref{eq:eqnU3var}\;
$\conj{\mv{U}_{3}} \mv{G}_{(3)} = \qr{\mv{X}_3, 0}$\;
$\mt{G} = \texttt{reconstruct} \left[\mv{G}_{(3)},~3\right]$\;
$\mt{M} = \mt{G} \times_1 \mv{U}_1  \times_2 \mv{U}_2  \times_3 \mv{U}_3$\;
\caption{Solve for the low-rank tensor decomposition of the solution $\mt{M}$ of a 3D Helmholtz problem with space-dependent wave number (version 3).}
\label{alg:var3dv3}
\end{algorithm}

\subsection{Projection operator for space-dependent wave number}
Consider a tensor $\mt{M}$ in Tucker tensor format and factorized as
$\mt{M} = \mt{G} \times_1 \mv{U}_1 \times_2 \mv{U}_2 \times_3 \mv{U}_3$,
with unknowns $\mt{G}, \mv{U}_1, \mv{U}_2$ and $\mv{U}_3$.
Discretization of \eqref{eq:driven3d} with a space-dependent wave number
leads to a linear operator $\mathcal{L}$ applied on tensors. Its matrix
representation $\mv{L}$ has again a structure as given in
\eqref{eq:operatorLasSumKroneckerProducts}.

Solving for the unknown factors $\mv{U}_1$, $\mv{U}_2$ or $\mv{U}_3$
(and the core-tensor $\mt{G}$) using \eqref{eq:eqnU1var},
\eqref{eq:eqnU2var} or \eqref{eq:eqnU3var} can, again, be interpreted as a
projection operator applied on the residual. For example,
\eqref{eq:eqnU1var} can be interpreted as
\begin{equation}
\conj{\left(\mv{U}_3^\H \otimes \mv{U}_2^\H \otimes \mv{I} \right) \mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)} \mvec{\conj{\mv{U}_1}\mv{G}_{(1)}} = \left( \mv{U}_3^\T \otimes \mv{U}_2^\T \otimes \mv{I} \right) \mvec{\mv{F}_{(1)}},
\end{equation}
The residual in tensor format is given by
\begin{equation}
\begin{split}
\mt{R} = \mt{F} &- \mathcal{L}\mt{M}, \\
= \mt{F}
&+ \mt{G} \times_1 \mv{D}_{xx} \mv{U}_1 \times_2 \mv{U}_2 \times_3 \mv{U}_3 \\
&+ \mt{G} \times_1 \mv{U}_1 \times_2 \mv{D}_{yy}\mv{U}_2 \times_3 \mv{U}_3 \\
&+ \mt{G} \times_1 \mv{U}_1 \times_2 \mv{U}_2 \times_3 \mv{D}_{zz}\mv{U}_3 \\
&+ \mt{K} \circ \left(\mt{G} \times_1 \mv{U}_1 \times_2 \mv{U}_2 \times_3 \mv{U}_3 \right).
\end{split}
\end{equation}

Writing this tensor equation in the first unfolding leads to the following matrix equation
\begin{equation}
\begin{split}
\mv{R}_{(1)} = \mv{F}_{(1)}
&+ \conj{\mv{D}_{xx}\mv{U}_1}\mv{G}_{(1)} \left(\mv{U}_3 \otimes \mv{U}_2 \right)^\H \\
&+ \conj{\mv{U}_1}\mv{G}_{(1)} \left(\mv{U}_3 \otimes \mv{D}_{yy}\mv{U}_2 \right)^\H \\
&+ \conj{\mv{U}_1}\mv{G}_{(1)} \left(\mv{D}_{zz}\mv{U}_3 \otimes \mv{U}_2 \right)^\H \\
&+ \mv{K}_{(1)} \circ \left( \conj{\mv{U}_1} \mv{G}_{(1)} \left( \mv{U}_3 \otimes \mv{U}_2 \right)^\H \right)
\end{split}
\end{equation}
which can be matricized as
\begin{equation*}
\begin{split}
\mvec{\mv{R}_{(1)}} &= \mvec{\mv{F}_{(1)}} \\
&- \left(-\conj{\left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{D}_{xx}\right)}
- \conj{\left(\mv{U}_3 \otimes \mv{D}_{yy}\mv{U}_2 \otimes \mv{I}\right)}
- \conj{\left(\mv{D}_{zz}\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)}
- \diag{\mvec{\mv{K}_{(1)}}} \conj{\left(\mv{U}_3 \otimes 	\mv{U}_2 \otimes \mv{I} \right)}\right)\mvec{\conj{\mv{U}_1}\mv{G}_{(1)}}.
\end{split}
\end{equation*}
Rewriting this results in exactly the same structure and projection operator as in the constant wave number case:
\begin{equation}
\begin{split}
\mvec{\mv{R}_{(1)}} &= \ldots \\
&= \mvec{\mv{F}_{(1)}} - \conj{\mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)}\mvec{\conj{\mv{U}_1}\mv{G}_{(1)}} \\
&= \mvec{\mv{F}_{(1)}} - \conj{\mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)} \left[\conj{\left(\mv{U}_3^\H \otimes \mv{U}_2^\H \otimes \mv{I} \right) \mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)}\right]^{-1}\left( \mv{U}_3^\T \otimes \mv{U}_2^\T \otimes \mv{I} \right) \mvec{\mv{F}_{(1)}} \\
&= P_{23}\mvec{\mv{F}_{(1)}}
\end{split}
\end{equation}
where projection operator $P_{23}$ is similar to the projector in the constant wave number case, see \eqref{eq:projector3dconst}, and now given by
\begin{equation}
\begin{split}
P_{23} &= \mv{I} - \conj{\mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)} \left[\conj{\left(\mv{U}_3^\H \otimes \mv{U}_2^\H \otimes \mv{I} \right) \mv{L} \left(\mv{U}_3 \otimes \mv{U}_2 \otimes \mv{I}\right)}\right]^{-1}\left( \mv{U}_3^\T \otimes \mv{U}_2^\T \otimes \mv{I} \right) \\
&= \mv{I} - \mv{X}.
\end{split}
\end{equation}

A similar derivation results in projection operators $P_{13}$ and $P_{12}$ for the updates in $\mv{U}_2$ and $\mv{U}_3$, respectively. Both are also the same as in the constant wave number case, as given in \eqref{eq:projectors3dconst}.

\section{Numerical results} \label{sec:NumericalResults}
In this section, we demonstrate the promising results of the derived
algorithms with some numerical experiments in two and three
dimensions. Furthermore, we consider discretizations of the Helmholtz
equation with constant and space-dependent wave numbers.

\subsection{2D Helmholtz problem with space-dependent wave number}
First, we consider a 2D Helmholtz problem with a space-dependent wave
number given by $k^2(x,y) = 2 + e^{-x^2 -y^2}$.

For this example the two dimensional domain $\Omega = [-10, ~10]^2$ is discretized with $M=1000$ equidistant mesh points per direction in the interior of the domain. Further it is extended with exterior complex scaling to implement the absorbing boundary conditions. In total, the number of discretization points per directions equals $n = n_1 = n_2 = 1668$. As external force $f(x, y) = -e^{-x^2-y^2}$ is applied.

In this space-dependent wave number example it is known that the matrix representation of the semi-exact solution of the Helmholtz equation on the full grid has a low rank. Indeed, approximating the semi-exact solution with a low-rank matrix with rank $r=17$ is in this case sufficient to obtain an error below the threshold $\tau = 10^{-6}$.

Starting with an random (orthonormalized) initial guess for $\mv{V}^{(0)} \in \C^{n \times r}$ only a small number of iterations of Algorithm \ref{alg:2d} is needed to obtain an error similar to the specified threshold $\tau$. As shown in Figure~\ref{fig:error+singularvalues-var2d-r=17} both the residual and the error with respect to the semi-exact solution decay in only a few iterations (i.e. in this example 4-8 iterations) to a level almost similar to the expected tolerance.

The singular values of the approximation $\mv{A}^{(k)} =
\mv{U}^{(k)}{\mv{R}^{(k)}}^\H {\mv{V}^{(k)}}^\H$ in iteration $i$ can
be computed and are shown for increasing iterations in
Figure~\ref{fig:error+singularvalues-var2d-r=17}. As expected the low-
rank approximations recover the singular values of the full grid
semi-exact solution. In fact $\mv{R}^{(k)}$ converges towards

\begin{figure}
\centering
\begin{tabular}{cc}
\input{error-variable-2d-r=17.tikz} &
\input{singularvalues-variable-2d.tikz}
\end{tabular}
\caption{Plot of error and residual (left) and singular values (right) per iteration for space-dependent wave number in 2D Helmholtz problem ($M=1000, r=17$).}
\label{fig:error+singularvalues-var2d-r=17}
\end{figure}

The numerical rank of the matrix representation of the solution of a
Helmholtz problem with a space-dependent wave number is unknown in
advance. But, the presented algorithm is stable with respect to over-
and underestimation of the numerical rank of the solution. Figure
\ref{fig:error+residuals-var2d-M=1000} shows both the error and
residual per iteration and illustrates this statement by approximating
the same semi-exact solution with increasing ranks $r \in \{ 12, 18, 24, 36 \}$.

\begin{figure}
\centering
\begin{tabular}{cc}
\input{error-variable-2d-M=1000.tikz} &
\input{residuals-variable-2d-M=1000.tikz}
\end{tabular}
\caption{Plot of errors (left) and residuals (right) per iteration for space-dependent wave number in 2D Helmholtz problem with increasing ranks ($M=1000$)}
\label{fig:error+residuals-var2d-M=1000}
\end{figure}

\begin{figure}
\centering
\begin{tabular}{cc}
\input{runtime-variable-2d-M=1000.tikz} &
\input{runtime-variable-2d-M=1000+loglog.tikz}
\end{tabular}
\caption{Plot of runtime for 10 iterations for space-dependent wave number in 2D Helmholtz problem with increasing ranks ($M=1000$). Left: linlin-scale, right: loglog-scale.}
\label{fig:runtime-variable-M=1000}
\end{figure}

In contrast to the constant coefficient wave number case the
convergence with space-dependent wave number depends also
on the maximal attainable rank. For increasing maximal attainable
ranks the number of needed iterations decreases. This is especially
observed when the error is considered, but it can also be seen in the
figure where the residuals are shown, Fig.~\ref{fig:error+residuals-var2d-M=1000}.

\subsection{3D Helmholtz problem with space-dependent wave number}
In this example we solve a 3D Helmholtz problem with a space-dependent
wave number discretized on a DVR-grid [20]. All
three versions of the 3D algorithm for space-dependent wave numbers
can successfully be applied.

First, to reduce computational cost of construction of the operators $K_1$, $K_2$ and $K_3$, see \eqref{eq:Koperators}, a CP-decomposition of the space-dependent wave number is constructed.
%Therefore a CP-decomposition using SPLATT is used.
As shown in Figure \ref{fig:cprank_wavenumber} the space-dependent wave number can be well-approximated by a small number of rank-1 tensors. For the examples discussed in this section we used a CP-rank $s = 32$ to approximate this space-dependent wave number. Hence, the error in approximating the wave number is approximately $\Oh{10^{-4}}$.

For all three versions of the algorithm we use 10 iterations of the algorithm to converge to the low-rank solution. For example if we compute the low-rank solution (with $r = r_x = r_y = r_z = 16$) the residual after each iteration for all algorithms is shown in Figure \ref{fig:residual-all-orderDvr=7-solrank=16-waverank=32}.

\begin{figure}
\centering
\begin{subfigure}{.45\textwidth}
\centering
\input{cprank-orderDvr=7-wavenumber.tikz} \\
\caption{CP-rank of space-dependent wave number for 3D Helmholtz problem.}
\label{fig:cprank_wavenumber}
\end{subfigure}
\quad
\begin{subfigure}{.45\textwidth}
\centering
\input{residual-all-orderDvr=7-solrank=16-waverank=32.tikz} \\
\caption{Residual per iteration ($r = 16, s = 32$).}
\label{fig:residual-all-orderDvr=7-solrank=16-waverank=32}
\end{subfigure}
\caption{Low rank approximation to space-dependent wave number and residuals for version 1, version 2 and version 3 of 3D Helmholtz problem with space-dependent wave number (orderDvr = 7).}
\label{fig:residual-all-orderDvr=7}
\end{figure}

\begin{figure}
\centering
\begin{subfigure}{.45\textwidth}
\centering
\input{residual-all-orderDvr=7-waverank=32.tikz} \\
\caption{orderDvr = 7}
\label{fig:residual-all-orderDvr=7-waverank=32}
\end{subfigure}
\quad
\begin{subfigure}{.45\textwidth}
\centering
\input{residual-all-orderDvr=14-waverank=32.tikz} \\
\caption{orderDvr = 14}
\label{fig:residual-all-orderDvr=14-waverank=32}
\end{subfigure}
\caption{Residual after iteration 10 iterations for all three versions of algorithm with $s = 32$.}
\label{fig:residual-all-waverank=32}
\end{figure}

If we increase the maximal attainable rank $r$ of the low-rank approximation, indeed the residual decreases as shown in Figure \ref{fig:residual-all-orderDvr=7-waverank=32}. The residual for version 1 and version 3 are good, while version 2 cannot compete with both other versions by reducing the residual as far as the other versions. Therefore version 1 or version 3, as given in Algorithm~\ref{alg:var3dv1} or Algorithm~\ref{alg:var3dv3} are preferred.

Considering the runtimes of the three versions, similar results as before are observed. In this experiment with \mbox{orderDvr=7} the number of gridpoints equals to $n = 41$. For version 1 and 3, again a runtime of $\Oh{nr^4}$ is observed. The runtime for version 2 splits into two parts: $\Oh{nr^2 + r^9}$. Due to the small rank $r$ and the large number of iterations in these examples algorithm 2 is the fastest version. The runtimes for version 1 and version 3 differ indeed approximately a factor $d$, which makes version 3 better then version 1. The runtimes with orderDvr = 7 (i.e. $n=41$) are shown in Figure \ref{fig:runtimes-all-orderDvr=7-waverank=32+loglog} and with orderDvr = 14 (i.e. $n=90$) are shown in Figure \ref{fig:runtimes-all-orderDvr=14-waverank=32+loglog}.

Comparing the runtimes for orderDvr=7 and orderDvr=14 we see for version 2 (when the rank gets larger) indeed approximately the same runtime independent of orderDvr. Also versions 1 and 3 consume approximately twice as much time which is as expected by the linear dependence on $n$ for both algorithms.

\begin{figure}
\centering
\begin{subfigure}{.45\textwidth}
\centering
\input{runtimes-all-orderDvr=7-waverank=32+loglog.tikz} \\
\caption{orderDvr = 7}
\label{fig:runtimes-all-orderDvr=7-waverank=32+loglog}
\end{subfigure}
\quad
\begin{subfigure}{.45\textwidth}
\centering
\input{runtimes-all-orderDvr=14-waverank=32+loglog.tikz} \\
\caption{orderDvr = 14}
\label{fig:runtimes-all-orderDvr=14-waverank=32+loglog}
\end{subfigure}
\caption{Runtime of 10 iterations for all three versions of algorithm for 3D Helmholtz with space-dependent wave number of low-rank $s = 32$.}
\label{fig:runtimes-all-waverank=32}
\end{figure}

An impression of the low-rank approximation to the wave function is shown in Figures \ref{fig:impression3d} and \ref{fig:impression3dv2}. In this impression the single, double and triple ionization are visible and can be represented by a low-rank wave function.
%For completeness also an impression for the 2D version of this problem is shown in Figure \ref{fig:impression2d}.

\begin{figure}
\centering
\begin{subfigure}{\textwidth}
\centering
\includegraphics[width=0.85\textwidth]{wave2d-LRimpression.png}
\caption{Impression of a low-rank matrix approximation in 2D.}
\label{fig:impression2d}
\end{subfigure}
\par
\begin{subfigure}{\textwidth}
\centering
\includegraphics[width=0.9\textwidth]{wave3d-LRimpression.png}
\caption{Impression of a low-rank tensor approximation in 3D.}
\label{fig:impression3d}
\end{subfigure}
\caption{Impressions of a low-rank approximation of a matrix and a Tucker tensor representing the wave function as solution to a 2D and 3D Helmholtz problem with a space-dependent wave number.}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.9\textwidth]{wave3d+nocube+c1+c2-v2.png}
\caption{Visualization of a 3D wave as low-rank approximation
to a 3D Helmholtz problem with space-dependent wave number
with single, double and triple ionization.}
\label{fig:impression3dv2}
\end{figure}

\section{Discussion and conclusions}\label{sec:discussion}
In this paper we have analyzed the scattering solutions of a driven
Schr\"odinger equations. These describe a  break-up reaction where a
quantum system is fragmented into multiple fragments.  These problems
are equivalent to solving a Helmholtz equations with space-dependent wave

We have shown, first in 2D and then in 3D, that the wave function of
multiple ionization can be well approximated by a low-rank solution.
In 2D, the waves can be represented as a product of two low-rank
matrices and 3D as a low-rank Tucker tensor decomposition.

We propose a method that determines these low-rank components
of the solution directly.  We write the solution as a product of low-rank
components and assume that a guess for all but one component is given.
We then write a linear system for the remaining unknown component.
This is the repeated until each of the components is updated.

This procedure can be interpreted as a series of projections of the
residual on a subspaces and a correction within that subspace.

In theory, the generalization for dimensions $d > 3$ is straightforward. But
for dimensions $d>3$ it starts to be beneficial to change to a Tensor
Train factorization [17]. It is expected that similar
strategy can also be applied to tensors in Tensor Train format.

As demonstrated by the numerical experiments, the presented algorithms
are able to exploit the low-rank structure of the solutions. This
gives the advantage to reduce the number of unknowns and shorten the
computational time to solve the Helmholtz equation.

In two dimensions, the low-rank representation of the solution can be
represented by only $2nr$ unknowns instead of the full grid of
$n^2$ unknowns. Also the linear systems to solve per iteration have
only $nr$ unknowns.

In high-dimensional Helmholtz equations, the low-rank Tucker tensor
decomposition represents the solution with $\Oh{r^d + dnr}$
unknowns. So, the total number of unknowns is reduced, but it is still
exponential in the dimension $d$. For increasing dimensions this leads,
again, to systems with a number of unknowns exponential in $d$. Maybe
other Tucker-like tensor decompositions with a number of unknowns only
polynomial in $d$ can resolve this problem and make the presented
algorithm also applicable for higher dimensions.

\bibliographystyle{plain}
\begin{thebibliography}{10}

[1]
Dominique Akoury, K~Kreidi, T~Jahnke, Th~Weber, A~Staudte, M~Schoffler,
N~Neumann, J~Titze, L~Ph~H Schmidt, A~Czasch, et~al.
\newblock The simplest double slit: interference and entanglement in double
photoionization of {H2}.
\newblock {\em Science}, 318(5852):949--952, 2007.

[2]
Jean-Pierre Berenger.
\newblock A perfectly matched layer for the absorption of electromagnetic
\newblock {\em Journal of computational physics}, 114(2):185--200, 1994.

[3]
JS~Briggs and V~Schmidt.
\newblock Differential cross sections for photo-double-ionization of the helium
\newblock {\em Journal of Physics B: Atomic, Molecular and Optical Physics},
33(1):R1, 2000.

[4]
Weng~Cho Chew and William~H Weedon.
\newblock A {3D} perfectly matched medium from modified {M}axwell's equations
with stretched coordinates.
\newblock {\em Microwave and optical technology letters}, 7(13):599--604, 1994.

[5]
Siegfried Cools and Wim Vanroose.
\newblock A fast and robust computational method for the ionization cross
sections of the driven {S}chr{\"o}dinger equation using an $\mathcal{O}(n)$
multigrid-based scheme.
\newblock {\em Journal of {C}omputational {P}hysics}, 308:20--39, 2016.

[6]
Sergey~V Dolgov and Dmitry~V Savostyanov.
\newblock Alternating minimal energy methods for linear systems in higher
\newblock {\em SIAM Journal on Scientific Computing}, 36(5):A2248--A2271, 2014.

[7]
Andr{\'e} Gaul and Nico Schl{\"o}mer.
\newblock Preconditioned recycling {K}rylov subspace methods for self-adjoint
\newblock {\em arXiv preprint arXiv:1208.0264}, 2012.

[8]
Lars Grasedyck, Daniel Kressner, and Christine Tobler.
\newblock A literature survey of low-rank tensor approximation techniques.
\newblock {\em GAMM-Mitteilungen}, 36(1):53--78, 2013.

[9]
Wolfgang Hackbusch.
\newblock {\em Tensor spaces and numerical tensor calculus}, volume~42.
\newblock Springer, 2012.

[10]
Sebastian Holtz, Thorsten Rohwedder, and Reinhold Schneider.
\newblock The alternating linear scheme for tensor optimization in the tensor
train format.
\newblock {\em SIAM Journal on Scientific Computing}, 34(2):A683--A713, 2012.

[11]
Thomas Huckle, Konrad Waldherr, and Thomas Schulte-Herbr{\"u}ggen.
\newblock Computations in quantum tensor networks.
\newblock {\em Linear Algebra and its Applications}, 438(2):750--781, 2013.

[12]
Tamara~G Kolda and Brett~W Bader.
\newblock Tensor decompositions and applications.
\newblock {\em SIAM review}, 51(3):455--500, 2009.

[13]
C~William McCurdy, M~Baertschy, and TN~Rescigno.
\newblock Solving the three-body coulomb breakup problem using exterior complex
\newblock {\em Journal of Physics B: Atomic, Molecular and Optical Physics},
37(17):R137, 2004.

[14]
Valentin Murg, Frank Verstraete, {\"O}rs Legeza, and Reinhard~M Noack.
\newblock Simulating strongly correlated quantum systems with tree tensor
\newblock {\em Physical Review B}, 82(20):205105, 2010.

[15]
\newblock Functional operators vol. ii. the geometry of orthogonal spaces.
annals of mathematical studies\# 22.
\newblock 1950.

[16]
Ivan Oseledets.
\newblock {DMRG} approach to fast linear algebra in the {TT}-format.
\newblock {\em Computational Methods in Applied Mathematics}, 11(3):382--393,

[17]
Ivan~V Oseledets.
\newblock Tensor-train decomposition.
\newblock {\em SIAM Journal on Scientific Computing}, 33(5):2295--2317, 2011.

[18]
\newblock The exact solution for a simplified model of electron scattering by
hydrogen atoms.
\newblock {\em Journal of {P}hysics {B}: {A}tomic and {M}olecular {P}hysics},
11(17):3081, 1978.

[19]
TN~Rescigno, M~Baertschy, WA~Isaacs, and CW~McCurdy.
\newblock Collisional breakup in a quantum system of three charged particles.
\newblock {\em Science}, 286(5449):2474--2479, 1999.

[20]
TN~Rescigno and CW~McCurdy.
\newblock Numerical grid methods for quantum-mechanical scattering problems.
\newblock {\em Physical review A}, 62(3):032706, 2000.

[21]
H~Schmidt-B{\"o}cking, R~D{\"o}rner, and J~Ullrich.
\newblock Coltrims.
\newblock {\em {E}urophysics {N}ews}, 33(6):210--211, 2002.

[22]
Barry Simon.
\newblock The definition of molecular resonance curves by the method of
exterior complex scaling.
\newblock {\em Physics Letters A}, 71(2-3):211--214, 1979.

[23]
Zachary~L Streeter, Frank~L Yip, Robert~R Lucchese, Benoit Gervais, Thomas~N
Rescigno, and C~William McCurdy.
\newblock Dissociation dynamics of the water dication following one-photon
double ionization. i. {T}heory.
\newblock {\em Physical Review A}, 98(5):053429, 2018.

[24]
Aaron Temkin.
\newblock Nonadiabatic theory of electron-hydrogen scattering.
\newblock {\em Physical Review}, 126(1):130, 1962.

[25]
Ledyard~R Tucker.
\newblock Some mathematical notes on three-mode factor analysis.
\newblock {\em Psychometrika}, 31(3):279--311, 1966.

[26]
Wim Vanroose, Daniel~A Horner, Fernando Martin, Thomas~N Rescigno, and
C~William McCurdy.
\newblock Double photoionization of aligned molecular hydrogen.
\newblock {\em Physical Review A}, 74(5):052702, 2006.

[27]
Wim Vanroose, Fernando Martin, Thomas~N Rescigno, and C~William McCurdy.
\newblock Complete photo-induced breakup of the {H2} molecule as a probe of
molecular electron correlation.
\newblock {\em Science}, 310(5755):1787--1789, 2005.

[28]
Jinchao Xu and Ludmil Zikatanov.
\newblock The method of alternating projections and the method of subspace
corrections in {H}ilbert space.
\newblock {\em Journal of the American Mathematical Society}, 15(3):573--597,

\end{thebibliography}

\end{document}