Density Matrix Reconstructions in Ultrafast Transmission Electron Microscopy: Uniqueness, Stability, and Convergence Rates

TL;DR

This paper provides a mathematical analysis of the inverse problem in ultrafast transmission electron microscopy, focusing on the uniqueness, stability, and convergence rates of density matrix reconstructions, complementing prior experimental work.

Contribution

It offers a theoretical framework analyzing the inverse problem's properties, including conditions for uniqueness and stability, and convergence rates, enhancing understanding of the reconstruction process.

Findings

Analysis of conditions for unique density matrix reconstruction

Stability estimates under various a-priori information

Convergence rates of reconstruction algorithms

Abstract

In the recent paper [17] the first experimental determination of the density matrix of a free electron beam has been reported. The employed method leads to a linear inverse problem with a positive semidefinite operator as unknown. The purpose of this paper is to complement the experimental and algorithmic results in the work mentioned above by a mathematical analysis of the inverse problem concerning uniqueness, stability, and rates of convergence under different types of a-priori information.

Equations

$$(U_{\omega }(\theta ){)}_{k,l}=e^{\mathrm{i}(k-l)\theta }{J}_{k-l}(2|g_{\omega }|),k,l\in \text{mathds}Z.$$

$${\rho }_{\text{out}}(\theta )=U_{\omega }(\theta )\rho U_{\omega }^{\ast }(\theta ),$$

$$T\rho =y$$

$$(T\rho )(l,\theta )=p(l,\theta )=(U_{\omega }(\theta )\rho U_{\omega }^{\ast }(\theta ){)}_{l,l}.$$

$${\hat{\rho }}_{\alpha }=\arg \underset{\rho }{\min }\Vert T\rho -y^{\mathrm{obs}}{\Vert }_{\mathbb{Y}}^2+\alpha \Vert \rho {\Vert }_{\mathbb{X}}^2\text{mbox}subjectto\rho \ge 0,\mathrm{trace}(\rho )=1$$

$${\Vert T{\rho }^{†}-y^{\mathrm{obs}}\Vert }_{\mathbb{Y}}<\delta$$

$$\delta \le {\Vert T{\hat{\rho }}_{\alpha }-y^{\mathrm{obs}}\Vert }_{\mathbb{Y}}\le \tau \delta$$

$${\rho }_{n+k,n}=0,\text{mbox}forall|k|>k_0$$

$$\Vert {\rho }^{(1)}-{\rho }^{(2)}\Vert \le C{\Vert T{\rho }^{(1)}-T{\rho }^{(2)}\Vert }^{\frac{1}{1+2k_0}}$$

$$\Vert {\hat{\rho }}_{\alpha }-{\rho }^{†}\Vert \le C{\delta }^{\frac{1}{1+2k_0}}$$

$$\displaystyle\sum_{n=-\infty}^{\infty}|\rho_{n+k,n}|\leq C_{\rho}|k|^{-\frac{1}{2}-2\mu}\quad\mbox{for some $\mu>0$ or}$$

$$\sum _{n=-\infty }^{\infty }|{\rho }_{n+k,n}|\le C_{\rho }{b}^{|k|}\text{mbox}forsomeb<1$$

$$\Vert {\rho }^{(1)}-{\rho }^{(2)}\Vert \le \Phi \left(\Vert T{\rho }^{(1)}-T{\rho }^{(2)}\Vert \right)$$

$$\Vert {\hat{\rho }}_{\alpha }-{\rho }^{†}\Vert \le 2(1+\tau )\Phi (\delta )$$

$$\Phi (\delta ):=\{\begin{array}{ll}C{\left(\frac{-\log \delta }{\log (-\log \delta )}\right)}^{-2\mu }& \text{mbox}incaseof\text{eqref}Log-Assum,\\ C\exp \left(-\sqrt{(-\log \delta )(-\log b)}\right)& \text{mbox}incaseof\text{eqref}Hol-Assum2,\end{array}$$

$$(\mathcal{F}f{)}_n$$

$$({\mathcal{F}}^{-1}a)(\theta )$$

$$2\pi \sum _{n\in \text{mathds}Z}{e}^{i\phi n}{a}_n{b}_n={\int }_{-\pi }^{\pi }({\mathcal{F}}^{\ast }a)(\tau )({\mathcal{F}}^{\ast }b)(\phi -\tau )d\tau ,\phi \in [-\pi ,\pi ]$$

$$\begin{array}{rl}& {\mathcal{F}}^{(\ast 1)}:{\ell }^2(\text{mathds}{Z}^2)\to L^2([-\pi ,\pi ]\times \text{mathds}Z)\\ & ({\mathcal{F}}^{\ast (1)}a)(\theta ,m):=\frac{1}{\sqrt{2\pi }}\sum _{n\in \text{mathds}Z}{e}^{\mathrm{i}n\theta }{a}_{n,m},\end{array}$$

$$\begin{array}{rl}& {\mathcal{F}}^{(1,\ast 2)}:L^2([-\pi ,\pi ]\times \text{mathds}Z)\to L^2(\text{mathds}Z\times [-\pi ,\pi ])\\ & ({\mathcal{F}}^{(1,\ast 2)}f)(n,\phi ):=\frac{1}{2\pi }\sum _{n\in \text{mathds}Z}{\int }_{-\pi }^{\pi }{e}^{\mathrm{i}n\theta }{e}^{\mathrm{i}\phi m}f(\theta ,m)d\theta .\end{array}$$

$$({\mathcal{F}}^{(1,\ast 2)}{)}^{-1}=({\mathcal{F}}^{(1,\ast 2)}{)}^{\ast }={\mathcal{F}}^{(\ast 1,2)}.$$

$$T={\mathcal{F}}^{(\ast 1,2)}{M}_m{\mathcal{F}}^{(1)}\mathcal{G}.$$

$$\begin{array}{rl}& M_{m}f(\phi ,k):=f(\phi ,k)m(\phi ,k)\text{mbox}with\\ & m(\phi ,k):=\sqrt{2\pi }{\mathrm{i}}^k{e}^{\frac{\mathrm{i}\phi k}{2}}{J}_k\left(4|g_{\omega }|\sin \frac{\phi }{2}\right),\end{array}$$

$$(\mathcal{G}\rho {)}_{n,k}:={\rho }_{n+k,n},$$

$$\begin{array}{rl}(T\rho )(l,\theta )& =(U_{\omega }(\theta )\rho U_{\omega }^{\ast }(\theta ){)}_{l,l}\\ & =\sum _{m,n}(U_{\omega }(\theta ){)}_{l,m}{\rho }_{m,n}(\overline{U_{\omega }^{\ast }(\theta )}{)}_{n,l}\\ & =\sum _{m,n}{e}^{\mathrm{i}(l-m)\theta }{J}_{l-m}(2|g_{\omega }|){\rho }_{m,n}{e}^{\mathrm{i}(n-l)\theta }\overline{J_{l-n}(2|g_{\omega }|)}.\end{array}$$

$$\left({\mathcal{F}}^{(1,\ast 2)}T\rho \right)(\phi ,k)=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }\sum _{l\in \text{mathds}Z}{e}^{\mathrm{i}\phi l}{e}^{\mathrm{i}\theta k}(T\rho )(l,\theta )d\theta$$

