<figure><img src="image_1.png" /></figure>

HAL
open science

# A generalized spectral concentration problem and the varying masks algorithm

Erwan Faou, Yoann Le Henaff

<figure><img src="image_2.png" /></figure>

## To cite this version:

Erwan Faou, Yoann Le Henaff. A generalized spectral concentration problem and the varying masks algorithm. 2024. hal-04718024

HAL Id: hal-04718024

https://hal.science/hal-04718024v1

Preprint submitted on 2 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A GENERALIZED SPECTRAL CONCENTRATION PROBLEM AND THE VARYING MASKS ALGORITHM

E. FAOU, Y. LE HÉNAFF

**ABSTRACT.** In this paper we generalize the spectral concentration problem as formulated by Slepian, Pollak and Landau in the 1960s. We show that a generalized version with arbitrary space and Fourier masks is well-posed, and we prove some new results concerning general quadratic domains and gaussian filters. We also propose a more general splitting representation of the spectral concentration operator allowing to construct quasi-modes in some situations. We then study its discretization and we illustrate the fact that standard eigen-algorithms are not robust because of a clustering of eigenvalues. We propose a new alternative algorithm that can be implemented in any dimension and for any domain shape, and that gives very efficient results in practice.

## 1. INTRODUCTION

The spectral concentration problem was posed by Slepian, Landau and Pollak in 1961 [22]. This problems stems from the well known fact that a function $f \in L^2$ with fixed $L^2$ norm cannot be both concentrated in space and in Fourier, owing to the celebrated Heisenberg inequality

$$ \|xf\|_{L^2}^2 \|\xi \hat{f}\|_{L^2}^2 \geq \frac{\|f\|_{L^2}^4}{16\pi^2} $$

where $\hat{f}$ is the Fourier transform of $f$. The spectral concentration problem studied by Slepian, Landau and Pollak was to find the function maximizing the $L^2([-1, 1])$ norm of a function for which its Fourier transform is supported in $[-c, c]$ for some given parameter $c > 0$:

$$ \arg \max \left\{ \frac{\int_{-1}^{1} |f|^2}{\int_{-\infty}^{\infty} |f|^2} \mid \text{supp} \hat{f} \subset [-c, c] \right\}. $$

It can be shown to be equivalent after some rescaling to the explicit eigenvalue problem of finding eigenpairs $(\lambda, \psi)$ of the problem

$$ \lambda \psi(x) = \int_{-1}^{1} \frac{\sin(c(x-y))}{x-y} \psi(y) dy =: (\mathcal{K}\psi)(x), \quad x \in [-1, 1]. $$

In a series of papers [15, 16, 21–23], they gave a very satisfying and efficient answer to the above setting by finding an explicit second order operator $\mathcal{P}$ commuting with $\mathcal{K}$:

$$ (1) \qquad \mathcal{P}\mathcal{K} = \mathcal{K}\mathcal{P}, \qquad \mathcal{P} = -\partial_x(1-x^2)\partial_x + c^2x^2. $$

The eigenvalues of $\mathcal{P}$ are special functions known as the prolate spheroidal wave functions or simply Slepian functions, giving automatically eigenfunctions for the operator $\mathcal{K}$. The analysis of the eigenvalue distribution $\lambda$ has also received a lot of attention, see [13, 17] and the references therein. Note that the previous one-dimensional analysis obviously extends to
higher dimension by tensorization, restricting however the analysis to cartesian products of intervals.

A few years later, Brander and DeFacio [4] showed that when the similar problem of space and Fourier Gaussian filtering is considered, it is also possible to find a commuting differential operator (the quantum harmonic oscillator operator, correctly rescaled) and eigenvectors made of scaled Hermite functions.

Apart from the interval and the Gaussian cases, very few results are known about the existence of a second-order commuting differential operator in other settings. The work of Grünbaum, Longhi and Perlstadt [9] even seems to point towards the fact that it does not exist in general except exceptional situations.

Since the spectral concentration problem has been formulated in 1961, it has found numerous applications in different fields of physics (see e.g. the review by Wang [24]). Karnik et al. even proposed a *Fast Slepian transform* [14], underlining the importance of this problem in applications.

For general domain in higher dimension, the description of the spectrum and eigenvectors remains largely unknown. In [20], Simons and Wang considered a spectral concentration problem that cannot be reduced to the historical problem on the interval or to the gaussian filter problem, and they had to resort to purely numerical solutions. Their numerical experiments seemed to indicate that solving the spectral concentration problem with arbitrary space and Fourier restrictions is not an issue, but this is essentially due to the favorable numerical parameters. In practice, simple experiments show that solving numerically a spectral concentration problem can be a source of dramatic instabilities. This is a purely numerical issue, already present in the interval case, and it boils down to finding eigenvectors of a matrix for which the eigenvalues are very close to each other, forming almost *large clusters* of eigenvalues preventing standard algorithms to catch physically relevant eigenvectors. Typically for domains, the most relevant eigenvectors are associated with eigenvalues extremely close to 1 and it is numerically extremely difficult to classify them at a reasonable cost. The usual method used to circumvent these instabilities is to use a very large number of discretization points, at a cost that becomes rapidly prohibitive in dimensions 2 or 3 – which are the most widely used cases for applications.

The main goals of this paper are the following:

* We define and study a generalized spectral concentration problem in any dimension covering the previous situations and we give some basic properties of the associated spectrum.

* We give some examples where the spectrum can be calculated and estimated, in particular for general *quadratic* domains. We also use a general representation of the spectral concentration operator as a *Strang splitting* operator allowing to have an exact expression of eigenvalues and eigenvectors for Gaussians filters. We also give a method for constructing *quasi-modes* using commutators of the Baker-Campbell-Hausdorff formula for general filters close to the identity.

* Eventually, we propose a new algorithm for approximating the eigenpairs of the spectral concentration operator. This methods that we call the *varying mask* algorithm allows to track relevant eigenvectors by letting the size of the domain vary. We report in one and two-dimensional examples the excellent behavior of this method compared to standard eigendecomposition algorithms.
**Acknowledgment.** The authors would like to thank Pierre Vernaz-Gris for several stimulating discussions of this subject. This work was conducted within the France 2030 program, Centre Henri Lebesgue ANR-11-LABX-0020-01.

1.1. **Notations.** The Fourier transform on $L^2(\mathbb{R}^d)$ and its inverse are chosen respectively as follows:

$$ (2) \qquad \mathcal{F}[f](\xi) := \int_{\mathbb{R}^d} f(x)e^{-i\xi \cdot x}dx, \quad \text{and} \quad \mathcal{F}^{-1}[f](x) := \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} f(\xi)e^{ix \cdot \xi}d\xi. $$

We may use the shorthand $\hat{f}$ to denote $\mathcal{F}[f]$. Some useful classical properties of the Fourier transform are the following: for $f, g \in L^2(\mathbb{R}^d)$,

$$ \bullet \quad (2\pi)^d \int_{\mathbb{R}^d} f(x) \overline{g(x)} dx = \int_{\mathbb{R}^d} \mathcal{F}[f](\xi) \overline{\mathcal{F}[g](\xi)} d\xi, $$

$$ \bullet \quad \mathcal{F}[f(\cdot - a)](\xi) = e^{-i\xi \cdot a} \mathcal{F}[f](\xi), \text{ for } \xi, a \in \mathbb{R}^d, $$

$$ \bullet \quad \mathcal{F}(f * g) = \mathcal{F}(f)\mathcal{F}(g), $$

where $*$ denotes the convolution operator: $(f * g)(x) = \int_{\mathbb{R}^d} f(y)g(x - y)dy$.

The $L^2$ inner product will be denoted by $(\cdot, \cdot)_{L^2}$ and it is taken anti-hermitian in its second argument. $\|f\|_{L^p}$ denote the standard $L^p$ norm of a given function $f$.

We write $\mathbf{A} \in \mathbb{K}^{m \times n}$ to denote a matrix with components in $\mathbb{K}$ with $m$ rows and $n$ columns. For a given function $f \in L^2(\mathbb{R}^d)$ and a given matrix $\mathbf{A} \in \mathbb{R}^{d \times d}$, we use the shorthand $f \circ \mathbf{A}(x) := f(\mathbf{A}x)$.

A ball centered at $c \in \mathbb{R}^d$ with radius $r > 0$ is denoted $B(c, r)$. For any $p \in \mathbb{R}^d$, $\tau_p$ denotes a translation by $p$: $\tau_p(x) := x - p$.

For two integers $m < n \in \mathbb{Z}$, we write $\llbracket m, n \rrbracket = [m, n] \cap \mathbb{Z}$.

Every index $s$ denotes a quantity that is related to the space domain, and every index $f$ denotes a quantity that is related to the Fourier domain.

## 2. DERIVATION OF THE GENERALIZED SPECTRAL CONCENTRATION PROBLEM

Let $m_S, \widehat{m_F} \in L^2(\mathbb{R}^d)$. The function $m_S$ will be called the *space filter*, or *space mask*, and $\widehat{m_F}$ will be called the *Fourier filter*, or *Fourier mask*. We define the following operators: for $f \in L^2(\mathbb{R}^d)$,

$$ (\mathcal{M}_S f)(x) := m_S(x)f(x), \quad \text{and} \quad (\mathcal{M}_F f)(x) := \mathcal{F}^{-1}[\widehat{m_F} \mathcal{F}[f]](x). $$

In this work we will consider a composition of these operators, more specifically $\mathcal{M}_F \mathcal{M}_S$. Most of what follows would also apply if we considered $\mathcal{M}_S \mathcal{M}_F$, with small modifications.

One has

$$ (\mathcal{M}_F \mathcal{M}_S f)(x) = \mathcal{F}^{-1}[\widehat{m_F} \mathcal{F}[m_S f]](x) = (m_F * (m_S f))(x). $$

Owing to the fact that the Fourier transform of the Dirac distribution $\delta_0$ is the identity function, i.e. $\int_{\mathbb{R}^d} e^{i\eta \cdot x} dx = (2\pi)^d \delta_0(\eta)$, direct computations yield

$$ \begin{align} \|\mathcal{M}_F \mathcal{M}_S f\|_{L^2}^2 &= \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} \int_{\mathbb{R}^d} m_S(x) f(x) \overline{m_S(y) f(y)} \left( \int_{\mathbb{R}^d} e^{i\xi \cdot (x-y)} |\widehat{m_F}(\xi)|^2 d\xi \right) dy dx \tag{3} \\ &= (\mathcal{K} f, f)_{L^2}. \end{align} $$

The operator $\mathcal{K}$ is the *generalized concentration operator* on $L^2(\mathbb{R}^d)$ that we have defined by the formula

$$ (4) \qquad \mathcal{K} = \mathcal{M}_S^* \mathcal{M}_F^* \mathcal{M}_F \mathcal{M}_S, \qquad 3 $$
where $\mathcal{M}_S^*$ and $\mathcal{M}_F^*$ are the $L^2$-adjoint of the operators $\mathcal{M}_S$ and $\mathcal{M}_F$ respectively. We thus have

$$ (5) \qquad \left| \begin{array}{l} (\mathcal{K}f)(x) := \overline{m_S(x)} \int_{\mathbb{R}^d} k(x-y)m_S(y)f(y)dy. \\ k(z) = \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} e^{i\xi \cdot z} |\widehat{m_F}(\xi)|^2 d\xi. \end{array} \right. $$

**Lemma 1.** If $m_S, \widehat{m_F} \in L^2(\mathbb{R}^d)$, then for all $f \in L^2(\mathbb{R}^d)$,

$$ (6) \qquad \|\mathcal{K}f\|_{L^2} \le \|\widehat{m_F}\|_{L^2}^2 \|m_S\|_{L^2}^2 \|f\|_{L^2}. $$

*Proof.* It is a consequence of Young and Hölder inequalities which yield

$$ \|\mathcal{M}_F \mathcal{M}_S f\|_{L^2} = \|m_F * (m_S f)\|_{L^2} \le \|m_F\|_{L^2} \|m_S f\|_{L^1} \le \|\widehat{m_F}\|_{L^2} \|m_S\|_{L^2} \|f\|_{L^2} $$

and the definition of $\mathcal{K}$.

<figure><img src="image_4.png" /></figure>

For any function $f \in L^2(\mathbb{R}^d)$, we define the associated concentration ratio:

$$ (7) \qquad \nu(f) := \frac{\|\mathcal{M}_F \mathcal{M}_S f\|_{L^2}^2}{\|f\|_{L^2}^2} = \frac{(\mathcal{K}f, f)_{L^2}}{(f, f)_{L^2}}. $$

It is clear that $\nu \ge 0$ and that $\nu \le \|m_S\|_{L^2}^2 \|\widehat{m_F}\|_{L^2}^2$ using Lemma 1. In this work we are interested in finding the functions that maximize the concentration ratio, i.e. solutions to the eigenvalue problem

$$ (8) \qquad \text{Find } (\lambda, \psi) \in \mathbb{R} \times L^2 \text{ s.t. } \mathcal{K}\psi = \lambda\psi. $$

and to compute efficiently the eigenvectors corresponding to the largest eigenvalues $\lambda$. This problem is equivalent to finding the singular values of the operator $\mathcal{M}_F \mathcal{M}_S$.

**2.1. Properties.** Using classical results (see e.g. [19, Chapter VI]), one gets the following properties:

**Proposition 1.** The concentration operator $\mathcal{K}$ defined in (5) enjoys the following properties:

a) $\mathcal{K}$ is a Hilbert-Schmidt operator, self-adjoint, compact, and positive semi-definite.

b) The countable family $\{\psi_i\}_{i=1}^\infty$ of eigenfunctions of $\mathcal{K}$ is orthonormal for the usual $L^2(\mathbb{R}^d)$ inner product and complete in $L^2(\mathbb{R}^d)$. The associated eigenvalues $\{\lambda_i\}_{i=1}^\infty$ are real, nonnegative, and we can order them so that $\lambda_i \ge \lambda_{i+1} \ge 0$ for $i \ge 1$.

c) The orthonormal basis of eigenfunctions $\{\psi_i\}_{i=1}^\infty$ are critical points for the concentration ratio (7), and can be obtained by the successive maximization problems

$$ (9) \qquad \lambda_n = \sup_{\substack{f \in L^2(\mathbb{R}^d) \\ f \in [\psi_1, \dots, \psi_{n-1}]^\perp}} \frac{(\mathcal{K}f, f)_{L^2}}{(f, f)_{L^2}}, $$

where $[\psi_1, \dots, \psi_{n-1}]^\perp := \{u \in L^2(\mathbb{R}^d) : (u, \psi_i)_{L^2} = 0, i = 1, \dots, n-1\}$.

d) For large $n$, $\lambda_n = o(n^{-1/2})$.

e) Suppose $|\widehat{m_F}|^2$ is even, and $m_S$ is real, then $\mathcal{K}$ is real-valued for real inputs.

*Proof.* The fact that $\mathcal{K}$ is a Hilbert-Schmidt operator is a consequence of Lemma 1 and the compactness is due to [19, Theorem VI.22]. The other properties of a) are easily derived from (4).
To obtain b), we use the Hilbert-Schmidt Theorem (see [19, Chapter VI]) which can be applied since we have just shown that $\mathcal{K}$ is a self-adjoint and compact operator.

To prove c) and (9), we recall the min-max theorem:

$$ \lambda_n = \min_{\varphi_1, \dots, \varphi_{n-1}} \max_{\substack{f \in L^2(\mathbb{R}^d) \\ f \in [\varphi_1, \dots, \varphi_{n-1}]^\perp}} \frac{(\mathcal{K}f, f)_{L^2}}{(f, f)_{L^2}}, $$

see for instance [5, Thm. 6.23] or [19, Thm. VI.15]. Now, if $f \in [\psi_1, \dots, \psi_{n-1}]^\perp$, we have $f = \sum_{i \ge n} a_i \psi_i$ where $\psi_i$ is the i-th eigenfunction of $\mathcal{K}$. For such a function $f$, we get

$$ \frac{(\mathcal{K}f, f)_{L^2}}{(f, f)_{L^2}} = \frac{\sum_{i=n}^{+\infty} |a_i|^2 \lambda_i(\mathcal{K})}{\sum_{i=n}^{+\infty} |a_i|^2} \le \lambda_n(\mathcal{K}), $$

and the fact that equality is attained with $f = \psi_n$ shows (9).

The point d) follows from [19, Theorem VI.22], which states that

$$ \sum_{n=1}^{\infty} \lambda_n^2 < \infty. $$

This series is convergent only if $\lambda_n = o(n^{-1/2})$ for large $n$.

About the point e), it suffices to show that the inverse Fourier transform of $|\widehat{m_F}|^2$ is real. We have

$$ \int_{\mathbb{R}^d} |\widehat{m_F}(\xi)|^2 e^{i\xi \cdot (y-x)} d\xi = \int_{\mathbb{R}^d} |\widehat{m_F}(\xi)|^2 (\cos(\xi \cdot (y-x)) + i \sin(\xi \cdot (y-x))) d\xi. $$

When $|\widehat{m_F}|^2$ is even, the complex part vanishes as the integral of an odd function. Thus, only the real part remains. $\square$

Some properties of the eigenfunctions $\{\psi_i\}_{i=1}^\infty$ can easily be obtained:

**Lemma 2 (Symmetries).** Suppose there is an orthogonal matrix $\mathbf{S} \in \mathbb{R}^{d \times d}$ and $\alpha \in \mathbb{C}$, $|\alpha| = 1$, such that $m_S \circ \mathbf{S} = \alpha m_S$ and $|\widehat{m_F} \circ \mathbf{S}| = |\widehat{m_F}|$. Then

$$ \mathcal{K}(\psi \circ \mathbf{S}) = (\mathcal{K}(\psi)) \circ \mathbf{S}. $$

In particular, if $\psi$ is an eigenfunction of $\mathcal{K}$ associated to an eigenvalue $\lambda$ of multiplicity one, then there exists $\beta \in \mathbb{C}$, $|\beta| = 1$, such that

$$ \psi \circ \mathbf{S} = \beta \psi. $$

*Proof.* It follows from straightforward computations. First of all, since $\mathbf{S}$ is orthogonal, $|\det \mathbf{S}| = 1$. We now compute

$$ (\mathcal{K}\psi)(\mathbf{S}x) = \int_{\mathbb{R}^d} m_S(y) \overline{m_S(\mathbf{S}x)} \mathcal{F}^{-1}[|\widehat{m_F}|^2] (\mathbf{S}x - y)\psi(y)dy. $$

The change of variables $y = \mathbf{S}\tilde{y}$ yields

$$ (\mathcal{K}\psi)(\mathbf{S}x) = \int_{\mathbb{R}^d} m_S(\mathbf{S}\tilde{y}) \overline{m_S(\mathbf{S}x)} \mathcal{F}^{-1}[|\widehat{m_F}|^2] (\mathbf{S}x - \mathbf{S}\tilde{y})\psi(\mathbf{S}\tilde{y})d\tilde{y}, $$

where we have used $|\det \mathbf{S}| = 1$. Owing to the assumption $m_S \circ \mathbf{S} = \alpha m_S$, $|\alpha| = 1$,

$$ (\mathcal{K}\psi)(\mathbf{S}x) = \int_{\mathbb{R}^d} m_S(\tilde{y}) \overline{m_S(x)} \mathcal{F}^{-1}[|\widehat{m_F}|^2] (\mathbf{S}x - \mathbf{S}\tilde{y})\psi(\mathbf{S}\tilde{y})d\tilde{y}. $$
It only remains to show that $\mathcal{F}^{-1}[|\widehat{m_F}|^2](\mathbf{S}x - \mathbf{S}\tilde{y}) = \mathcal{F}^{-1}[|\widehat{m_F}|^2](x - \tilde{y})$. Letting $\xi = \mathbf{S}\tilde{\xi}$,

$$ \begin{aligned} \mathcal{F}^{-1}[|\widehat{m_F}|^2](\mathbf{S}x - \mathbf{S}\tilde{y}) &= \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} |\widehat{m_F}(\xi)|^2 e^{i\xi \cdot \mathbf{S}(x-\tilde{y})} d\xi \\ &= \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} |\widehat{m_F}(\mathbf{S}\tilde{\xi})|^2 e^{i(\mathbf{S}\tilde{\xi}) \cdot \mathbf{S}(x-\tilde{y})} d\tilde{\xi}. \end{aligned} $$

We have again used the fact that $|\det \mathbf{S}| = 1$. Since $\mathbf{S}$ is orthogonal, $\mathbf{S}^T \mathbf{S} = I$, therefore $(\mathbf{S}\tilde{\xi}) \cdot \mathbf{S}(x - \tilde{y}) = \tilde{\xi} \cdot (x - \tilde{y})$ and

$$ \mathcal{F}^{-1}[|\widehat{m_F}|^2](\mathbf{S}x - \mathbf{S}\tilde{y}) = \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} |\widehat{m_F}(\tilde{\xi})|^2 e^{i\tilde{\xi} \cdot (x-\tilde{y})} d\tilde{\xi}, $$

where we have used the assumption $|\widehat{m_F} \circ \mathbf{S}| = |\widehat{m_F}|$. We finally obtain

$$ (K\psi)(\mathbf{S}x) = \int_{\mathbb{R}^d} m_S(\tilde{y}) \overline{m_S(x)} \mathcal{F}^{-1}[|\widehat{m_F}|^2](x - \tilde{y}) \psi(\mathbf{S}\tilde{y}) d\tilde{y}, $$

which is exactly $(K\psi) \circ \mathbf{S} = K(\psi \circ \mathbf{S})$. If $\psi$ is an eigenfunction associated to an eigenvalue $\lambda$ of multiplicity one, so is $\psi \circ \mathbf{S}$. Therefore, they must agree up to some multiplicative constant. Due to the orthogonality of $\mathbf{S}$ they have the same $L^2(\mathbb{R}^d)$ norm, so that constant must have modulus one. $\square$

Let us consider the case where $\lambda$ is a multiple eigenvalue of multiplicity $p \in \mathbb{N}^*$. Write $\phi_1, \dots, \phi_p$ the eigenfunctions of $K$ associated to $\lambda$. For any $i = 1, \dots, p$, the same computations as above yield that $\phi_i \circ \mathbf{S}$ is an eigenfunction associated to $\lambda$. Therefore, we can only decompose

$$ \phi_i \circ \mathbf{S} = \sum_{j=1}^{p} b_{ij} \phi_j, \quad b_{ij} \in \mathbb{C} \text{ with } \sum_{j=1}^{p} |b_{ij}| = 1. $$

This formula can have some applications, for example in presence of rotational or asymmetric invariances. Indeed, assume for instance that $\mathbf{S}^k = \text{Id}$ for some $k$ (which is the case for discs or polygons with axis of symmetries in 2D or for 3D axisymmetric domains, with invariance by angular rotation of angle $\frac{2\pi}{k}$). Then the previous relation implies that $\mathbf{B}^k = \text{Id}$ which implies that the eigenfunctions can be sorted with respect to the eigenvalues $\alpha_n = e^{\frac{2i\pi n}{k}}$, $n = 1, \dots, k$ of $\mathbf{B} = (b_{i,j})_{i,j}$. This will be easily observed in the 2D examples below (the disc and the cat-head).

## 2.2. Special case of binary masks.
The situation with binary masks has some interesting properties. For two (smooth enough) domains $\Omega_S, \Omega_F \subset \mathbb{R}^d$ if we take $m_S(x) = \mathbf{1}_{\Omega_S}$ and $\widehat{m_F} = \mathbf{1}_{\Omega_F}$, the eigenvalue problem associated with the operator (5) can be written (10)

$$ \text{Find } (\lambda, \psi), \text{ s.t. } \lambda\psi(x) = \int_{\Omega_S} k(x-y)\psi(y)dy, \quad x \in \Omega_S, \quad k(z) = \frac{1}{(2\pi)^d} \int_{\Omega_F} e^{i\lambda z} d\lambda. $$

Note that the study of the function $k(z)$ for general domains is a difficult question, see for instance [8].

**Lemma 3 (Translations with binary masks).** Let $\Omega_S, \Omega_F \subset \mathbb{R}^d, p \in \mathbb{R}^d$. The following equivalences hold:
• $(\lambda, \psi)$ is an eigenpair of the concentration operator associated to masks $m_S = \mathbf{1}_{\Omega_S+p}$ and $\widehat{m_F} = \mathbf{1}_{\Omega_F}$ iff $(\lambda, \psi \circ \tau_{-p})$ is an eigenpair of the concentration operator associated to masks $m_S = \mathbf{1}_{\Omega_S}$ and $\widehat{m_F} = \mathbf{1}_{\Omega_F}$;

• $(\lambda, \psi)$ is an eigenpair of the concentration operator associated to masks $m_S = \mathbf{1}_{\Omega_S}$ and $\widehat{m_F} = \mathbf{1}_{\Omega_F}$ iff $(\lambda, x \mapsto \psi(x)e^{-ip \cdot x})$ is an eigenpair of the concentration operator associated to masks $m_S = \mathbf{1}_{\Omega_S}$ and $\widehat{m_F} = \mathbf{1}_{\Omega_F}$.

*Proof.* Let us start with the first claim, and consider $(\lambda, \psi)$ an eigenpair of the concentration operator associated to masks $m_S = \mathbf{1}_{\Omega_S+p}$ and $\widehat{m_F} = \mathbf{1}_{\Omega_F}$. We have

$$ \lambda\psi(x) = (2\pi)^{-d}\mathbf{1}_{\Omega_S+p}(x) \int_{\Omega_S+p} \psi(y) \int_{\Omega_F} e^{i\xi \cdot (x-y)} d\xi dy \\ \iff \lambda\psi(x) = (2\pi)^{-d}\mathbf{1}_{\Omega_S}(\tau_p(x)) \int_{\Omega_S+p} \psi(y) \int_{\Omega_F} e^{i\xi \cdot (\tau_p(x)-\tau_p(y))} d\xi dy. $$

Write $u := \tau_p(x)$ and $v := \tau_p(y)$, then $dv = dy$ and

$$ \lambda\psi(\tau_{-p}(u)) = (2\pi)^{-d}\mathbf{1}_{\Omega_S}(u) \int_{\Omega_S} \psi(\tau_{-p}(v)) \int_{\Omega_F} e^{i\xi \cdot (u-v)} d\xi dv. $$