$$\left({\mathcal{F}}^{(1,\ast 2)}T\rho \right)(\phi ,k)$$

$$=\sum _{n,l\in \text{mathds}Z}{e}^{\mathrm{i}\phi l}{J}_{l-n-k}(2|g_{\omega }|)\overline{J_{l-n}(2|g_{\omega }|)}{\rho }_{n+k,n}$$

$$=\left(\sum _{n^{\prime }\in \text{mathds}Z}{e}^{\mathrm{i}\phi n^{\prime }}{J}_{n^{\prime }-k}(2|g_{\omega }|)\overline{J_{n^{\prime }}(2|g_{\omega }|)}\right)\left(\sum _{n\in \text{mathds}Z}{e}^{\mathrm{i}\phi n}{\rho }_{n+k,n}\right)$$

$$\sum _{n\in \text{mathds}Z}{e}^{\mathrm{i}\tau n}{\rho }_{n+k,n}=\sum _{n\in \text{mathds}Z}{e}^{\mathrm{i}\tau n}(\mathcal{G}\rho {)}_{n,k}=\sqrt{2\pi }({\mathcal{F}}^{(1)}(\mathcal{G}\rho ))(\phi ,k)$$

$$\tilde{m}(\phi ,k):=\sqrt{2\pi }\sum _{n^{\prime }\in \text{mathds}Z}{e}^{\mathrm{i}\phi n^{\prime }}{J}_{n^{\prime }-k}(2|g_{\omega }|)\overline{J_{n^{\prime }}(2|g_{\omega }|)}$$

$$\left({\mathcal{F}}^{(1,\ast 2)}T\rho \right)(\phi ,k)=\tilde{m}(\phi ,k)\left({\mathcal{F}}^{(1)}(\mathcal{G}\rho )\right)(\phi ,k).$$

Full text

Keywords: Tikhonov regularization, variational source conditions, uniqueness, stability, electron microscopy, SQUIRRELS

1 Introduction

The density matrix is a fundamental notion in quantum statistics which describes the statistical state of an ensemble of identical single or many body quantum systems. It is a positive semidefinite operator of trace $1$ on the Hilbert space describing the state of a single quantum system. In the area of quantum optics, there are well-established techniques for characterizing the quantum state of the electromagnetic field in terms of its density matrix [14, 19]. Such ‘quantum state tomography’ facilitates the discrimination of, for example, coherent states, squeezed states, thermal states or photon number (Fock) states. In contrast, the reconstruction of the quantum state of a beam of free electrons has only recently been established using inelastic electron-light scattering [17]. The reconstruction technique, termed ‘SQUIRRELS’ for ‘Spectral Quantum Interference for the Regularized Reconstruction of free ELectron States’), is experimentally based on the principle of ‘Photon-Induced Near Field Electron Microscopy’ (PINEM) [2]. In PINEM, a beam of electrons is passed through the near field of laser-illuminated nanostructures or thin films, leading to the formation of sidebands in the electron energy spectrum, spaced by the photon energy [15, 6]. The spatially varying number of created sidebands yields the optical field strength with very high resolution on the nanometer scale [2, 16]. However, the quantum coherent nature of the electron-light interaction has also led to the observation of other fundamental quantum effects, such as multilevel Rabi oscillations [4, 7] or Ramsey-type phase interference in spatially separated fields [3]. Moreover, it has recently been shown that the interaction can be used to temporally structure electron beams into a train of attosecond pulses, the duration of which was determined by SQUIRRELS [17]. The various existing and future applications of inelastic electron-light scattering and the relevance of the specific electron state resulting from such interactions calls for a solid mathematical basis underlying the quantum state reconstruction scheme. The mathematical aspects of SQUIRRELS involve a linear inverse problem with the density matrix as unknown. In [17] this inverse problem was solved by Tikhonov regularization with positive semidefiniteness and trace constraints using quadratic semi-definite programming. The purpose of the present paper is to provide mathematical foundations of the SQUIRRELS method.

In the experiment in [17], light reflection from a thin graphite film mediates the interaction of free electrons with laser photons of two frequencies $\omega$ and $2\omega$ and a controllable relative phase $\theta$. As the interactions of electrons with laser photons lead to a comb-type energy spectrum of the electrons separated by the photon energy, the Hilbert space describing the state of an electron may be chosen as $l^{2}(\mathds{Z})$. The effect of the interaction of an $\omega$-photon with a single electron is described by a unitary operator $U_{\omega}(\theta):l^{2}(\mathds{Z})\to l^{2}(\mathds{Z})$ given in matrix representation by

[TABLE]

Here $J_{l}(2|g_{\omega}|)$ denotes the Bessel function of the first kind of order $l$, and $g_{\omega}$ is a coupling constant associated with the laser. The effect of the photon-electron interaction on the free-electron density matrix $\rho=(\rho_{j,k})_{j,k\in\mathds{Z}}$ is then described by

[TABLE]

where $\rho$ and $\rho_{\text{out}}(\theta)$ are the density matrices before and after the interaction, respectively. However, only the diagonal values $\rho_{\text{out}}(\theta)$ are observable. On the other hand, since the phase parameter $\theta\in[-\pi,\pi]$ is experimentally controllable, we may observe a spectrogram $y(\theta,l)=\left(\rho_{\text{out}}(\theta)_{ll}\right)_{l\in\mathds{Z}}$ for each value of $\theta$. The inverse problem to find the electron density matrix $\rho$ from the measured data $y$ is then described by the operator equation

[TABLE]

with a bounded linear forward operator $T:\mathbb{X}\to\mathbb{Y}$ between Hilbert spaces $\mathbb{X}:=l^{2}(\mathds{Z}\times\mathds{Z})$ and $\mathbb{Y}:=L^{2}(\mathds{Z}\times[-\pi,\pi])$ given by

[TABLE]

The aim of this paper is to analyze this inverse problem mathematically concerning uniqueness, stability, and rates of convergence.

The plan of the remainder of this paper is organized as follows: In section 2 we formulate our main results. Section 3 is devoted to the proof that $T$ is injective. It is based on a factorization of $T$ which is also fundamental for the rest of this paper. Our main tool for the proofs of stability estimates and convergence rates are variational source conditions, which will be treated in section 4.

2 Main results

Our first main result asserts that the unknown density matrix is in fact uniquely determined by the experimental data in the absence of noise and modelling errors:

Theorem 1** (uniqueness).**

The operator $T$ defined in (3) is injective.

It follows from our analysis (see Corollary 3.2) and has been observed in numerical experiments that the inverse problem (2) is ill-posed. Therefore, a natural question concerns the degree of ill-posedness or the degree of stability that can be obtained under certain types of a-priori information on the true density matrix $\rho^{\dagger}$. This will be addressed in the following three theorems.

Moreover, to obtain stable reconstruction for noisy experimental data $y^{\mathrm{obs}}$, some kind of regularization has to be employed. In [17] constrained Tikhonov regularization of the following form has been used:

[TABLE]

Here we minimize only over density matrices, i.e. positive semidefinte operators of trace $1$. The regularization parameter $\alpha>0$ is chosen by the discrepancy principle as follows: Let $\delta>0$ be the deterministic noise level, i.e.

[TABLE]

for the true density matrix $\rho^{\dagger}$. Then $\alpha>0$ is chosen such that

[TABLE]

for some $\tau>1$. We will also derive error bounds for Tikhonov regularization described by (4) and (5).

We first consider band-limited density matrices:

Theorem 2** (Hölder-type estimates for band limited $\rho$).**

Suppose the density matrices $\rho^{(1)}$, $\rho^{(2)}$, and $\rho^{\dagger}$ are band-limited, i.e.

[TABLE]

for $\rho\in\{\rho^{\dagger},\rho^{(1)},\rho^{(2)}\}$ and some $k_{0}>0$. Then the stability estimate

[TABLE]

holds true for some constant $C$ depending only on $k_{0}$ and $|g_{\omega}|$. Moreover, for the Tikhonov regularized solution $\widehat{\rho}_{\alpha}$ given by (4) with parameter choice rule $\alpha\sim\delta^{\frac{2+2k_{0}}{1+2k_{0}}}$, or with $\alpha$ chosen by the discrepancy principle (5), the error bound

[TABLE]

is satisfied for all $\delta\in(0,1]$ with $C$ independent of $\rho^{\dagger}$.

If we relax band-limitation to a polynomial or exponential decay condition, we only obtain slower than Hölder stability estimates and convergence rates:

Theorem 3** (Sub-Hölder rates under decay conditions).**

Suppose there exists $C_{\rho}>0$ such that the off-diagonal entries of the density matrices $\rho^{(1)}$, $\rho^{(2)}$, and $\rho^{\dagger}$ satisfy either a exponential or a polynomial decay condition:

[TABLE]

and for all $k\in\mathbb{Z}$ and $\rho\in\{\rho^{(1)},\rho^{(2)},\rho^{\dagger}\}$. Then the stability estimate

[TABLE]

holds true for $\left\|T\rho^{(1)}-T\rho^{(2)}\right\|\leq\frac{1}{2}$, and for the Tikhonov regularized solution $\widehat{\rho}_{\alpha}$ given by (4) with $\alpha$ chosen by the discrepancy principle (5), the error is bounded by

[TABLE]

where the function $\Phi$ is given by

[TABLE]

for $\delta\in(0,1/2]$ with some constant $C$ depending only on $C_{\rho}$ and $|g_{\omega}|$, $\mu$, and $b$.

Note that the logarithmic rate $\left(\frac{-\log\delta}{\log(-\log\delta)}\right)^{-\mu}$ for the polynomial decay condition (8a) is slower than $\mathcal{O}((-\log\delta)^{-\mu})$ as $\delta\to 0$, but faster than $\mathcal{O}((-\log\delta)^{-\mu^{\prime}})$ for any $\mu^{\prime}<\mu$. On the other hand, the rate for the exponential decay condition (8b) is slower that any Hölder rate $\mathcal{O}(\delta^{\nu})$ for $\nu>0$, but faster that any logarithmic rate. Such rates of convergence and stability estimates occur much less frequently than Hölder and logarithmic rates, but similar rates have been derived e.g. in scattering theory for the reconstruction of a near field data from far field data (see [9, Lemma 4.2]).