In other words, $(\lambda, \psi \circ \tau_{-p})$ is an eigenpair of the concentration operator associated to masks $m_S = \mathbf{1}_{\Omega_S}$ and $\widehat{m_F} = \mathbf{1}_{\Omega_F}$.

We now proceed to proving the second claim. Consider $(\lambda, \psi)$ an eigenpair of the concentration operator associated to masks $m_S = \mathbf{1}_{\Omega_S}$ and $\widehat{m_F} = \mathbf{1}_{\Omega_F+p}$, then

$$ \lambda\psi(x) = (2\pi)^{-d}\mathbf{1}_{\Omega_S}(x) \int_{\Omega_S} \psi(y) \int_{\Omega_F+p} e^{i\xi \cdot (x-y)} d\xi dy. $$

Write $\xi = \zeta + p$ for $\zeta \in \Omega_F$, then $d\xi = d\zeta$ and

$$ \lambda\psi(x) = (2\pi)^{-d}\mathbf{1}_{\Omega_S}(x) \int_{\Omega_S} \psi(y) \int_{\Omega_F} e^{i(\zeta+p) \cdot (x-y)} d\zeta dy \\ \iff \lambda\psi(x)e^{-ip \cdot x} = (2\pi)^{-d}\mathbf{1}_{\Omega_S}(x) \int_{\Omega_S} \psi(y)e^{-ip \cdot y} \int_{\Omega_F} e^{i\zeta \cdot (x-y)} d\zeta dy. $$

In other words, $(\lambda, x \mapsto \psi(x)e^{-ip \cdot x})$ is an eigenpair of the concentration operator associated to masks $m_S = \mathbf{1}_{\Omega_S}$ and $\widehat{m_F} = \mathbf{1}_{\Omega_F}$. $\square$

The following result now generalizes the scaling invariance of the problem, implying in 1D the fact that the eigenfunctions depend only on the product of the size of the space and Fourier intervals:

**Lemma 4** (Affine transformations with binary masks). Let $\mathbf{A} \in \mathbb{R}^{d \times d}$ an invertible matrix, $\Omega_S, \Omega_F \subset \mathbb{R}^d$, and write

$$ \mathbf{A}\Omega_F := \{\mathbf{A}z : z \in \Omega_F\}. $$

Let $(\lambda, \psi)$ an eigenpair of the concentration operator associated to binary masks $m_S = \mathbf{1}_{\Omega_S}$ and $\widehat{m_F} = \mathbf{1}_{\mathbf{A}\Omega_F}$. Then, $(\lambda, \psi \circ \mathbf{A}^{-T})$ is an eigenpair of the concentration operator associated to binary masks $m_S = \mathbf{1}_{\mathbf{A}^T\Omega_S}$ and $\widehat{m_F} = \mathbf{1}_{\Omega_F}$. The converse is also true.
Proof. An eigenpair $(\lambda, \psi)$ of the concentration operator associated to masks $m_S = \mathbf{1}_{\Omega_S}$ and $\widehat{m_F} = \mathbf{1}_{\mathbf{A}\Omega_F}$ satisfies the following equality:

$$ \lambda\psi(x) = (2\pi)^{-d}\mathbf{1}_{\Omega_S}(x) \int_{\Omega_S} \psi(y) \int_{\mathbf{A}\Omega_F} e^{i\xi\cdot(x-y)}d\xi dy. $$

Let $\xi := \mathbf{A}\zeta$ for $\zeta \in \Omega_F$, then

$$ \begin{aligned} \lambda\psi(x) &= (2\pi)^{-d}\mathbf{1}_{\Omega_S}(x) \int_{\Omega_S} \psi(y) \int_{\Omega_F} e^{i(\mathbf{A}\zeta)\cdot(x-y)} |\det \mathbf{A}| d\zeta dy \\ &= (2\pi)^{-d}\mathbf{1}_{\Omega_S}(x) \int_{\Omega_S} \psi(y) \int_{\Omega_F} e^{i\zeta\cdot(\mathbf{A}^T(x-y))} |\det \mathbf{A}| d\zeta dy. \end{aligned} $$

Let $u := \mathbf{A}^T x$ and $v := \mathbf{A}^T y$, then

$$ \lambda\psi(\mathbf{A}^{-T}u) = (2\pi)^{-d}\mathbf{1}_{\mathbf{A}^T\Omega_S}(u) \int_{\mathbf{A}^T\Omega_S} \psi(\mathbf{A}^{-T}v) \int_{\Omega_F} e^{i\zeta\cdot(u-v)}d\zeta dv. $$

Hence, $(\lambda, \psi \circ \mathbf{A}^{-T})$ is an eigenpair of the concentration operator associated to masks $m_S = \mathbf{1}_{\mathbf{A}^T\Omega_S}$ and $\widehat{m_F} = \mathbf{1}_{\Omega_F}$. $\square$

### 3. NEW EXAMPLES

**3.1. Quadratic domains.** In the literature, the only examples of quadratic operator $\mathcal{P}$ commuting with $\mathcal{K}$ are the operator (1) in the case of interval, as well as the case of balls obtained by Slepian. We give below a new general result for domains delimited by general quadrics. We define the family of $d$-dimensional quadratic domains as follows:

$$ (11) \qquad Q(c, a, b) := \left\{ x \in \mathbb{R}^d : \sum_{m=1}^{d} (x_m - c_m)^2 a_m \le b \right\}, \quad a, c \in \mathbb{R}^d, b \in \mathbb{R}. $$

Note that in full generality, we allow negative coefficients, but that in the case where $a_m > 0$, these domains are ellipsoidal domains. In the following we will assume these ellipses to be centered, i.e. $c = 0$, to simplify the calculations. This is done without loss of generality using Lemma 3.

**Proposition 2.** Consider the $d$-dimensional concentration problem where the space domain is restricted to $\Omega_S := Q(0, a, b)$ and the Fourier domain to $\Omega_F := Q(0, \alpha, \beta)$, for some $a, \alpha \in \mathbb{R}^d$ and $b, \beta \in \mathbb{R}$. Let $\mathcal{K}$ be the concentration operator associated to masks $m_S = \mathbf{1}_{\Omega_S}$ and $\widehat{m_F} = \mathbf{1}_{\Omega_F}$, see (10). Then, there exists a second-order differential operator $\mathcal{P}$, self-adjoint on $L^2(\Omega_S)$, which commutes with the concentration operator $\mathcal{K}$. It is given by

$$ (12) \qquad \mathcal{P}(x, \nabla_x) = \nabla_x \cdot (\mathbf{A}(x)\nabla_x) + C(x), $$

where

$$ (13) \qquad \mathbf{A}(x) := \operatorname{diag} \left\{ \alpha_m \left( \sum_{n=1}^{d} a_n x_n^2 - b \right) \right\}_{\substack{m=1 \\ 8}}^{d} \quad \text{and} \quad C(x) = \beta x \cdot \operatorname{diag} \{a_n\}_{n=1}^{d} x. $$
Proof. Let $\mathcal{P}$ be a self-adjoint second-order differential operator of the form (12) where $x \mapsto \mathbf{A}(x)$ is a matrix-valued real function such that $\mathbf{A}(x) = \mathbf{A}(x)^T$ on $\Omega_S$, and $x \mapsto C(x)$ is a real scalar function. Note that the condition

$$ (14) \qquad \mathbf{A} = 0 \quad \text{on } \partial\Omega_S = \left\{ x \in \mathbb{R}^d : \sum_{n=1}^d a_n x_n^2 = b \right\}, $$

ensures that $\mathcal{P}$ is self-adjoint on $L^2(\Omega_S)$. Hence we have with the notation (10),

$$ (\mathcal{K}\mathcal{P}f)(x) = \int_{\Omega_S} [\mathcal{P}(y, \nabla_y)f(y)] k(x-y)dy = \int_{\Omega_S} f(y)[\mathcal{P}(y, \nabla_y)k(x-y)] dy. $$

The commutation relation $\mathcal{P}\mathcal{K} = \mathcal{K}\mathcal{P}$ is thus implied by the condition

$$ (15) \qquad \forall (x, y) \in \Omega_S^2, \quad \mathcal{P}(x, \nabla_x)k(x-y) = \mathcal{P}(y, \nabla_y)k(x-y). $$

Letting $\varphi(x, \xi) := e^{i\xi \cdot x}$, we have

$$ k(x-y) = \frac{1}{(2\pi)^d} \int_{\Omega_F} \varphi(x, \xi) \overline{\varphi(y, \xi)} d\xi. $$

Owing to

$$ (16) \qquad \nabla_x \varphi(x, \xi) = i\xi \varphi(x, \xi) \quad \text{and} \quad x\varphi(x, \xi) = -i\nabla_x \varphi(x, \xi), $$

we get

$$ \mathcal{P}(x, \nabla_x)\varphi(x, \xi) = \mathcal{P}(-i\nabla_x, i\xi)\varphi(x, \xi). $$

The differential operator $\mathcal{P}$ having real coefficients, we get

$$ \mathcal{P}(y, \nabla_y)\overline{\varphi(y, \xi)} = \overline{\mathcal{P}(y, \nabla_y)\varphi(y, \xi)} = \overline{\mathcal{P}(-i\nabla_x, i\xi)\varphi(y, \xi)}. $$

The commutation relation (15) then writes

$$ (17) \qquad \int_{\Omega_F} \overline{\varphi(y, \xi)} \mathcal{P}(-i \nabla_\xi, i \xi) \varphi(x, \xi) d\xi = \int_{\Omega_F} \varphi(x, \xi) \overline{\mathcal{P}(-i \nabla_\xi, i \xi) \varphi(y, \xi)} d\xi. $$

In other words, we want the differential operator

$$ \mathcal{P}(-i\nabla_x, i\xi) = i\xi \cdot (\mathbf{A}(-i\nabla_x)(i\xi)) + C(-i\nabla_x) $$

to be self-adjoint on $L^2(\Omega_F)$ for the Hermitian inner product.

Let $\mathbf{A}(x)$ be the symmetric matrix defined by the equation (13), and which satisfies the boundary condition (14). One has

$$ \begin{aligned} \mathcal{P}(-i\nabla_\xi, i\xi) &= -\xi \cdot \text{diag} \left\{ \alpha_m \left( -\sum_{n=1}^d a_n \partial_{\xi_n}^2 - b \right) \right\}_{m=1}^d \xi + C(-i\nabla_\xi) \\ &= \sum_{m=1}^d \xi_m \alpha_m \left( \sum_{n=1}^d a_n \partial_{\xi_n}^2 \xi_m + b\xi_m \right) + C(-i\nabla_\xi) \\ &= \sum_{m=1}^d \alpha_m \left( \sum_{n=1}^d a_n \xi_m \partial_{\xi_n}^2 \xi_m + b\xi_m^2 \right) + C(-i\nabla_\xi). \end{aligned} $$

We emphasize that $\partial_{\xi_m}\xi_n$ is actually the operator given for $h \in C^\infty(\mathbb{R}^d)$ by

$$ \partial_{\xi_m} \xi_n h = \partial_{\xi_m} (\xi_n h). $$
In particular, since $\partial_{\xi_m}(\xi_m h) = h + \xi_m \partial_{\xi_m} h$, we have

$$ \partial_{\xi_m} \xi_m - \xi_m \partial_{\xi_m} = 1 \quad \text{and} \quad \partial_{\xi_m} \xi_n = \xi_n \partial_{\xi_m}, \quad m \neq n. $$

Therefore,

$$ \begin{aligned} \sum_{n=1}^{d} a_n \xi_m \partial_{\xi_n}^2 \xi_m &= a_m \xi_m \partial_{\xi_m}^2 \xi_m + \sum_{n \neq m} a_n \xi_m \partial_{\xi_n}^2 \xi_m = a_m \xi_m \partial_{\xi_m}^2 \xi_m + \sum_{n \neq m} a_n \partial_{\xi_n} \xi_m^2 \partial_{\xi_n} \\ &= a_m (\partial_{\xi_m} \xi_m - 1) (\xi_m \partial_{\xi_m} + 1) + \sum_{n \neq m} a_n \partial_{\xi_n} \xi_m^2 \partial_{\xi_n} \\ &= a_m (\partial_{\xi_m} \xi_m^2 \partial_{\xi_m} - \xi_m \partial_{\xi_m} + \partial_{\xi_m} \xi_m - 1) + \sum_{n \neq m} a_n \partial_{\xi_n} \xi_m^2 \partial_{\xi_n} \\ &= \sum_{n=1}^{d} a_n \partial_{\xi_n} \xi_m^2 \partial_{\xi_n}. \end{aligned} $$

Hence we have

$$ \mathcal{P}(-i\nabla_{\xi}, i\xi) = \sum_{n=1}^{d} a_n \sum_{m=1}^{d} \alpha_m \partial_{\xi_n} \xi_m^2 \partial_{\xi_n} + b \sum_{k=1}^{d} \alpha_m \xi_m^2 + C(-i\nabla_{\xi}) $$

We recognize here that