3 Uniqueness

The aim of this section is to prove Theorem 1. Throughout this paper we use the following notations:

Let $\mathcal{F}:L^{2}([-\pi,\pi])\to\ell^{2}(\mathds{Z})$,

[TABLE]

denote the periodic Fourier transform. Since the scaling factor is chosen such that $\mathcal{F}$ is unitary, we have

[TABLE]

Recall the periodic Fourier convolution theorem

[TABLE]

for $a,b\in\ell^{2}(\mathds{Z})$. Fourier transforms on spaces of multi-variate functions will be labeled by superscript(s) indicating the position of the variable(s) on which they act. For example,

[TABLE]

It is easy to see that these operators are again unitary. (This can either be proved directly or by noting that they are tensor products of unitary operator, $\mathcal{F}^{*1}=\mathcal{F}^{*}\otimes I$ and $\mathcal{F}^{(1,*2)}=\mathcal{F}\otimes\mathcal{F}^{*}$.) In particular,

[TABLE]

Our basic tool for the uniqueness proof and also for the following sections is the following factorization of the operator $T$:

Proposition 3.1**.**

The operator $T$ defined in (3) has the factorization

[TABLE]

Here $\mathcal{F}^{(1)}$ and $\mathcal{F}^{(*1,2)}$ are defined in (LABEL:FourierTransformF1)–(12), $M_{m}:L^{2}([-\pi,\pi]\times\mathds{Z})\to L^{2}([-\pi,\pi]\times\mathds{Z})$ is the multiplication operator

[TABLE]

and $\mathcal{G}:\ell^{2}(\mathds{Z}^{2})\to\ell^{2}(\mathds{Z}^{2})$ is a matrix shift operator defined by

[TABLE]

i.e. the $k$-th column of $\mathcal{G}\rho$ is the $k$-th diagonal of $\rho$.

Proof.

By plugging the definition of $U_{\omega}$ into the definition of $T$ and using $U^{*}_{n,l}=\overline{U_{l,n}}$, we obtain

[TABLE]

Applying $\mathcal{F}^{(1,*2)}$ to both sides of this equation yields

[TABLE]

As $\int_{-\pi}^{\pi}e^{\mathrm{i}\theta(n-m+k)}d\theta=2\pi\delta_{m-n.k}$, this simplifies to

[TABLE]

where we have used the substitution $l=n^{\prime}+n$ in the last line. Using the identity

[TABLE]

and setting

[TABLE]

leads to the formula

[TABLE]

In view of (12) it remains to show that $\tilde{m}=m$. To this end we apply the Fourier convolution theorem (9) to obtain

[TABLE]

With the help of the identity $\sum_{m}e^{\mathrm{i}m\theta}J_{m}(z)=e^{\mathrm{i}z\sin\theta}$ for Bessel functions (see [20, eq. (9.20)]), we obtain

[TABLE]

Now the identity

[TABLE]

(see [20, eq. (9.19)]) shows that $\widetilde{m}=m$ and completes the proof in view of (16). ∎

The factorization (13) leads us to a proof of the injectivity of $T$.

Proof of Theorem 1:.

Since the Fourier-type transforms ${\mathcal{F}^{(1)}}$ and ${\mathcal{F}^{(*1,2)}}$ and the operator $\mathcal{G}$ in the factorization (13) are all unitary, it suffices to prove the multiplication operator $M$ is injective. To this end, we notice that the multiplier $\sqrt{2\pi}\mathrm{i}^{k}e^{\frac{\mathrm{i}\varphi k}{2}}J_{k}(4|g_{\omega}|\sin\frac{\varphi}{2})$ is a holomorphic function with respect to $\varphi$, so it only has isolated zeros. Therefore, if $Mf=0$ for some $f\in L^{2}([-\pi,\pi]\times\mathds{Z})$, we are able to infer that $f$ vanishes almost everywhere. This means $f=0$ in $L^{2}([-\pi,\pi]\times\mathds{Z})$, which concludes the proof. ∎

Corollary 3.2** (ill-posedness).**

The inverse of $T$ is unbounded.

Proof.

It follows from the factorization (13) that

[TABLE]

All operators on the right hand side are unitary except the multiplication operator $M_{m}^{-1}=M_{1/m}$. A multiplication operator on $L^{2}$ is bounded if and only if the multiplier function is essentially bounded. Since the Bessel function $J_{k}$ has zeroes at [math] for $k\geq 1$, $1/m$ is not essentially bounded. ∎

4 Variational source conditions

In this section we review variational source conditions and verify conditions of this type, which then imply both the stability estimates and the convergence rates in Theorems 2 and 3.

4.1 Basic theory

We first recall some standard regularization theory for inverse problems. We consider a general linear ill-posed inverse problem $Tf^{\dagger}=y^{\mathrm{obs}}$ where $T:\mathbb{X}\to\mathbb{Y}$ is a bounded linear operator between Hilbert spaces, which does not have a bounded inverse. Let $y^{\mathrm{obs}}$ is the noisy data with noise level $\|y^{\mathrm{obs}}-Tf^{\dagger}\|_{\mathbb{Y}}<\delta$. Tikhonov regularization with constraint set $\mathcal{C}\subset\mathbb{X}$ and regularization parameter $\alpha>0$ is given by

[TABLE]

To obtain bounds on the reconstruction error $\|\widehat{f}_{\alpha}-f^{\dagger}\|$, which tend to [math] as the noise level tends to [math], we need to impose conditions on $f^{\dagger}$. Such conditions are usually referred to as source condition. Classically, the source condition is of the form

[TABLE]

where $\omega\in\mathbb{X}$ is some “source” and $h:[0,\|T^{*}T\|]\to[0,\infty)$ is some increasing function, which determines the convergence rate. The usefulness of such spectral source conditions is linked to the applicability of tools from spectral theory which is mostly restricted to linear reconstruction procedure. Even for linear reconstruction procedures such as unconstrained Tikhonov regularization they are only sufficient, but not quite necessary for certain convergence rates. Both of these shortcomings can be overcome by the use of source conditions in the form of variational inequalities (see [10, 18, 5, 21]):

Definition 4.1** (Variational Source Condition).**

A function $\psi:[0,\infty)\to[0,\infty)$ is called an index function if it is continuous, strictly increasing, and $\psi(0)=0$. We say that $f^{\dagger}$ satisfies a variational source condition with index function $\psi$ if

[TABLE]

Conditions of the form (18) lead to the following error bounds in terms of the index function $\psi$ (see [8, 5]):

Proposition 4.2** (convergence rates with a-priori choice of $\alpha$).**

Consider Tikhonov regularization given by eq. (17). If $f^{\dagger}\in\mathcal{C}$ satisfies a variational source condition (18) with some concave, differentiable index function $\psi$ and if the regularization parameter $\alpha$ is chosen by $\alpha=1/\psi^{\prime}(4\delta^{2})$, then the following error bound holds true:

[TABLE]

Actually eq. (18) only needs to hold for all $f\in\mathcal{C}$. The same convergence rate can also be achieved by the discrepancy principle, which does not require prior knowledge of the index function $\psi$ encoding properties of the unknown solution.

Proposition 4.3** (convergence rates with discrepancy principle).**

Consider Tikhonov regularization given by eq. (17). If $f^{\dagger}\in\mathcal{C}$ satisfies a variational source condition (18) with some concave, differentiable index function $\psi$ and if the regularization parameter $\alpha$ is chosen according to the discrepancy principle (5), then the following error bound holds true:

[TABLE]

Proof.

The proof is adapted from [11, Theorem 4.3(iii)], see also [5]. We first notice that since $\widehat{f}_{\alpha}$ is defined to be the minimizer of the Tikhonov functional (17), it follows from (5) that

[TABLE]

Combining this inequality with (18) yields

[TABLE]

As

[TABLE]

and $\psi$ is monotonically increasing, we obtain $\|\widehat{f}_{\alpha}-f^{\dagger}\|^{2}\leq\psi\left((1+\tau)^{2}\delta^{2}\right)$. Now the proof is completed by noting that $\psi\left((1+\tau)^{2}\delta^{2}\right)\leq(1+\tau)^{2}\psi\left(\delta^{2}\right)$ as $\psi$ is concave with $\psi(0)=0$. ∎

Here we do not address the question whether a parameter $\alpha>0$ satisfying (5) exists. We refer to [1] for the so-called sequential discrepancy principle, which determines $\alpha$ by an explicit algorithm for which similar error bounds can be shown.

Remark 4.4** ([12, Eq. 6]).**

If the variational source condition (18) is satisfied for all $f^{\dagger}$ in some subset $\mathcal{K}\subset\mathbb{X}$, then the conditional stability estimate

[TABLE]

holds true for all $f_{1},f_{2}\in\mathcal{K}$.

To verify variational source conditions for our problem, we will check the sufficient conditions in the following lemma, which is a special case of [13, Theorem 2.1] (where a different scaling is used such that the index function $\psi$ in [13] is twice the Bessel function $\psi$ here):

Lemma 4.5** (Verification of VSCs).**

Let $\mathbb{X}$ and $\mathbb{Y}$ be Hilbert spaces and $V_{\varepsilon}\in\mathbb{X}$ be a family of subspaces. Suppose that there exists a family of orthogonal projection operators $P_{\varepsilon}:\mathbb{X}\rightarrow V_{\varepsilon}$ for $\varepsilon$ in some index set $\mathcal{I}\subset(0,\infty)$, and there exist families $(\kappa_{\varepsilon})_{\varepsilon\in\mathcal{I}}$, $(\sigma_{\varepsilon})_{\varepsilon\in\mathcal{I}}$ of positive numbers, such that the following conditions hold true:

• •

$\|(I-P_{\varepsilon})f^{\dagger}\|_{\mathbb{X}}\leq\kappa_{\varepsilon}$* for all $\varepsilon\in\mathcal{I}$;*

• •

$\inf_{\varepsilon\in\mathcal{I}}\kappa_{\varepsilon}=0$;

• •

For all $f\in\mathbb{X}$ and all $\varepsilon\in\mathcal{I}$ we have

[TABLE]

Then the true solution $f^{\dagger}$ satisfies the variational source condition (18) with the index function

[TABLE]

4.2 Tools for the verification of variational source conditions for $T$

In order to verify the VSC for our forward operator $T$, we first look at the decomposition in Proposition 3.1. The fact that ${\mathcal{F}^{(1)}},{\mathcal{F}^{(1,*2)}}$ and $\mathcal{G}$ are all unitary operators implies that, in order to analyze the properties of the operator $T$, it suffices to analyze the properties of the multiplication operator $M_{m}$ as defined in (14). We write the forward problem $T\rho^{\dagger}=y^{\mathrm{obs}}$ as

[TABLE]

where $f^{\dagger}$ is defined as

[TABLE]

and will verify the three conditions in Lemma 4.5 for the multiplication operator $M_{m}$ and the true solution $f^{\dagger}$.

For any $\varepsilon>0$, we define the sublevel sets

[TABLE]

with $\Omega:=[-\pi,\pi]\times\mathds{Z}$ such that $|I_{\varepsilon}|=\sum_{k\in\mathds{Z}}|I_{k,\varepsilon}|$. We choose $P_{\varepsilon}$ as orthogonal projections from $L^{2}(\Omega)$ to $L^{2}(\Omega\setminus I_{\varepsilon})$. Obviously, $P_{\varepsilon}$ can be written as a multiplication operator