$$ \mathcal{P}(-i\nabla_{\xi}, i\xi) = \nabla_{\xi} \cdot \text{diag} \left\{ a_n \sum_{m=1}^{d} \alpha_m \xi_m^2 \right\}_{n=1}^{d} \nabla_{\xi} + b\xi \cdot \text{diag} \{\alpha_m\}_{m=1}^{d} \xi + C(-i\nabla_{\xi}). $$

By choosing

$$ C(-i\nabla_{\xi}) = -\beta\nabla_{\xi} \cdot \text{diag} \{a_n\}_{n=1}^{d} = -\beta \sum_{n=1}^{d} a_n \partial_{\xi_n}^2, $$

we obtain

$$ \mathcal{P}(-i\nabla_{\xi}, i\xi) = \nabla_{\xi} \cdot \text{diag} \left\{ a_n \left( \sum_{m=1}^{d} \alpha_m \xi_m^2 - \beta \right) \right\}_{n=1}^{d} \nabla_{\xi} + b\xi \cdot \text{diag} \{\alpha_m\}_{m=1}^{d} \xi. $$

It is a self-adjoint operator on $L^2(\Omega_F)$ for the Hermitian inner product since

$$ \text{diag} \left\{ a_n \left( \sum_{m=1}^{d} \alpha_m \xi_m^2 - \beta \right) \right\}_{n=1}^{d} = 0 \quad \text{on} \quad \partial \Omega_F = \left\{ \xi \in \mathbb{R}^d : \sum_{m=1}^{d} \alpha_m \xi_m^2 = \beta \right\}. $$

This shows that (17) is satisfied, and a posteriori that $\mathcal{P}$ commutes with $\mathcal{K}$.

<figure><img src="image_5.png" /></figure>

We can generalize the previous result by considering domains of the form

$$ D(\mathbf{M}, v, c) := \{x \in \mathbb{R}^d : x^T \mathbf{M} x + v^T x + c \le 0\}. $$

**Theorem 1.** Let $\Omega_S = D(\mathbf{M}_S, v_S, c_S)$ and $\Omega_F = D(\mathbf{M}_F, v_F, c_F)$, for some symmetric, diagonalizable and invertibles matrices $\mathbf{M}_S, \mathbf{M}_F \in \mathbb{R}^{d \times d}$, vectors $v_S, v_F \in \mathbb{R}^d$ and scalars $c_S, c_F \in \mathbb{R}$. Let $\mathcal{K}$ the concentration operator associated to masks $m_S = \mathbf{1}_{\Omega_S}$ and $\widehat{m}_F = \mathbf{1}_{\Omega_F}$. Then, there exists a second-order differential operator $\mathcal{P}$ that commutes with $\mathcal{K}$.
Let $\mathbf{M}_S = U_S \Lambda_S U_S^T$ and $\mathbf{M}_F = U_F \Lambda_F U_F^T$, where $U_S, U_F \in \mathbb{R}^d$ are orthogonal matrices and $\Lambda_S, \Lambda_F$ are diagonal, and write $a_n = (\Lambda_S)_{n,n}$, $\alpha_n = (\Lambda_F)_{n,n}$. Let $w_S = -\frac{1}{2}\Lambda_S^{-1}U_S^T v_S$, $w_F = -\frac{1}{2}\Lambda_F^{-1}U_F^T v_F$, and define $b = w_S^T \Lambda_S w_S - c_S$, $\beta = w_F^T \Lambda_F w_F - c_F$. Then,

$$ \mathcal{P}(x, \nabla_x) = \nabla_x^T U^T \mathbf{A} \left( U_S^T x + \frac{1}{2} \Lambda_S^{-1} U_S^T v_S \right) U^T \nabla_x + C \left( U_S^T x + \frac{1}{2} \Lambda_S^{-1} U_S^T v_S \right). $$

where

$$ \mathbf{A}(y) := \operatorname{diag} \left\{ \alpha_m \left( \sum_{n=1}^{d} a_n y_n^2 - b \right) \right\}_{m=1}^{d} \quad \text{and} \quad C(y) := \beta y \cdot \operatorname{diag} \{a_n\}_{n=1}^{d} y. $$

*Proof.* We simply show that it is possible via an affine change of variables to recover the case of quadratic domains of Proposition 2.

Let us drop the indices $S$ and $F$ since we do the exact same computations. We assume $\mathbf{M}$ to be diagonalizable and symmetric, so there exist matrices $U, \Lambda \in \mathbb{R}^{d \times d}$ such that $\mathbf{M} = U\Lambda U^T$, where $U$ is an orthogonal matrix and $\Lambda$ is diagonal. We then have

$$ \Omega = \{x \in \mathbb{R}^d : (U^T x)^T \Lambda (U^T x) + v^T x + c \le 0\}. $$

This hints for the change of variable $\tilde{y} = U^T x$, and $\tilde{w} := U^T v$. Then,

$$ \Omega = \{ \tilde{y} \in \mathbb{R}^d : \tilde{y}^T \Lambda \tilde{y} + \tilde{w}^T \tilde{y} + c \le 0 \}. $$

Since $\mathbf{M}$ is invertible, so is $\Lambda$. Define $w = -\frac{1}{2}\Lambda^{-1}\tilde{w}$ and $y = \tilde{y} - w$, then

$$ \begin{aligned} y^T \Lambda y &= (\tilde{y} - w)^T \Lambda (\tilde{y} - w) = \tilde{y}^T \Lambda \tilde{y} - 2w^T \Lambda \tilde{y} + w^T \Lambda w \\ &= \tilde{y}^T \Lambda \tilde{y} + \tilde{w}^T \tilde{y} + w^T \Lambda w, \end{aligned} $$

and we get

$$ \Omega = \{y \in \mathbb{R}^d : y^T \Lambda y - w^T \Lambda w + c \le 0\}. $$

We observe that $\Omega$ is of the form (11).

By doing similar changes of variables for the space and Fourier domains $\Omega_S$ and $\Omega_F$, we can apply Proposition 2 and obtain a commuting differential operator $\mathcal{P}$ in variables $(y, \nabla_y)$.

In order to obtain the explicit expression of $\mathcal{P}$ in variables $(x, \nabla_x)$, we have to define a few quantities. Let $a_n := (\Lambda_S)_{n,n}$, $\alpha_n := (\Lambda_F)_{n,n}$, $b := w_S^T \Lambda_S w_S - c_S$, $\beta := w_F^T \Lambda_F w_F - c_F$, and recall $y = U_S^T x + \frac{1}{2}\Lambda_S^{-1}U_S^T v_S$. Also, define

$$ \mathbf{A}(y) := \operatorname{diag} \left\{ \alpha_m \left( \sum_{n=1}^{d} a_n y_n^2 - b \right) \right\}_{m=1}^{d} \quad \text{and} \quad C(y) := \beta y \cdot \operatorname{diag} \{a_n\}_{n=1}^{d} y. $$

Then

$$ (18) \qquad \mathcal{P}(y, \nabla_y) = \nabla_y \cdot (\mathbf{A}(y)\nabla_y) + C(y), $$

We have $\nabla_y = U^T \nabla_x$, thus

$$ \mathcal{P}(x, \nabla_x) = \nabla_x^T U^T \mathbf{A} \left( U_S^T x + \frac{1}{2} \Lambda_S^{-1} U_S^T v_S \right) U^T \nabla_x + C \left( U_S^T x + \frac{1}{2} \Lambda_S^{-1} U_S^T v_S \right). $$

<figure><img src="image_6.png" /></figure>

By Lemma 5, this commutation property allows us to look for eigenfunctions of $\mathcal{P}$ in order to know the eigenfunctions of $\mathcal{K}$.
**Lemma 5.** Let $d \ge 1$, and $Q, \mathcal{R}: L^2(\mathbb{R}^d) \to L^2(\mathbb{R}^d)$ two commuting operators acting on $L^2(\mathbb{R}^d)$. In other words, $QR = \mathcal{R}Q$. Suppose that each eigenfunction $\varphi_i$ of $Q$ is associated to an eigenvalue $\kappa_i$ of multiplicity one, and that $\{\varphi_i\}_{i \in \mathbb{N}}$ is a complete family in $L^2(\mathbb{R}^d)$. Then $Q$ and $\mathcal{R}$ have the same eigenfunctions.

*Proof.* Using the commutation relation between $Q$ and $\mathcal{R}$, one obtains

$$Q\varphi_i = \kappa_i\varphi_i \implies \mathcal{R}Q\varphi_i = \kappa_i\mathcal{R}\varphi_i \implies \mathcal{R}\mathcal{R}\varphi_i = \kappa_i\mathcal{R}\varphi_i.$$

This means that $\mathcal{R}\varphi_i$ is also an eigenfunction of $Q$, and since $\kappa_i$ is an eigenvalue of multiplicity one we must have $\mathcal{R}\varphi_i = c\varphi_i$ for some constant $c \in \mathbb{C}$. In other words, all eigenfunctions of $Q$ are eigenfunctions of $\mathcal{R}$. Moreover, since the eigenfunctions of $Q$ form a complete family of $L^2(\mathbb{R}^d)$, we deduce that the eigenfunctions of $\mathcal{R}$ are exactly the eigenfunctions of $Q$. $\square$

**3.2. A splitting approach.** Assume that $1 \ge m_S > 0$ and $1 \ge \widehat{m_F} > 0$. Then there exists functions $V(x) \ge 0$ and $H(\xi) \ge 0$ such that

$$m_S(x) = e^{-\frac{1}{2}V(x)} \quad \text{and} \quad \widehat{m_F}(\xi) = e^{-\frac{1}{2}H(\xi)}.$$

Owing to the fact that $\xi \widehat{f}(\xi) = -i \partial_x f(\xi)$, the operator $\mathcal{K}$ can be expressed under the splitting form

$$(19) \qquad \mathcal{K} = e^{-\frac{1}{2}V(x)}e^{-H(-i\partial_x)}e^{-\frac{1}{2}V(x)}$$

and we recognize a *Strang splitting* decomposition of the pseudo-differential operator $H(-i\partial_x) + V(x)$ (see [10, 11]). In particular, we can write formally the Baker-Campbell-Hausdorff (BCH) formula [3, 10]

$$e^{-\frac{1}{2}V(x)}e^{-H(-i\partial_x)}e^{-\frac{1}{2}V(x)} = e^{-Z(x, -i\partial_x)}$$

where

$$(20) \qquad Z(x, -i\partial_x) = H + V + \frac{1}{12}[H, [H, V]] - \frac{1}{24}[V, [V, H]] + \sum_{k \ge 2} Z_{2k+1}$$

where the $Z_{2k+1}$ are made of nested commutators between the operators $H$ and $V$. Note that in general, if $H$ and $V$ are polynomials, then the operators $Z_{2k+1}$ can be expressed as polynomials of higher degrees making the previous series non convergent and a general singularly perturbed problem. Before giving more precise example, we first give a case of convergence:

**Theorem 2.** Let $\alpha, \beta > 0$, $m_S = e^{-\frac{\alpha}{2}x^2}$ and $\widehat{m_F} = e^{-\frac{\beta}{2}\xi^2}$. Then we have

$$(21) \qquad \mathcal{K} = e^{-\frac{1}{2}\alpha x^2} e^{\beta \Delta} e^{-\frac{1}{2}\alpha x^2} = e^{-\operatorname{argsh}(\sqrt{\alpha\beta})} \left[ \sqrt{\frac{\alpha(1+\alpha\beta)}{\beta}} x^2 - \sqrt{\frac{\beta}{\alpha(1+\alpha\beta)}} \Delta \right].$$

The eigenpairs $(\lambda_n, \psi_n)_n$ of the operator $\mathcal{K}$ are given by

$$\psi_n = \left( \frac{\alpha(1 + \alpha\beta)}{\beta} \right)^{\frac{1}{8}} \varphi_n \left( \left( \frac{\alpha(1 + \alpha\beta)}{\beta} \right)^{\frac{1}{4}} x \right) \\ \lambda_n = e^{-\operatorname{argsh}(\sqrt{\alpha\beta})(2n+1)}$$

where $\varphi_n$ are Hermite functions.
Proof. We start with the following formula (see [1]), for $z > 0$,

$$e^{-\frac{1}{2} \tanh(z) x^2} e^{\frac{1}{2} \sinh(2z) \Delta} e^{-\frac{1}{2} \tanh(z) x^2} = e^{-z(x^2 - \Delta)},$$

from which we deduce by a scaling by $\sqrt{\cosh z}$ in $x$:

$$e^{-\frac{1}{2} \sinh(z) x^2} e^{\sinh(z) \Delta} e^{-\frac{1}{2} \sinh(z) x^2} = e^{-z \left( \cosh(z) x^2 - \frac{1}{\cosh(z)} \Delta \right)}$$

and hence for all $c > 0$

$$e^{-\frac{1}{2}cx^2} e^{c\Delta} e^{-\frac{1}{2}cx^2} = e^{-\operatorname{argsh}(c) \left[ \sqrt{1+c^2}x^2 - \frac{1}{\sqrt{1+c^2}} \Delta \right]}.$$

Now by scaling again we obtain