[TABLE]

with the characteristic function $\chi_{I_{\varepsilon}}$ of $I_{\varepsilon}$.

To bound the sizes of the sublevel sets $I_{k,\varepsilon}$ of $m(\cdot,k)$, we recall some properties of the Bessel functions $J_{k}$ which can be found in [20, Ch. 9]:

$J_{-k}(z)=(-1)^{k}J_{k}(z)$. Hence, it suffices to look at the case $k\geq 0$. 2. 2.

For any $k>0$, $J_{k}$ has a zero of order $k$ at $z=0$. Around this zero it has the asymptotic behavior $J_{k}(z)=\frac{z^{k}}{2^{k}k!}(1+\mathcal{O}(|z|))$ as $|z|\to 0$. 3. 3.

$J_{k}$ also possesses an infinite number of simple zeroes $0<j_{k,1}<j_{k,2}<\dots$ on the positive real axis, and $\lim_{l\to\infty}j_{k,l}=\infty$ for all $k\in\mathds{Z}$. 4. 4.

The positions of the positive zeroes tend to infinity as $k\to+\infty$, i.e. $\lim_{k\to+\infty}j_{k,1}=+\infty$.

These properties translate into the following facts on $m(\cdot,k)$ and its zero set $I_{k,0}$:

Lemma 4.6**.**

Let $m(\varphi,k)$ be defined as in (14). Then we have the following:

$|m(\varphi,-k)|=|m(\varphi,k)|$* and $I_{k,0}=I_{-k,0}$ for all $k\in\mathds{Z}$ and $\varphi\in[-\pi,\pi]$.* 2. 2.

$0\in I_{k,0}$* for all $k>0$, and $|m(\varphi,k)|=\frac{|g_{\omega}|^{k}\varphi^{k}}{k!}(1+\mathcal{O}(|\varphi|))$ as $|\varphi|\to 0$.* 3. 3.

$I_{k,0}\setminus\{0\}$* is a finite set of simple zeros of $m(\cdot,k)$ for all $k\in\mathds{Z}$.* 4. 4.

There exists $K(g_{\omega})\in\mathds{N}$ such that $I_{k,0}=\{0\}$ for all $|k|\geq K(g_{\omega})$.

Proof.

In view of the expression (14) for $m$, the first three statements are immediate consequences of the first three properties of $J_{k}(z)$. For the last statement, note that $I_{k,0}=\{0\}\Leftrightarrow j_{k,1}\leq 4|g_{\omega}|$. As $\lim_{k\to+\infty}j_{k,1}=+\infty$, there exists $K(g_{\omega})$ such that $j_{k,1}>4|g_{\omega}|$ for all $|k|\geq K(g_{\omega})$. ∎

We further have the following properties of $I_{k,\varepsilon}$ as $\varepsilon\to 0$:

Lemma 4.7**.**

For all $k\in\mathds{Z}$ there exists $\varepsilon_{0}(k)>0$ such that for all $0<\varepsilon\leq\varepsilon_{0}(k)$ each connected component of $I_{k,\varepsilon}$ contains exactly one point of $I_{k,0}$. 2. 2.

For $k\neq 0$, the connected component of $I_{k,\varepsilon}$ containing [math] has size $\frac{k}{e|g_{\omega}|}\varepsilon^{1/|k|}(1+\mathcal{O}(1))$ as $\varepsilon\to 0$. 3. 3.

For all $|k|<K(g_{\omega})$, the finitely many connected components of $I_{k,\varepsilon}$ not containing [math] are of size $\mathcal{O}(\varepsilon)$ as $\varepsilon\to 0$. 4. 4.

There exists a constant $C_{I}>0$ depending only on $|g_{\omega}|$ such that

[TABLE]

Proof.

1.) Since $m(\cdot,k)$ has only finitely many zeros of finite order, there exists $\varepsilon_{0}(k)$ such that $m(\cdot,k)$ is monotonic (or even and monotonic on $[0,\pi]$ in case of $m(\cdot,2k)$ around $\varphi=0$) on all connected components of $\Omega_{\varepsilon_{0}(k)}$ viewed as subset of $\mathds{R}/(2\pi\mathds{Z})$. Since $[-\pi,\pi]\setminus\Omega_{\varepsilon_{0}(k)}$ is bounded, $|m(\cdot,k)|$ attains its infimum on the closure of this set, and the infimum is positive. Thus we obtain the claim by possibly reducing $\varepsilon_{0}(k)$.

2,3) This follows from Lemma 4.6, parts 2 and 3, respectively. In part 2 we also use the Stirling approximation $(k!)^{1/k}=\frac{k}{e}\left(1+o(1)\right)$ as $k\to\infty$.

4. This follows from parts 1–3. ∎

Moreover, we need a bound on the supremum norm of $f^{\dagger}$:

Lemma 4.8**.**

Let $\rho$ be an arbitrary density matrix, i.e. $\rho$ is self-adjoint and positive semi-definite, with trace equal to 1. Let the operators ${\mathcal{F}^{(1)}}$ and $\mathcal{G}$ be defined as in (LABEL:FourierTransformF1) and (15), respectively. Then

[TABLE]

Proof.

Using the definitions of ${\mathcal{F}^{(1)}}$ and $\mathcal{G}$, we can see that

[TABLE]

for all $\varphi\in[-\pi,\pi]$ and all $k\in\mathds{Z}$. The fact that $\rho$ is positive semi-definite implies that for all $n,k\in\mathds{Z}$, the principal submatrix

[TABLE]

is also positive semi-definite. Therefore, calculating the determinant of this submatrix yields that

[TABLE]

so we conclude using Young’s inequality that

[TABLE]

holds for all $n,k\in\mathds{Z}$. This in turn implies that

[TABLE]

which concludes the proof. ∎

Note that the third condition in Lemma 4.5 reduces to $\langle P_{\varepsilon}f^{\dagger},\tilde{f}\rangle\leq\sigma_{\varepsilon}\|M_{m}\tilde{f}\|$ for all $\tilde{f}\in L^{2}(\Omega)$. Using the Cauchy-Schwarz inequality and elementary estimates this condition can easily be verified with $\sigma_{\varepsilon}=\varepsilon^{-1}\|f^{\dagger}\|_{L^{2}}$. However, the following more elaborate argument along the lines of [13, Theorem 3.1] gives a sharper bound, which is optimal in some sense (see *****):

Lemma 4.9**.**

Suppose the first condition in Lemma 4.5 holds true for the projection operators defined in (23) and that $\varepsilon\mapsto\kappa_{\varepsilon}\varepsilon^{\nu-1}$ is decreasing for some $\nu\in(0,1)$. Then the third condition holds true with

[TABLE]

Proof.

Due to the inequality

[TABLE]

we have to show that

[TABLE]

Introducing the function

[TABLE]

which by assumption is bounded by $\mu_{f^{\dagger}}(\varepsilon)\leq\kappa_{\varepsilon}^{2}$, we obtain

[TABLE]

using a partial integration in the last line and the fact that $\mu_{f^{\dagger}}(t)=\|f^{\dagger}\|_{L^{2}}^{2}$ for $t>\|m\|_{L^{\infty}}$ such that $\lim_{t\to\infty}t^{-2}\mu_{f^{\dagger}}(t)=0$. Now we use the assumption that $t\mapsto\kappa_{t}^{2}t^{2\nu-2}$ is decreasing to estimate

[TABLE]

This completes the proof of (25). ∎

4.3 Proofs of Theorems 2 and 3

Proposition 4.10** (Hölder VSC for Theorem 2).**

Under the assumptions of Theorem 2 the matrix $\rho^{\dagger}$ satisfies a variational source condition (18) with the index function $\psi(\tau)=C\min\left(\tau^{\frac{1}{1+2k_{0}}},\sqrt{\tau}\right)$ and some $C>0$ depending only on $k_{0}$ and $|g_{\omega}|$.

Proof.

From the band-limited assumption of $\rho^{\dagger}$, we know $f^{\dagger}(\varphi,k)=0$ for $|k|>k_{0}$. Therefore, it follows from Lemmas 4.8 and 4.7 that

[TABLE]

for all $\varepsilon\leq\varepsilon_{0}:=\min\{\varepsilon_{0}(0),\dots,\varepsilon_{0}(k_{0})\}$ and some constant $\tilde{C}>0$ depending only on $k_{0}$ and $|g_{\omega}|$.

Thus the first and second conditions of Lemma 4.5 are satisfied for $\kappa_{\varepsilon}=\tilde{C}^{1/2}\varepsilon^{\frac{1}{2k_{0}}}$. Note that the function $\varepsilon\mapsto\kappa_{\varepsilon}\varepsilon^{-2/3}$ is decreasing as $k_{0}\geq 1$. Using Lemma 4.9 with $\nu=1/3$ we deduce that the third condition holds true for $\sigma_{\varepsilon}=\frac{1}{\sqrt{2}}\varepsilon^{-1}\kappa_{\varepsilon}$. This implies a variational source condition with index function

[TABLE]

If we choose \varepsilon=\min\Big{(}\tau^{k_{0}/(1+2k_{0})},\varepsilon_{0}\Big{)}, we obtain $\psi(\tau)=C_{1}\tau^{1/(1+2k_{0})}$ for $\tau<\varepsilon_{0}^{(1+2k_{0})/k_{0}}$ and $\psi(\tau)=C_{2}\sqrt{\tau}+C_{3}$ else with positive constants $C_{1}$, $C_{2}$, $C_{3}$ depending only on $|g_{\omega}|$ and $k_{0}$. This yields the assertion. ∎

Proof of Theorem 2.

The statement follows from the Hölder-type variational source condition in Proposition 4.10. In particular for the stability estimate (6) we note that $\operatorname{trace}(\rho)=1$ and $\rho\geq 0$ imply that all eigenvalues $\lambda_{j}$ of $\rho$ lie in the interval $[0,1]$, and hence $\|\rho\|_{l^{2}(\mathbb{Z}^{2})}^{2}=\sum_{j=0}^{\infty}\lambda_{j}^{2}\leq\sum_{j=0}^{\infty}\lambda_{j}=\operatorname{trace}(\rho)=1$. Hence, the stability estimate (6) follows from Remark 4.4 as only the behavior of $\psi(\tau)$ of small $\tau$ is relevant. The Hölder-type convergence rate (7) follows from Proposition 4.2. ∎

Proposition 4.11** (VSC for polynomial decay).**

Under the assumptions of Theorem 3, case (8a) the matrix $\rho^{\dagger}$ satisfies a logarithmic variational source condition (18) with

[TABLE]

for some $C>0$ depending only on $C_{\rho}$, $\mu$ and $|g_{\omega}|$.

Proof.

The decay condition on $\rho^{\dagger}$ implies corresponding bounds on $f^{\dagger}$:

[TABLE]

Therefore, considering that $|I_{0,\varepsilon}|=\mathcal{O}(\varepsilon)$, we have

[TABLE]

for some $\tilde{C_{\rho}}>0$.

We will utilize both upper bounds of $\left|I_{k,\varepsilon}\right|$ from Lemma 4.7. As $\lim_{k\to\infty}k\varepsilon^{1/k}=\infty$, we get the bound $\left|I_{k,\varepsilon}\right|\leq 2\pi$ for large $k$. For some cut-off index $k_{0}(\varepsilon)$ to be determined later, we bound the sum by

[TABLE]