$$e^{-\frac{1}{2}c\lambda^2 x^2} e^{\frac{c}{\lambda^2}\Delta} e^{-\frac{1}{2}c\lambda^2 x^2} = e^{-\operatorname{argsh}(c) \left[ \sqrt{1+c^2}\lambda^2 x^2 - \frac{1}{\lambda^2\sqrt{1+c^2}} \Delta \right]}.$$

and by taking $\lambda^2 = \sqrt{\frac{\alpha}{\beta}}$ and $c = \sqrt{\alpha\beta}$ we obtain (21). Remark that this formula implies that the operator

$$\mathcal{H} = -\frac{1}{\mu^2}\Delta + \mu^2 x^2, \quad \mu^2 = \sqrt{\frac{\alpha(1+\alpha\beta)}{\beta}}$$

commutes with $\mathcal{K}$. Moreover, the spectrum of the operator $\mathcal{H}$ is given by the normalized functions $\sqrt{\mu}\varphi_n(\mu x)$ with eigenvalues $2n + 1$, where the $\varphi_n$ are the normalized Hermite functions. This implies the result. $\square$

**Remark 1.** Using the framework of [1], this result can be easily extended to more general masks with quadratic function $H(\xi)$ and $V(x)$, and in higher dimension. The details of the exact formula are left to the reader.

We conclude this section by giving another consequence of Formula (19) which is a possible construction of quasimodes for the operator $\mathcal{K}$ for masks that are close to the identity. Again, we give only a simple example of application. Let us define the spaces

$$\mathcal{V}^s = \{f \in L^2(\mathbb{R}^2), \ |\langle x \rangle^s f \in L^2(\mathbb{R}^d) \text{ and } \langle \xi \rangle^s \hat{f} \in L^2(\mathbb{R}^2)\},$$

where $\langle y \rangle^2 = 1 + |y|^2$ for $y \in \mathbb{R}^d$.

**Theorem 3.** Assume that $m_S(x) = e^{-\frac{1}{2}\varepsilon V(x)}$ and $\widehat{m}_F(\xi) = e^{-\frac{1}{2}\varepsilon H(\xi)}$ with $\varepsilon > 0$, and $H(\xi) \ge 0$, $V(x) \ge 0$ two given smooth functions with polynomial growth i.e. there exists $C$ and $r, p \ge 0$ such that

$$|V(x)| \le C \langle x \rangle^r, \quad |H(\xi)| \le C \langle \xi \rangle^p, \quad \forall (x, \xi) \in \mathbb{R}^{2d}.$$

Assume that $(\omega, \psi) \in \mathbb{R} \times L^2(\mathbb{R}^d)$ is an eigenpair of the self adjoint operator

$$T = H(-i\partial_x) + V(x), \quad \text{i.e.} \quad T\psi = \omega\psi.$$

Assume moreover that for all $s \ge 0$, $\psi \in \mathcal{V}^s$, Then there exists $C$ and $\varepsilon_0$ such that for $\varepsilon \le \varepsilon_0$ we have

$$\|\mathcal{K}\psi - e^{-\varepsilon\omega}\psi\|_{L^2} \le C\varepsilon^3.$$
*Proof.* The proof is a consequence of the classical bounds for operator splitting as in [11]. Indeed, we have for all $\varphi$ smooth enough

$$ \|e^{-\frac{1}{2}\varepsilon V} e^{-\varepsilon H} e^{-\varepsilon \frac{1}{2}V} \varphi - e^{-\varepsilon T} \varphi\|_{L^2} \leq C\varepsilon^3 \|\varphi\|_{\mathcal{V}^{N(p,r)}} $$

where the exponent $N(p, r)$ depends on commutator bounds of the operators $[H, [H, V]]$ and $[V, [V, H]]$ which depend themselves of $p, r$ and the algebraic structure of the commutators (see formula (2.4) of [11]). Note that the fact that $H, V$ and thus $T$ are positive operators allow to define the semi group actions $e^{-tH}$, $e^{-tV}$ and $e^{-tT}$ in the Sobolev spaces $\mathcal{V}^s$. By applying the previous estimate to the function $\psi$ we obtain the result. $\square$

**Remark 2.** We could imagine that the previous theorem could be extended to the construction of quasi-mode of arbitrary order (i.e. with a precision $\mathcal{O}(\varepsilon^N)$ for all $N$), but this would require more elaborated mathematical techniques. Works in this direction can be found in [6] for a related analysis in the same “parabolic” situation, and in [7] for a more general treatment of the BCH formula for PDEs.

## 4. THE VARYING MASKS ALGORITHM

We describe now a new algorithm for computing the spectrum of $\mathcal{K}$, which gives very good and promising results in known unstable situations. Note however that no rigorous analysis of this numerical method is performed in this paper. Numerical examples will be given in the next section.

**4.1. Discretization of the generalized concentration operator.** For simplicity we assume that the space and Fourier masks are compactly supported in the intervals $[-1, 1]^d$ and $[-\pi, \pi]^d$ respectively. The space domain is discretized using the uniform midpoint rule, and the stepsize is the same for all dimensions. If we let $N$ the number of discretization points in each dimension, then the stepsize is given by $\Delta x = \frac{2}{N}$. Nodes on this grid write $x^{(k)} = -1 + (k + \frac{1}{2}) \Delta x$ for $k \in [0, N-1]^d$.

We are interested in discretizing the concentration operator (5), which can be written in this case

$$ (\mathcal{K}f)(x) = \int_{[-1,1]^d} f(y)k(y,x)dy, \quad x \in [-1,1]^d $$

where

$$ k(y, x) = m_S(y) \overline{m_S(x)} \mathcal{F}^{-1} [|\widehat{m_F}(\xi)|^2] (x - y). $$

For every node $x^{(j)}$ in the space grid, we obtain

$$ \int_{[-1,1]^d} f(y)k(y, x^{(j)})dy \approx (\Delta x)^d \sum_{k \in [0, N-1]^d} f(y^{(k)})k(y^{(k)}, x^{(j)}), $$

where the variable $y$ is assumed to be discretized exactly as the $x$ variable. It remains to compute $k(y^{(k)}, x^{(j)})$ by approximating

$$ \mathcal{F}^{-1} [|\widehat{m_F}(\xi)|^2] (x^{(j)} - y^{(k)}). $$

Note that the difference $z^{(j,k)} := x^{(j)} - y^{(k)}$ belongs to a uniform discretization centered around origin, with the same stepsize $\Delta x = \frac{2}{N}$ and with $2N - 1$ uniform points in each dimension.
We want to approximate the inverse Fourier transform by an inverse discrete Fourier transform, let us show that this approximation holds.

Conceptually, since $\widehat{m}_F$ is assumed to be compactly supported on $[-\pi, \pi]^d$, we can consider the restriction $\widehat{m}_F|_{[-\pi, \pi]^d}$ and extend it periodically to $\mathbb{R}^d$. Therefore, we are looking at the inverse Fourier transform of a periodic function, which is a discrete function. Moreover, since we are discretizing the inverse Fourier transform, the function $\widehat{m}_F$ is only evaluated on a grid over $[-\pi, \pi]^d$. For simplicity, we assume this grid to be uniform with $M$ points in each dimension, so the resulting function is not only discrete but also periodic with period $M$. The Fourier integral will be discretized using the midpoint quadrature rule, so the stepsize is given by $\Delta\xi = \frac{2\pi}{M}$, and a node of this grid discretization writes $\xi^{(l)} = -\pi + (l + \frac{1}{2})\Delta\xi$ for $l \in [0, M-1]^d$. For $u \in [0, M-1]^d$, we get

$$ 
\begin{aligned}
\mathcal{F}^{-1}[|\widehat{m}_F(\xi)|^2](u) &= \frac{1}{(2\pi)^d} \int_{[-\pi, \pi]^d} |\widehat{m}_F(\xi)|^2 e^{i\xi \cdot u} d\xi \\
&\approx \frac{(\Delta\xi)^d}{(2\pi)^d} \sum_{l \in [0, M-1]^d} |\widehat{m}_F(\xi^{(l)})|^2 e^{i\xi^{(l)} \cdot u} d\xi \\
&\approx \frac{e^{i(-\pi + \frac{\Delta\xi}{2}) \mathbf{1} \cdot u}}{M^d} \sum_{l \in [0, M-1]^d} |\widehat{m}_F(\xi^{(l)})|^2 e^{i\frac{2\pi}{M} l \cdot u} d\xi,
\end{aligned}
$$

where $\mathbf{1} = (1, \dots, 1)$.

We recognize here the inverse Discrete Fourier transform (IDFT) of the function $|\widehat{m}_F|^2$ evaluated on the $\xi$-grid. Note that we actually want to compute this IDFT on the grid $x-y$ which has $2N-1$ points, so we choose $M = 2N-1$. Moreover, $u$ is the multi-index of some node in the $x-y$ grid, which means that if we are at the node $x^{(j)} - y^{(k)}$ for some $j, k \in [0, N-1]^d$, then $u = j-k$.

We obtain

$$ (Kf)(x^{(j)}) \approx (\Delta x)^d \sum_{k \in [0, N-1]^d} f(y^{(k)}) m_S(y^{(k)}) \overline{m_S(x^{(j)})} \frac{e^{i(-\pi + \frac{\Delta\xi}{2}) \mathbf{1} \cdot (j-k)}}{(2N-1)^d} \times \sum_{l \in [0, 2N-2]^d} |\widehat{m}_F(\xi^{(l)})|^2 e^{i\frac{2\pi}{2N-1} l \cdot (j-k)}. $$

For numerical reasons, it is advisable not to manipulate quantities that are too small. Therefore, we consider a normalized eigenproblem, normalized by $(\Delta x)^d$. Due to the scaling properties of the operator, this normalization is equivalent to work on a cube of order $(\Delta x)^{-d}[-1, 1]^d$ but does not change the original problem. To simplify further, we can multiply the $j$-th component of $v$ by $e^{i(-\pi + \frac{\Delta\xi}{2})\mathbf{1} \cdot j}$, where $v$ is the grid discretization of a function $f$. We then obtain the following eigenproblem:

$$ \lambda v = K v, $$

where the matrix $K$ is the concentration matrix, and it is given by

$$ (22) \qquad \mathbf{K}_{j,k} = \frac{m_S(y^{(k)}) \overline{m_S(x^{(j)})}}{(2N-1)^d} \sum_{l \in [0, 2N-2]^d} |\widehat{m}_F(\xi^{(l)})|^2 e^{i \frac{2\pi}{2N-1} l \cdot (j-k)}, $$
for $j, k \in [0, N-1]^d$. We recognize the Discrete Fourier transform, and we emphasize that the matrix defined by (22) allows to compute the matrix $\mathbf{K}$ efficiently using the Fast Fourier Transform. Moreover, $\mathbf{K}$ can be written as a block matrix. For simplicity, we give details below in the two-dimensional case:

$$ (23) \qquad \mathbf{K} = \begin{pmatrix} \mathbf{K}^{(0,0)} & \mathbf{K}^{(0,1)} & \dots & \mathbf{K}^{(0,N-2)} & \mathbf{K}^{(0,N-1)} \\ \mathbf{K}^{(1,0)} & \mathbf{K}^{(1,1)} & \dots & \mathbf{K}^{(1,N-2)} & \mathbf{K}^{(1,N-1)} \\ \vdots & \vdots & \dots & \vdots & \vdots \\ \mathbf{K}^{(N-2,0)} & \mathbf{K}^{(N-2,1)} & \dots & \mathbf{K}^{(N-2,N-2)} & \mathbf{K}^{(N-2,N-1)} \\ \mathbf{K}^{(N-1,0)} & \mathbf{K}^{(N-1,1)} & \dots & \mathbf{K}^{(N-1,N-2)} & \mathbf{K}^{(N-1,N-1)} \end{pmatrix}, $$

where each block $\mathbf{K}^{(r,c)}$ is defined component-wise by

$$ [\mathbf{K}^{(r,c)}]_{m,n} := \mathbf{K}_{j=(r,m), k=(c,n)}. $$

Let us now explain why this indexing convention is particularly efficient. In (22), we can write $j-k = (j_1-k_1, j_2-k_2) = (r-c, m-n)$. This means that, along each diagonal of $\mathbf{K}^{(r,c)}$, the value of $(j-k)$ is constant: indeed, $r-c$ is constant in $\mathbf{K}^{(r,c)}$, and $m-n$ is constant along each diagonal of $\mathbf{K}^{(r,c)}$. In other words, each submatrix $\mathbf{K}^{(r,c)}$ is a Toeplitz matrix, multiplied row-wise by the function $\overline{m_S}$ and column-wise by $m_S$. The Toeplitz nature of each block $\mathbf{K}^{(r,c)}$ allows for efficient computational storage and complexity.

This two-dimensional discussion easily generalizes to the multi-dimensional case when only the last index varies within each block $\mathbf{K}^{(r_1, \dots, r_{d-1})}$ of $\mathbf{K}$. Therefore, each block can be expressed as a Toeplitz matrix and component-wise multiplications.

We have the following easy result, which is more or less a discrete version of Proposition 1:

**Proposition 3.** The matrix $\mathbf{K}$ enjoys the following properties:

a) Hermitian character: $\mathbf{K}^* = \mathbf{K}$.

b) Structure: $\mathbf{K} = D^*BD$, with $B$ a block matrix where each block is Toeplitz, and $D$ is a diagonal matrix.

c) Its eigenvalues are real.

d) Its eigenvectors form a unitary basis of $\mathbb{C}^{N_1 \cdots N_d}$.

*Proof.* For the first point, we use (22) combined with the fact that $|\overline{m_F}|$ is a real-valued function.

The second point has already been mentioned earlier: the diagonal matrix $D$ follows from (22), and it corresponds to the component-wise multiplication by the function $\overline{m_S}(x^{(j)})$ for each row $j$ and by $m_S(y^{(k)})$ for each column $k$ of $\mathbf{K}$. Hence $\mathbf{K}$ is of the form $D^*BD$, where $B$ is some matrix. The block nature of $B$, where each block is a Toeplitz matrix, follows from (23).

The third point is due to the Hermitian character of $\mathbf{K}$. The fourth point is a classical result in linear algebra: for any normal matrix, there exists an orthonormal basis of eigenvectors. See for instance [2, Theorem 7.31]. $\square$

We end this section with a remark concerning notation: we are interested in eigenpairs of the matrix $\mathbf{K}$ of finite dimension. With a slight abuse of notation, we will denote $\psi_i$ the eigenvectors, which is the same notation as used for the eigenfunctions of $\mathcal{K}$. Whether we are talking about an eigenvector or an eigenfunction will always be clear from the context: if we are talking about the continuous concentration operator $\mathcal{K}$, $\psi_i$ will denote an eigenfunction,
and if we are talking about the discretized version of $\mathcal{K}$ (i.e. the matrix $\mathbf{K}$), then $\psi_i$ will denote an eigenvector.

4.2. **Approximating eigenvectors.** We focus now on the case of binary compact filters, i.e.

$$m_S = \mathbf{1}_{\Omega_S} \quad \text{and} \quad \widehat{m_F} = \mathbf{1}_{\Omega_F},$$

for two compact subsets $\Omega_S, \Omega_F \subset \mathbb{R}^d$.

Let $\Omega_S : \mathbb{R}_+ \to [0, 1]$ and $\Omega_F : \mathbb{R}_+ \to [0, 1]$ two set-valued functions that depend on a parameter $\varepsilon$. Note the slight abuse of notation, where $\Omega_i(\cdot)$ denotes a set-valued function while $\Omega_i$ denotes a subset of $\mathbb{R}^d$, $i \in \{S, F\}$. The functions $\Omega_i(\cdot)$ are chosen such that $\Omega_i(0) = \Omega_i$, and such that $\Omega_i(\varepsilon)$ reduces to a set $Z_i$ of null measure as $\varepsilon \to \infty$, $i \in \{S, F\}$. To these set-valued functions we associate modified masks $m_S(\varepsilon, \cdot) := \mathbf{1}_{\Omega_S(\varepsilon)}$ and $\widehat{m_F}(\varepsilon, \cdot) = \mathbf{1}_{\Omega_F(\varepsilon)}$. We may use the shorthands $m_S = m_S(0, \cdot)$ and $\widehat{m_F} = \widehat{m_F}(0, \cdot)$.

We denote $\mathcal{K}(\varepsilon)$ the spectral concentration operator with space mask $m_S(\varepsilon, \cdot)$ and Fourier mask $\widehat{m_F}(\varepsilon, \cdot)$, and $\mathbf{K}(\varepsilon)$ its discretization as described in Section 4.1. We may use the shorthand $\mathbf{K} = \mathbf{K}(0)$. Moreover, we write $(\lambda_n(\varepsilon), \psi_n(\varepsilon))$ the n-th eigenpair of the modified concentration matrix $\mathbf{K}(\varepsilon)$.

We define the concentration ratio of a vector similarly to the continuous setting:

$$\nu(v) = \frac{v^*\mathbf{K}v}{v^*v}.$$

Note that it is computed with respect to the initial concentration matrix $\mathbf{K}(0)$.

We now make an assumption, which will be crucial in the following. It has always been observed to hold in practice during our experiments, so we believe this assumption does not impose too much restriction.

**Assumption 1.** It is assumed that, as $\varepsilon \to 0$, the eigenvalues corresponding to the concentration operator with masks $m_S(\varepsilon, \cdot)$ in space and $\widehat{m_F}(\varepsilon, \cdot)$ in Fourier are such that the first eigenvalue $\lambda_1(\varepsilon)$ reaches $\lambda_1(0)$ before $\lambda_2(\varepsilon)$ reaches $\lambda_2(0)$, and so on for the next eigenvalues $\lambda_n(\varepsilon)$. In other words, we assume that eigenvalues corresponding to the modified masks are all (numerically!) distinct for $\varepsilon > 0$, and that the order $\lambda_i(\varepsilon) > \lambda_{i+1}(\varepsilon)$ is preserved for all $\varepsilon$.

To illustrate this assumption we go back to the one-dimensional historical example by Slepian, but instead of considering $m_S = \mathbf{1}_{[-1,1]}$ and $\widehat{m_F} = \mathbf{1}_{[-c,c]}$, we consider their modified versions $m_S(\varepsilon, \cdot) = \mathbf{1}_{[-\mu(\varepsilon),\mu(\varepsilon)]}$ and $\widehat{m_F}(\varepsilon, \cdot) = \mathbf{1}_{[-c\mu(\varepsilon),c\mu(\varepsilon)]}$, where

$$(24) \qquad \mu(\varepsilon) := \frac{1}{(1 + \varepsilon^4)^{1/4}}.$$

There is no numerical issue in looking for the eigenvalues, even if they are close to each other. Issues arise when we look for eigenvectors. We give in Figure 1 the 30 first eigenvalues obtained for several values of $\varepsilon$. In Figure 2, we give the first 16 eigenvectors obtained with $\varepsilon = 0$ (solid blue curve), with $\varepsilon = 2$ (dashed orange curve), as well as the exact ones, in the case of a concentration matrix corresponding to $m_S(\varepsilon, \cdot) = \mathbf{1}_{[-\mu(\varepsilon),\mu(\varepsilon)]}$ and $\widehat{m_F}(\varepsilon, \cdot) = \mathbf{1}_{[-0.1\cdot 2\pi\mu(\varepsilon),0.1\cdot 2\pi\mu(\varepsilon)]}$. The exact eigenvectors in this particular case are known as the eigenvectors of a tridiagonal matrix and they are called Discrete Prolate Spheroidal Sequences (DPSS), see e.g. [21].

We can observe that $\varepsilon = 0$ yields eigenvectors that are not symmetric (neither even nor odd), while $\varepsilon = 2$ does. It is known in this particular situation that the eigenvectors are
<figure><img src="image_7.png" /><figcaption>FIGURE 1. 30 first eigenvalues of the matrix $\mathbf{K}(\varepsilon)$, with $N = 150$ points of discretization, $\Omega = 0.1 \cdot 2\pi$.</figcaption></figure>

symmetric. Moreover, $\varepsilon = 2$ yields some vectors that may happen to be close to the desired ones: the first vectors are very close to the exact eigenvectors, while the next vectors are very different from the exact eigenvectors.

The key takeaway from this simple experiment is that when the concentration problem is scaled down (or “shrinked”), the first eigenvectors are the same but the next ones have a smaller concentration ratio. This means that it is now possible and easy to obtain approximate eigenvectors, by looking at different values of $\varepsilon > 0$.

This motivates the following idea for an approximate eigendecomposition of the spectral concentration matrix $\mathbf{K}$: start from very narrow modified masks (i.e. $\varepsilon$ large), then it is easy to obtain the first eigenvector $\psi_1(\varepsilon)$ of $\mathbf{K}(\varepsilon)$ since the gap between the first and second eigenvalues is large according to Assumption 1. However, this vector $\psi_1(\varepsilon)$ is probably not a satisfying approximation of the first eigenvector $\psi_1(0)$ of $\mathbf{K}(0)$, in the sense that its concentration ratio may not be close to the exact first eigenvalue $\lambda_1(0)$ of $\mathbf{K}(0)$. This is the situation occuring for $\varepsilon = 10$ or $\varepsilon \approx 5.6$ in Figure 1. When this happens, we can take $\varepsilon$ smaller and check again if $\psi_1(\varepsilon)$ is a satisfying approximation of $\psi_1(0)$. For some $\varepsilon_1$ small enough, the concentration ratio of $\psi_1(\varepsilon_1)$ becomes close enough (down to some prescribed tolerance $\eta$) to the first eigenvalue $\lambda_1(0)$ of $\mathbf{K}(0)$. In this case, we consider that $\psi_1(\varepsilon_1)$ is a good enough approximation of $\psi_1(0)$, and we record this vector $\psi_1(\varepsilon_1)$. We then repeat this process in order to find an approximation of the second eigenvector $\psi_2(0)$ of $\mathbf{K}(0)$. Note that, when looking for $\psi_2(\varepsilon)$ an approximation of $\psi_2(0)$, we can look for it in the orthogonal of $\text{Span}\{\psi_1(\varepsilon_1)\}$. This will guarantee that the set of vector we obtain at the end of the procedure is indeed a basis of $\mathbb{C}^N$.

The orthogonalization $\tilde{w}$ of a vector $w$ with respect to a vector $v$ is simply given by
$$ \tilde{w} = \frac{u}{\sqrt{u^*u}} \text{ where } u = (1 - \frac{vv^*}{v^*v})w. $$
In practice, we use an iterative process to approximate eigenvectors, so the orthogonalization is required at each iteration.

By construction of this basis of approximate eigenvectors, the concentration ratio of the $n$-th vector $\psi_n(\varepsilon_n)$ is $\eta$-close to the true $n$-th eigenvalue $\lambda_n(0)$. Another advantage of this procedure is that we do not need to know beforehand all the exact eigenvalues, since it is possible to work one eigenvector at a time.
<figure><img src="image_8.png" /><figcaption>FIGURE 2. First sixteen eigenvectors obtained for the initial discretized concentration problem ($\epsilon = 0$, solid blue), for the modified discretized concentration problem ($\epsilon = 2$, orange dashes), and the exact eigenvectors given by the DPSS (green dots).</figcaption></figure>
The procedure is described above in the one-dimensional case, but the same ideas apply in the multi-dimensional setting. Note that some care has to be taken when choosing the set-valued function $\Omega_i(\cdot)$, whose role is to mimic the scaling of the interval when $d = 1$. For example, if $\Omega_i$ has “holes”, then we don’t want the holes to move, only possibly to get bigger or smaller. This explains why, in the case $\Omega_S = \text{cat-head}$ given in Section 5.4, the set-valued function $\Omega_S(\cdot)$ cannot simply be a scaling of $\Omega_S$. Otherwise, the holes within the cat-head shape would move, and the convergence would require smaller values of $\varepsilon$.

The algorithm is given by Algorithm 1.

**Remark 3.** When looking for an eigenvector $\psi(\varepsilon)$ of $\mathbf{K}(\varepsilon)$, one can take advantage of the eigenvector obtained for a previous (larger) value of $\varepsilon$, in order to start the search of an approximate eigenvector and thus speed up computations.

**Remark 4.** In Algorithm 1, we allow considering approximate eigenvalues of the matrix $\mathbf{K}(0)$. This is due to the fact that, with high-dimensional matrices, it is often difficult to obtain the eigenvalues precisely because of the computational cost. Thus we allow approximate eigenvalues so that one can speed up the computations of $\lambda_n(0)$ by only looking for an $\eta$-close approximation.

Under our current assumptions, the fact that the procedure described above yields a correct approximation of the eigenvectors is shown by Lemma 6.

**Lemma 6.** Let $\mathbf{A} \in \mathbb{R}^{N \times N}$ an Hermitian matrix, and denote $\{\lambda_i\}_{i=1}^N$ its eigenvalues (they are all real), ordered so that $\lambda_i \ge \lambda_{i+1}$. Let $\{v_1, \dots, v_N\}$ an orthonormal basis of $\mathbb{C}^N$, where $v_i \in \mathbb{C}^N$ is an eigenvector of $\mathbf{A}$ associated to $\lambda_i$. Let $\eta > 0$ and $w \in \mathbb{C}^N$, such that $|w|_2 = \sqrt{w^*w} = 1$ and

$$ (25) \qquad |w^*\mathbf{A}w - \lambda_1| \le \eta. $$

For $m \le N$, we have

$$ |w - \text{Proj}_{\text{Span}\{v_1, \dots, v_m\}} w|_2^2 \le \frac{\eta}{\lambda_1 - \lambda_{m+1}}. $$

*Proof.* Decompose $w$ into the $\{v_i\}_{i=1}^N$ basis: $w = \sum_{i=1}^N c_i v_i$, for some coefficients $c_i \in \mathbb{C}$ such that $\sum_{i=1}^N |c_i|^2 = 1$. Owing to

$$ w^* \mathbf{A} w = \left( \sum_{i=1}^{N} c_i v_i \right)^* \left( \sum_{j=1}^{N} c_j \lambda_j v_j \right) = \sum_{i=1}^{N} |c_i|^2 \lambda_i, $$
**Algorithm 1** Varying masks method

**Require:**

* $\Omega_S, \Omega_F$: two finite-volume domains of $\mathbb{R}^d$.