Using the logarithmic derivative $\frac{d}{dk}\ln(k^{-\mu}\varepsilon^{1/k})=-4k^{-1}\mu-k^{-2}\ln\varepsilon$, it can be seen that $k\mapsto k^{-\mu}\varepsilon^{1/k}$ is increasing on the interval $\left[0,k_{1}(\varepsilon)\right]$ with $k_{1}(\varepsilon):=\frac{1}{4\mu}\ln\frac{1}{\varepsilon}$. Therefore, as long as $k_{0}(\varepsilon)\leq k_{1}(\varepsilon)$, there exist constants $\varepsilon_{0}$ and $C>0$ depending only on $C_{\rho}$, $|g_{\omega}|$, and $\mu$ such that

[TABLE]

for $\varepsilon\leq\varepsilon_{0}$. We choose $k_{0}$ such that both terms on the right hand side are approximately equal, i.e. $k_{0}\varepsilon^{1/k_{0}}\approx 1$ or equivalently $k_{0}\ln k_{0}\approx\ln\frac{1}{\varepsilon}$. Solving $k_{0}\ln k_{0}=y$ for $k_{0}$ yields the asymptotic relation

[TABLE]

Therefore, we set $k_{0}(\varepsilon):=\lfloor\frac{-\log\varepsilon}{\log(-\log\varepsilon)}\rfloor$. Note that $k_{0}(\varepsilon)\leq k_{1}(\varepsilon)$ for $\varepsilon$ sufficiently small. This yields $\|(I-P_{\varepsilon})f^{\dagger}\|_{L^{2}}\leq\kappa_{\varepsilon}$ with

[TABLE]

and the first and second conditions of Lemma 4.5 are satisfied. By Lemma 4.9 the third condition holds true with $\sigma_{\varepsilon}=\kappa_{\varepsilon}/\varepsilon$. Therefore, Lemma 4.5 yields a VSC with $\psi(\tau)=\inf_{0<\varepsilon\leq\varepsilon_{0}}[\kappa_{\varepsilon}\sqrt{\tau}/\varepsilon+\varepsilon^{2}]$. Choosing $\varepsilon=\min(\tau^{1/3},\varepsilon_{0})$, the first term is asymptotically neglectible against the second, and we obtain a VSC with $\psi(\tau)=\mathcal{O}\left(\kappa_{\tau}^{2}\right)$. This completes the proof. ∎

Proposition 4.12** (VSC for exponential decay).**

Under the assumptions of Theorem 3, case (8b) the matrix $\rho^{\dagger}$ satisfies a logarithmic variational source condition (18) with $\psi(\tau)=Ce^{-2\sqrt{(-\log b)(-\log\tau)}}$ for some $C>0$ depending only on $C_{\rho}$, $b$ and $|g_{\omega}|$.

Proof.

The decay condition on $\rho^{\dagger}$ implies corresponding bounds on $f^{\dagger}$:

[TABLE]

Therefore, again noting that $|I_{0,\varepsilon}|=\mathcal{O}(\varepsilon)$, we have

[TABLE]

for some $\tilde{C_{\rho}}>0$. We will utilize both upper bounds of $\left|I_{k,\varepsilon}\right|$ from Lemma 4.7. As $\lim_{k\to\infty}k\varepsilon^{1/k}=\infty$, we use the trivial bound $\left|I_{k,\varepsilon}\right|\leq 2\pi$ for large $k$. We choose a cut-off $k_{0}(\varepsilon)=\lfloor\sqrt{-\log\varepsilon/(-2\log b)}\rfloor$ for $0<\epsilon<1$ and obtain

[TABLE]

for some generic constant $C>0$ that depends only on $b$, $C_{\rho}$ and $|g_{\omega}|$. Taking the logarithm of the first two terms shows that for our choice of $k_{0}(\varepsilon)$ both logarithms equal $-\sqrt{2(-\log b)(-\log\varepsilon)}$. This shows that

[TABLE]

Therefore, the first and second conditions of Lemma 4.5 are satisfied for $\kappa_{\varepsilon}=C\exp\left(-\sqrt{\frac{1}{2}(-\log b)(-\log\varepsilon)}\right)$. Note that the function $\epsilon\mapsto\varepsilon^{-1/3}\kappa^{\epsilon}$ is decreasing, Lemma 4.9 with $\nu=\frac{2}{3}$ allows us to deduce that the third condition holds true for $\sigma_{\varepsilon}=\frac{\kappa_{\varepsilon}}{2\varepsilon}$. This implies a variational source condition with index function

[TABLE]

If we choose $\varepsilon=\tau^{\frac{1}{4}}$ for $0<\tau<1$, we obtain

[TABLE]

as $\tau\searrow 0$. This yields the assertion. ∎

Proof of Theorem 3.

We set $\Phi(t)=2\sqrt{\psi(t^{2})}$ with the functions $\psi$ in the variational source conditions of Propositions 4.11 and 4.12. Then the statement follows from Proposition 4.2 and Remark 4.4. ∎

5 Conclusions

We have shown that the data acquired in the SQUIRRELS method (without noise and modelling errors) are indeed sufficient to uniquely determine the unknown electron density matrix. Moreover, we have estimated the intrinsic difficulty (or degree of ill-posedness) of the inverse problem to reconstruct a density matrix from these data under noise. As expected, the answer strongly depends on the type of available a-priori information on the unknown density matrix. If this matrix is band-limited, we obtain Hölder rates, whereas under polynomial decay conditions only logarithmic rates can be shown. For the most realistic exponential decay conditions the rates are in between Hölder and logarithmic rates.

We conjecture that both the stability estimates and the convergence rates are of optimal order under the given a-priori information if $T$ is considered as an operator defined on all bounded, Hermitian matrices. However, it is possible that the positive semidefiniteness constraint, which has a strong regularizing effect in numerical experiments, may be further explointed to improve rates.

Another topic of further research in this direction may be to extend the analysis of this paper to a model involving a continuum of energy states.