* $\varepsilon \mapsto \Omega_i(\varepsilon)$: a set-valued function, decreasing for the relation of set inclusion, such that $\Omega_i(0) = \Omega_i$ and $\Omega_i(+\infty) = Z_i$ with meas$(Z_i) = 0$, $i \in \{S, F\}$.

* $m_S(\varepsilon, \cdot) = \mathbf{1}_{\Omega_S(\varepsilon)}$ and $\widehat{m}_F(\varepsilon, \cdot) = \mathbf{1}_{\Omega_F(\varepsilon)}$: the modified masks.

* $N$: the number of discretization points for each dimension.

* $\varepsilon \mapsto \mathbf{K}(\varepsilon)$: the modified concentration matrix. It is a square matrix of size $N^d$, corresponding to the modified masks $m_S(\varepsilon, \cdot)$ and $\widehat{m}_F(\varepsilon, \cdot)$.

* $M$: number of eigenvectors we are looking for, $M \le N$.

* $\varepsilon_{\max}, \varepsilon_{\min}$: maximum and minimum value of the parameter $\varepsilon$.

* $\{\varepsilon_T, \dots, \varepsilon_1\} \subset [\varepsilon_{\min}, \varepsilon_{\max}]$: $T$ discretization points of the interval $[\varepsilon_{\min}, \varepsilon_{\max}]$ (can be a uniform discretization, log-uniform, ...). They are assumed to be such that $\varepsilon_t > \varepsilon_{t-1}$, $t \in [1, T]$.

* $\eta$: numerical tolerance to compare two eigenvalues.

$q := 0$: this is the number of recorded eigenvectors yet.

$\alpha_{saved} = (\alpha_{saved,1}, \dots, \alpha_{saved,M})$: the vector to hold all the concentration ratios.

$v_{saved} = (v_{saved,1}, \dots, v_{saved,M})$: the matrix to hold all the recorded eigenvectors

**for** $\varepsilon = \varepsilon_T, \dots, \varepsilon_0$ **do**

Find $\tilde{\lambda}_q(0)$, an $\eta$-approximation of the exact eigenvalue $\lambda_q(0)$. A standard eigenalgorithm can be used here.

Find $(\kappa, u)$ the most significant eigenpair of $\mathbf{K}(\varepsilon)$ (i.e. the one associated to the eigenvalue of highest magnitude), where $u \perp \text{Span}\{v_{saved,1}, \dots, v_{saved,q}\}$ and $|u|_2 = 1$. A standard eigenalgorithm can be used here.

Compute the concentration ratio with respect to the unperturbed problem:

$$ \nu := \frac{u^*\mathbf{K}(0)u}{u^*u}. $$

if $|\nu - \tilde{\lambda}_q(0)| \le \eta$ then
$v_{saved,q} \leftarrow u$.
$\alpha_{saved,q} \leftarrow \nu$.
$q \leftarrow q + 1$.

end if

Stop if $q \ge M$.

end for

we get

$$
\begin{align*}
\lambda_1 - w^* \mathbf{A} w &= \sum_{i=1}^{N} |c_i|^2 (\lambda_1 - \lambda_i) = \sum_{i=1}^{m} |c_i|^2 (\lambda_1 - \lambda_i) + \sum_{i=m+1}^{N} |c_i|^2 (\lambda_1 - \lambda_i) \\
&\geq \sum_{i=1}^{m} |c_i|^2 (\lambda_1 - \lambda_i) + \sum_{i=m+1}^{N} |c_i|^2 (\lambda_1 - \lambda_{m+1}) \\
&= \sum_{i=1}^{m} |c_i|^2 (\lambda_1 - \lambda_i) + \left(1 - \sum_{i=1}^{m} |c_i|^2\right) (\lambda_1 - \lambda_{m+1}) \\
&= \sum_{i=1}^{m} |c_i|^2 (\lambda_{m+1} - \lambda_i) + \lambda_1 - \lambda_{m+1}.
\end{align*}
$$
Use the ordering of eigenvalues to obtain

$$ \eta \geq \lambda_1 - w^*Aw > (\lambda_1 - \lambda_{m+1}) \sum_{i=m+1}^{N} |c_i|^2 $$

Therefore,

$$ |w - \text{Proj}_{\text{Span}\{v_1, \dots, v_m\}} w|_2^2 = \left| \sum_{i=m+1}^{N} c_i v_i \right|_2^2 = \sum_{i=m+1}^{N} |c_i|^2 \leq \frac{\eta}{\lambda_1 - \lambda_{m+1}}. $$

□

**Corollary 1.** Under the same assumptions as Lemma 6, we have

$$ |Aw - \lambda_1 w|_2 \leq \frac{\lambda_1 - \lambda_n}{\lambda_1 - \lambda_2} \eta. $$

*Proof.* By Lemma 6, we can write $w = r + c_1v_1$, where $c_i \in \mathbb{C}$ and $r \in \mathbb{C}^d$ such that $|r|_2 \leq \frac{\eta}{\lambda_1 - \lambda_2}$. Thus,

$$ Aw = Ar + \lambda_1 c_1 v_1 $$

and

$$ Aw - \lambda_1 w = Ar - \lambda_1 r. $$

Finally,

$$ |Aw - \lambda w|_2 \leq |||A - \lambda_1 I|||_2 |r|_2 \leq \frac{\lambda_1 - \lambda_n}{\lambda_1 - \lambda_2} \eta, $$

where we have used the estimate of $|r|_2$.

□

**Remark 5.** Lemma 6 is useful numerically: if one considers the criterion (25) to select the i-th eigenvector $v_i$, then the numerical tolerance $\eta$ has to be small compared the difference $\lambda_i - \lambda_{i+1}$ in order to have a good approximation of $v_i$. In particular, if (25) is satisfied with $\eta \ll \lambda_1 - \lambda_2$, then we have a good approximation of $v_1$.

**Remark 6 (Continuation methods).** Since we are considering a perturbed problem for $\varepsilon > 0$, it seems natural to have in mind continuation methods. Here, the clusters of eigenvalues at zero and one force us to consider continuation methods where eigenvalues have to be distinguished using their derivatives because their value is not enough. This type of method has for instance been introduced and used in [12, 18]. However, unreported numerical experiments showed that the approximate eigenvectors obtained are not significantly different from those presented here with the varying masks procedure. Moreover, continuation methods have some drawbacks that are naturally fixed using the varying masks procedure:

* we cannot control how close the concentration ratios are to the true eigenvalues;

* in order to compute the eigenvector derivative, all eigenvectors are needed, which is very costly in two- and higher-dimensional settings;

* continuation methods require a concentration matrix that is differentiable with respect to the perturbation parameter $\varepsilon$, but the indicator masks considered are not differentiable. Thus, applying continuation methods would only help us solve an approximate concentration problem where the masks would be smooth (with respect to $\varepsilon$);

* the orthogonalization of eigenvectors is not guaranteed with eigenvector continuation (though some work could probably be done in that regard to alleviate this).
## 5. NUMERICAL EXAMPLES

In this section we will compare the eigenvectors obtained from a standard eigendecomposition of the concentration matrix $K(0)$, with the approximate eigenvectors obtained via the varying masks procedure described in Algorithm 1.

Note that the eigenvectors are a priori complex vectors, but the Fourier restrictions we consider in the following numerical examples are chosen even so that the eigenvectors can actually be chosen real. This is why we only present the real part of the (approximate) eigenvectors and don't mention their imaginary part.

We present in this section the results obtained using Algorithm 1, where both space and Fourier masks are varying, with the same variation $\mu(\varepsilon)$. During our numerical experiments we have also tried to let only one of the two masks vary. However, the results were always worse than with the simultaneously varying masks. We emphasize that Algorithm 1 is not rigorously justified, and we propose it only because it yielded interesting results that were better than a standard eigenalgorithm. For instance, we don't know if the choice of $\mu(\varepsilon)$ is the best, nor if it is better to have different variations $\mu_S$ and $\mu_F$ for the space and Fourier masks.

### 5.1. Interval masks – Moderate Fourier restriction.
This is a one-dimensional example, and we choose the number of discretization points $N = 150$. This means that the matrix is rather small, so we can have a very small numerical tolerance, hence we choose $\eta = 10^{-10}$.

The $\varepsilon$-discretization is a log-uniform discretization of the interval $[10^{-1}, 10^2]$, with $T = 250$ discretization points. The space mask is $m_S(\varepsilon, \cdot) = \mathbf{1}_{[-\mu(\varepsilon), \mu(\varepsilon)]}$ where $\mu(\varepsilon)$ is given by (24), and the Fourier mask is $\widehat{m}_F(\varepsilon, \cdot) = \mathbf{1}_{[-\omega\mu(\varepsilon), \omega\mu(\varepsilon)]}$ where $\omega = 0.3 \cdot 2\pi$.

We plot in Figure 3 the eigenvectors obtained with a standard eigendecomposition of $K(0)$ (solid blue curve) as well as the exact eigenvectors obtained via $\mathbf{P}$ (dashed orange curve). We can note that the eigendecomposition yields eigenvectors that are not localized, and they are also not symmetric (neither even nor odd).

The approximate eigenvectors obtained via the varying masks procedure are given in Figure 4. They are much closer to the exact eigenvectors, and exhibit the expected localization and symmetry properties. Note that we have not postprocessed the approximate eigenvectors, and in particular we have not applied any normalization convention. This explains why some approximate eigenvectors are the opposite of the desired ones. Moreover, the concentration ratios of the approximate eigenvectors are guaranteed to be $\eta$-close to the true eigenvalues, so the varying masks procedure is a satisfying robust alternative to the eigendecomposition.

### 5.2. Interval masks – Strong Fourier restriction.
For this second example, all the parameters of the concentration problem are the same as in Section 5.1, except for the Fourier restriction which is now $\widehat{m}_F = \mathbf{1}_{[-\omega\mu(\varepsilon), \omega\mu(\varepsilon)]}$ with $\omega = 0.49 \cdot 2\pi$. This is an important example, since almost all eigenvalues are clustered at one.

Figure 5 shows the eigendecomposition results, and Figure 6 shows the results obtained with the varying masks procedure. Once again, the varying masks procedure allows to recover approximate eigenvectors with a concentration ratio $\eta$-close to the exact eigenvalues of $K(0)$, and these vectors exhibit localization and symmetry properties that are lacking from the
<figure><img src="image_9.png" /><figcaption>FIGURE 3. 16 first eigenvectors obtained with a standard eigendecomposition (solid blue curve) of $K(0)$, compared to the exact eigenvectors (dashed orange curve). Here, $\hat{m}_F = \mathbf{1}_{[-0.3 \cdot 2\pi\mu(\epsilon), 0.3 \cdot 2\pi\mu(\epsilon)]}$. $N_1 = 150$, $\eta = 10^{-10}$.</figcaption></figure>

eigenvectors obtained via a standard eigendecomposition. They are, however, slightly more localized than expected.

5.3. **Two-dimensional – Centered balls.** We now consider the classical two-dimensional example of centered balls. Since the matrix $K$ is a square matrix of size $N^2$, we have to keep this product reasonable for computational purposes. We recall that $N$ is the number of discretization points in each dimension. We choose $N = 60$, and $\eta = 10^{-6}$.
<figure><img src="image_10.png" /><figcaption>FIGURE 4. 16 first eigenvectors obtained with the varying masks procedure (solid blue curve) of $K(0)$, compared to the exact eigenvectors (dashed orange curve). Here, $\widehat{m}_F = \mathbf{1}_{[-0.3 \cdot 2\pi\mu(\varepsilon), 0.3 \cdot 2\pi\mu(\varepsilon)]}$. $N_1 = 150$, $\eta = 10^{-10}$.</figcaption></figure>

The $\varepsilon$-discretization is a log-uniform discretization of the interval $[10^{-1}, 10^1]$, with $T = 250$ discretization points. We choose $m_S(\varepsilon, \cdot) = \mathbf{1}_{B(0,0.8\mu(\varepsilon))}$, where $\mu(\varepsilon)$ is again given by (24), and $\widehat{m}_F(\varepsilon, \cdot) = \mathbf{1}_{B(0,0.3 \cdot 2\pi\mu(\varepsilon))}$.

The results obtained using a standard eigendecomposition are given in Figure 7, and the Fourier transform of each eigenvector is given in Figure 8. We can note that the eigenvectors indeed have a compact support in space and that they are concentrated in the Fourier domain,
<figure><img src="image_11.png" /><figcaption>FIGURE 5. 16 first eigenvectors obtained with a standard eigendecomposition (solid blue curve) of $K(0)$, compared to the exact eigenvectors (dashed orange curve). Here, $\hat{m}_F = \mathbf{1}_{[-0.49 \cdot 2\pi\mu(\varepsilon), 0.49 \cdot 2\pi\mu(\varepsilon)]}$. $N_1 = 150$, $\eta = 10^{-10}$.</figcaption></figure>

but they do not exhibit the expected localization and symmetry properties. For instance, the first eigenvector is not more localized than the sixteenth

The results obtained using the varying masks procedure are given in Figure 9, and the Fourier transform of each vector is given in Figure 10. These approximate eigenvectors do exhibit the expected localization and symmetry properties, illustrating once again that the varying masks procedure is a robust approximate alternative to the eigendecomposition of the initial concentration matrix.
<figure><img src="image_12.png" /><figcaption>FIGURE 6. 16 first eigenvectors obtained with the varying masks procedure (solid blue curve) of $K(0)$, compared to the exact eigenvectors (dashed orange curve). Here, $\hat{m}_F = \mathbf{1}_{[-0.49 \cdot 2\pi\mu(\varepsilon), 0.49 \cdot 2\pi\mu(\varepsilon)]}$. $N_1 = 150$, $\eta = 10^{-10}$.</figcaption></figure>

5.4. **Two-dimensional – cat-head shape.** This example is studied because of its sharp corners, for having a unique symmetry, and for having holes. The set-valued function $\Omega_S(\varepsilon)$ is shown in Figure 11. Note that it is not simply a scaling of $\Omega_S = \Omega_S(0)$, because we have to take care of the holes in the domain. Indeed, we want these holes to not move, only to possibly grow or shrink. This is to guarantee that an eigenvector $v(\varepsilon)$ of $K(\varepsilon)$ with eigenvalue $\lambda(\varepsilon)$ will have a concentration ratio $\nu(v(\varepsilon))$ with respect to $K(0)$ that is larger than $\lambda(\varepsilon)$.
<figure><img src="image_13.png" /><figcaption>FIGURE 7. 16 first eigenvectors obtained with a standard eigendecomposition of $\mathbf{K}(0)$, where $m_S(\varepsilon, \cdot) = \mathbf{1}_{B(0,0.8\mu(\varepsilon))}$ and $\tilde{m}_F(\varepsilon, \cdot) = \mathbf{1}_{B(0,0.3\cdot2\pi\mu(\varepsilon))}$. Here, $m_S$ is outlined in gray. $N = 60$, $\eta = 10^{-6}$.</figcaption></figure>

The results obtained via a standard eigendecomposition are given in Figure 12, and the Fourier transform of each eigenvector is given in Figure 13. The results obtained via the varying masks procedure are given in Figure 14, and the Fourier transform of each approximate eigenvector is given in Figure 15.

We can draw the same conclusions as those of Section 5.3: the eigendecomposition results in eigenvectors that are not localized as expected, and they also do not have the expected symmetry. The varying masks procedure yields approximate eigenvectors with a concentration
<figure><img src="image_14.png" /><figcaption>FIGURE 8. Absolute value of the Fourier transform of the 16 first eigenvectors obtained with a standard eigendecomposition of $\mathbf{K}(0)$, where $m_S(\varepsilon, \cdot) = \mathbf{1}_{B(0,0.8\mu(\varepsilon))}$ and $\widehat{m_F}(\varepsilon, \cdot) = \mathbf{1}_{B(0,0.3 \cdot 2\pi\mu(\varepsilon))}$. Here, $\widehat{m_F}$ is outlined in gray. $N = 60, \eta = 10^{-6}$.</figcaption></figure>

ratio that is 2$\eta$-close to the exact eigenvalues (or $\eta$-close if we consider $\lambda_q$ instead of $\tilde{\lambda}_q$ in Algorithm 1), they exhibit the expected symmetry, and they also have the expected localization properties.

## 6. PERSPECTIVES AND UNANSWERED INTERROGATIONS

In this work we have extended the framework laid by Slepian, Pollak, and Landau in the 1960s and 1970s, when they considered the problem of finding $L^2$ functions with the best
<figure><img src="image_15.png" /><figcaption>FIGURE 9. 16 first eigenvectors obtained with the varying masks procedure of $\mathbf{K}(0)$, where $m_S(\varepsilon, \cdot) = \mathbf{1}_{B(0,0.8\mu(\varepsilon))}$ and $\widehat{m_F}(\varepsilon, \cdot) = \mathbf{1}_{B(0,0.3\cdot2\pi\mu(\varepsilon))}$. Here, $m_S$ is outlined in gray. $N = 60$, $\eta = 10^{-6}$.</figcaption></figure>

simultaneous space and Fourier localization. They solved the problem in a very elegant manner by looking at particular situations, but their ideas do not seem to be applicable in most cases. We extended their work from balls in $\mathbb{R}^d$ to general quadric domains in $\mathbb{R}^d$, and considered a more general framework inspired by the natural *splitting representation* of the concentration operator (19) which allows to obtain exact formulas in the gaussian case and quasi-modes for masks close to the identity.

We then introduced a new numerical algorithm for computing the eigenmodes of the concentration operator, named *varying masks method* which retains the idea of scaling
<figure><img src="image_16.png" /><figcaption>FIGURE 10. Absolute value of the Fourier transform of the 16 first eigenvectors obtained with the varying masks procedure of $K(0)$, where $m_S(\varepsilon, \cdot) = \mathbf{1}_{B(0,0.8\mu(\varepsilon))}$ and $\widehat{m}_F(\varepsilon, \cdot) = \mathbf{1}_{B(0,0.3\cdot2\pi\mu(\varepsilon))}$. Here, $\widehat{m}_F$ is outlined in gray. $N = 60$, $\eta = 10^{-6}$.</figcaption></figure>

progressively the masks from a situation where no relevant mode can be found (when both masks are too concentrated) to the targeted final situation. This simultaneous modification of the space and Fourier masks can be easily implemented in practice, and we reported excellent behavior of this algorithm. In particular, we showed on several numerical examples that the given procedure seems to be more robust than a standard eigendecomposition algorithm, and the results using this procedure exhibit better localization and symmetry properties than standard eigendecomposition algorithms.
<figure><img src="image_17.png" /><figcaption>FIGURE 11. Set-valued function $\Omega_S(\varepsilon)$ when $\Omega_S = \Omega_S(0) = \text{cat-head}$, for $\varepsilon \in \{0, 1, 3/2\}$.</figcaption></figure>

However, the main idea of this work is based upon Assumption 1, and some detailed analysis is required to determine if it actually holds. Some analysis is also required to give the proposed algorithm some rigorous foundations, as well as study how different scalings in space and Fourier would affect the results. Finally, we believe that the splitting approach (19) can be useful for more general understanding of the spectrum of the spectral concentration operator in very general situation.

## REFERENCES

[1] Paul Alphonse and Joackim Bernier, *Polar decomposition of semigroups generated by non-selfadjoint quadratic differential operators and regularizing effects*, Annales scientifiques de l'École Normale Supérieure **56** (2023), no. 2, 323-382.

[2] Sheldon Axler, *Linear algebra done right*, Undergraduate Texts in Mathematics, Springer International Publishing, Cham, 2024.

[3] H. F. Baker, *Alternants and continuous groups*, Proceedings of the London Mathematical Society **s2-3** (1905), no. 1, 24-47.

[4] O. Brander and B. DeFacio, *A generalisation of Slepian's solution for the singular value decomposition of filtered Fourier transforms*, Inverse Problems **2** (May 1986), no. 2, L9-L14.

[5] Christophe Cheverry and Nicolas Raymond, *A Guide to Spectral Theory: Applications and Exercises*, Birkhäuser Advanced Texts Basler Lehrbücher, Springer International Publishing, Cham, 2021.

[6] Arnaud Debussche and Erwan Faou, *Weak Backward Error Analysis for SDEs*, SIAM Journal on Numerical Analysis **50** (January 2012), no. 3, 1735-1752.

[7] Erwan Faou, *Geometric numerical integration and Schrödinger equations*, Zurich Lectures in Advanced Mathematics, European Mathematical Society, Zürich, Switzerland, 2012.

[8] Charles Fefferman, *The Multiplier Problem for the Ball*, The Annals of Mathematics **94** (September 1971), no. 2, 330-336.

[9] F. Alberto Grünbaum, L Longhi, and M Perlstadt, *Differential operators commuting with finite convolution integral operators: Some non-Abelian examples*, SIAM Journal on Applied Mathematics **42** (1982), no. 5, 941-955.
<figure><img src="image_18.png" /><figcaption>FIGURE 12. 16 first eigenvectors obtained with a standard eigendecomposition of $\mathbf{K}(0)$, where $m_S(\varepsilon, \cdot) = \mathbf{1}_{\text{cat-head}(\varepsilon)}$ and $\widehat{m_F}(\varepsilon, \cdot) = \mathbf{1}_{B(0,0.3\cdot2\pi\mu(\varepsilon))}$. Here, $m_S$ is outlined in gray. $N = 60$, $\eta = 10^{-6}$.</figcaption></figure>

[10] Ernst Hairer, Christian Lubich, and Gerhard Wanner, *Geometric numerical integration: Structure-preserving algorithms for ordinary differential equations*, 2nd ed, Springer Series in Computational Mathematics, Springer, Berlin; New York, 2006.

[11] Tobias Jahnke and Christian Lubich, *Error Bounds for Exponential Operator Splittings*, Bit Numerical Mathematics **40** (2000), no. 4, 735–744.

[12] Jer-Nan Juang, Peiman Ghaemmaghami, and Kyong Been Lim, *Eigenvalue and eigenvector derivatives of a nondefective matrix*, Journal of Guidance, Control, and Dynamics **12** (July 1989), no. 4, 480–486.
<figure><img src="image_19.png" /><figcaption>FIGURE 13. Absolute value of the Fourier transform of the 16 first eigenvectors obtained with a standard eigendecomposition of $\mathbf{K}(0)$, where $m_S(\varepsilon, \cdot) = \mathbf{1}_{\text{cat-head}(\varepsilon)}$ and $\tilde{m}_F(\varepsilon, \cdot) = \mathbf{1}_{B(0,0.3\cdot2\pi\mu(\varepsilon))}$. Here, $\tilde{m}_F$ is outlined in gray. $N = 60, \eta = 10^{-6}$.</figcaption></figure>

[13] Santhosh Karnik, Justin Romberg, and Mark A. Davenport, *Improved bounds for the eigenvalues of prolate spheroidal wave functions and discrete prolate spheroidal sequences*, Applied and Computational Harmonic Analysis **55** (November 2021), 97–128.

[14] Santhosh Karnik, Zhihui Zhu, Michael B. Wakin, Justin Romberg, and Mark A. Davenport, *The fast Slepian transform*, Applied and Computational Harmonic Analysis **46** (May 2019), no. 3, 624–652.

[15] H. J. Landau and H. O. Pollak, *Prolate spheroidal wave functions, fourier analysis and uncertainty - II*, Bell System Technical Journal **40** (January 1961), no. 1, 65–84.
<figure><img src="image_20.png" /><figcaption>FIGURE 14. 16 first eigenvectors obtained with the varying masks procedure of $\mathbf{K}(0)$, where $m_S(\varepsilon, \cdot) = \mathbf{1}_{\text{cat-head}(\varepsilon)}$ and $\widehat{m_F}(\varepsilon, \cdot) = \mathbf{1}_{B(0,0.3\cdot2\pi\mu(\varepsilon))}$. Here, $m_S$ is outlined in gray. $N = 60$, $\eta = 10^{-6}$.</figcaption></figure>

[16] __________, *Prolate spheroidal wave functions, fourier analysis and uncertainty-III: The dimension of the space of essentially time- and band-limited signals*, Bell System Technical Journal **41** (July 1962), no. 4, 1295-1336.

[17] H.J Landau and H Widom, *Eigenvalue distribution of time and frequency limiting*, Journal of Mathematical Analysis and Applications **77** (October 1980), no. 2, 469-481.

[18] Kyong B. Lim, Jer-Nan Juang, and Peiman Ghaemmaghami, *Eigenvector derivatives of repeated eigenvalues using singular value decomposition*, Journal of Guidance, Control, and Dynamics **12** (March 1989), no. 2, 282-283.
<figure><img src="image_21.png" /><figcaption>FIGURE 15. Absolute value of the Fourier transform of the 16 first eigenvectors obtained with the varying masks procedure of $K(0)$, where $m_S(\varepsilon, \cdot) = \mathbf{1}_{\text{cat-head}(\varepsilon)}$ and $\widehat{m_F}(\varepsilon, \cdot) = \mathbf{1}_{B(0,0.3 \cdot 2\pi\mu(\varepsilon))}$. Here, $\widehat{m_F}$ is outlined in gray. $N = 60$, $\eta = 10^{-6}$.</figcaption></figure>

[19] Michael Reed and Barry Simon, *Methods of modern mathematical physics*, Rev. and enl. ed, Academic Press, New York, 1980.

[20] Frederik J. Simons and Dong V. Wang, *Spatiospectral concentration in the Cartesian plane*, GEM - International Journal on Geomathematics **2** (June 2011), no. 1, 1–36.

[21] D. Slepian, *Prolate spheroidal wave functions, fourier analysis, and uncertainty-V: The discrete case*, Bell System Technical Journal **57** (May 1978), no. 5, 1371–1430.

[22] D. Slepian and H. O. Pollak, *Prolate spheroidal wave functions, fourier analysis and uncertainty - I*, Bell System Technical Journal **40** (January 1961), no. 1, 43–63.
[23] David Slepian, *Prolate spheroidal wave functions, fourier analysis and uncertainty - IV: Extensions to many dimensions; generalized prolate spheroidal functions*, Bell System Technical Journal **43** (November 1964), no. 6, 3009–3057.

[24] L. Wang, *A review of prolate spheroidal wave functions from the perspective of spectral methods*, Journal of Mathematical Study **50** (June 2017), no. 2, 101–143.