# Keller and Lieb-Thirring estimates of the eigenvalues in the gap of Dirac operators

Jean Dolbeault, David Gontier, Fabio Pizzichillo, Hanne Van Den Bosch

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

## To cite this version:

Jean Dolbeault, David Gontier, Fabio Pizzichillo, Hanne Van Den Bosch. Keller and Lieb-Thirring estimates of the eigenvalues in the gap of Dirac operators. Revista Matemática Iberoamericana, inPress, 10.4171/RMI/1443. hal-03803758v2

HAL Id: hal-03803758

https://hal.science/hal-03803758v2

Submitted on 24 Jul 2023

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Keller and Lieb-Thirring estimates of the eigenvalues in the gap of Dirac operators

Jean Dolbeault, David Gontier, Fabio Pizzichillo and Hanne Van Den Bosch

**Abstract.** We estimate the lowest eigenvalue in the gap of the essential spectrum of a Dirac operator with mass in terms of a Lebesgue norm of the potential. Such a bound is the counterpart for Dirac operators of the *Keller estimates* for the Schrödinger operator, which are equivalent to Gagliardo-Nirenberg-Sobolev interpolation inequalities. Domain, self-adjointness, optimality and critical values of the norms are addressed, while the optimal potential is given by a Dirac equation with a Kerr nonlinearity. A new *critical* bound appears, which is the smallest value of the norm of the potential for which eigenvalues may reach the bottom of the gap in the essential spectrum. The Keller estimate is then extended to a Lieb-Thirring inequality for the eigenvalues in the gap. Most of our result are established in the *Birman-Schwinger* reformulation.

<table><tr><td>1</td><td>Introduction and main results</td><td>2</td></tr><tr><td>2</td><td>Properties of Dirac operators</td><td>8</td></tr><tr><td>2.1</td><td>A self-adjoint realization</td><td>8</td></tr><tr><td>2.2</td><td>The Birman-Schwinger operator</td><td>9</td></tr><tr><td>2.3</td><td>Proofs of Proposition 2.1, Lemma 2.3 and Proposition 2.4</td><td>11</td></tr><tr><td>3</td><td>The variational problem</td><td>15</td></tr><tr><td>3.1</td><td>An auxiliary maximization problem</td><td>15</td></tr><tr><td>3.2</td><td>Proof of Theorem 3.1</td><td>17</td></tr><tr><td>3.3</td><td>Regularity of the solutions of the non-linear Dirac equation</td><td>23</td></tr><tr><td>4</td><td>Lieb-Thirring inequality</td><td>24</td></tr><tr><td>5</td><td>Explicit computations</td><td>27</td></tr><tr><td>5.1</td><td>The case d = 1: proof of Theorem 1.3</td><td>27</td></tr><tr><td>5.2</td><td>The radial case in dimension d = 2</td><td>29</td></tr><tr><td>5.3</td><td>The radial case in dimension d = 3</td><td>30</td></tr><tr><td>5.4</td><td>An explicit bound in the radial case in dimensions d = 2 or d = 3</td><td>32</td></tr><tr><td>A</td><td>Open questions</td><td>33</td></tr><tr><td>B</td><td>Is the optimal potential radial? A numerical answer</td><td>34</td></tr><tr><td>C</td><td>A nonlinear interpolation inequality for the Dirac operator</td><td>34</td></tr><tr><td>C.1</td><td>Non-relativistic limit and Keller-Lieb-Thirring inequalities</td><td>34</td></tr><tr><td>C.2</td><td>An interpolation inequality for the Dirac operator</td><td>36</td></tr><tr><td>D</td><td>The case p = d = 1</td><td>38</td></tr><tr><td></td><td>References</td><td>41</td></tr></table>

*Mathematics Subject Classification (2020): Primary 81Q10; Secondary 49R05, 49J35, 47A75, 47B25.*

*Keywords:* Dirac operators, potential, spectral gap, eigenvalues, ground state, min-max principle, Birman-Schwinger operator, domain, self-adjoint operators, Keller estimate, Lieb-Thirring inequality, interpolation, Gagliardo-Nirenberg-Sobolev inequality, Kerr nonlinearity.
# 1. Introduction and main results

In 1961, J.B. Keller established in [45] the expression of the potential which minimizes the lowest eigenvalue, or ground state, $\lambda_S(V)$ of the Schrödinger operator $-\Delta - V$ in dimension $d = 1$, under a constraint on the Lebesgue norm

$$ \|V\|_p = \left( \int_{\mathbb{R}^d} |V|^p dx \right)^{1/p} $$

of exponent $p$ of $V$. This estimate was later extended in [53] by E.H. Lieb and W. Thirring to higher dimensions and to a sum of the lowest eigenvalues. During the last forty years, various refinements were published. As an example, we quote stability results for $\lambda_S(V)$ proved in [12] by E.A. Carlen, R.L. Frank, and E.H. Lieb. Although Dirac operators inherit many qualitative properties of Schrödinger operators, dealing with Dirac operators turns out to be a delicate issue.

If $\mathcal{D}_m$ denotes the free Dirac operator and $V$ is a non-negative valued function, $\mathcal{D}_m - V$ is not bounded from below. One is actually interested in the lowest eigenvalue $\lambda_D(V)$ in the essential gap $(-mc^2, mc^2)$, where $m$ denotes the mass and $c$ the speed of light. We shall speak of $\lambda_D(V)$ as the ground state energy of $\mathcal{D}_m - V$. In the standard setting, it is expected that $\lambda_D(V) - mc^2$ converges to $\lambda_S(V)$ in the non-relativistic limit, i.e., as $c \to +\infty$. It is therefore a natural question to estimate $\lambda_D(V)$ in terms of $\|V\|_p$ and identify the corresponding optimal potential. This question is the main purpose of our paper. A new critical value appears, which corresponds to the smallest value of $\|V\|_p$ for which $\lambda_D(V)$ reaches, for some potential $V \ge 0$, the lower end of the essential gap $-mc^2$. In a linear setting, a similar question has been raised in [34,35], where the authors find a critical value $\nu_1$ so that $\lambda_D(\mu * |\cdot|^{-1}) > -mc^2$ for all positive measures $\mu$ with $\mu(\mathbb{R}^3) < \nu_1$, with $2/(\pi/2 + 2/\pi) < \nu_1 \le 1$. Going back to [21, 27, 28], it is known that *Hardy inequalities* play an essential role in the analysis of the spectrum of Dirac-Coulomb operators. In the present article, except for the case $p = d = 1$, we rather find a nonlinear functional inequality of Gagliardo-Nirenberg-Sobolev nature, instead of a Hardy inequality (see comments in Appendix C.2).

It is possible to characterize the eigenvalues of $\mathcal{D}_m - V$ in the gap by a *min-max principle* according to [28–30] but this raises delicate issues involving the domain of the operator and its self-adjoint extensions addressed respectively in [30, 33, 36, 37, 63]. Applied with a Coulombian potential $V$, the method gives rise, after the maximising step in the min-max method, to a lower bounded quadratic form which amounts to a kind of Hardy inequality for the upper component: see [10, 21, 27] for details. The same strategy applies to a general potential $V$ under a constraint on $\|V\|_p$, except that the Keller type bound on $\lambda_D(V)$ is given by an implicit condition: see Appendix C. The optimal potential solves a nonlinear Dirac equation with Kerr-type nonlinearity. For the two-dimensional case, this equation has been studied in [5–8] by W. Borrelli. In the one-dimensional case, the solution is explicit, which allows us to identify it as in the case of the Schrödinger operator studied in [45]. Alternatively to the min-max principle, the properties of the *Birman-Schwinger operator* corresponding to $\mathcal{D}_m - V$ allows us to characterize $\lambda_D(V)$ and, except in Appendix C, we will adopt this point of view.
The Keller-Lieb-Thirring inequality for a Schrödinger operator goes as follows. Let us assume that $q > 2$, with $q < 2^* := 2d/(d-2)$ if $d \ge 3$, and let $\vartheta = d(q-2)/(2q)$. For any function $u \in H^1(\mathbb{R}^d)$, the Gagliardo-Nirenberg-Sobolev inequality

$$ \|\nabla u\|_2^{\vartheta} \|u\|_2^{1-\vartheta} \ge \mathcal{C}_q \|u\|_q $$

can be rewritten in the non-scale invariant form as

$$ (1.1) \quad \forall (\lambda, u) \in (0, +\infty) \times H^1(\mathbb{R}^d), \quad \|\nabla u\|_2^2 + \lambda \|u\|_2^2 \ge \mathcal{C}_q \lambda^{1-\vartheta} \|u\|_q^2 $$

with an optimal constant $C_q$ such that $\mathcal{C}_q^2 = \vartheta^{\vartheta}(1-\vartheta)^{1-\vartheta} C_q$. The equivalence of the two forms can be recovered by optimizing on $\lambda$ in (1.1). There is also an inequality which is dual of (1.1) and goes as follows. Consider a potential $V \in L^p(\mathbb{R}^d)$. Using Hölder's inequality with exponents $p$ and $q$ such that $1/p + 2/q = 1$ and $p > d/2$, and taking $\lambda$ so that $C_q \lambda^{1-\vartheta} = \|V\|_p$, we deduce from (1.1) that

$$ \int_{\mathbb{R}^d} |\nabla u|^2 dx - \int_{\mathbb{R}^d} V |u|^2 dx \ge \|\nabla u\|_2^2 - \|V\|_p \|u\|_q^2 \ge -(C_q^{-1} \|V\|_p)^{\frac{1}{1-\vartheta}} \|u\|_2^2 $$

This is the Keller-Lieb-Thirring estimate for $-\Delta - V$, i.e.,

$$ (1.2) \quad \forall V \in L^p(\mathbb{R}^d), \quad 0 \le \lambda_S^-(V) \le K_p \|V\|_p^\eta $$

where $\eta := 1/(1-\vartheta) = 2p/(2p-d)$ and $\lambda^- := \max(0, -\lambda)$ denotes the negative part of $\lambda$. See [22–24] for details. An optimization on $V$ shows that (1.1) and (1.2) are equivalent. The optimal constant in (1.2) is $K_p = C_q^{-\eta}$. In addition, for all $\lambda > 0$, if $u$ is a radial positive solution of

$$ (1.3) \quad -\Delta u - u^{\frac{p+1}{p-1}} = -\lambda u, $$

then $(u, \lambda)$ is an optimal pair for (1.1), and $V := u^{q-2} = u^{2/(p-1)}$ is an optimal potential for (1.2), which moreover satisfies $\lambda_S(V) = -\lambda$. It turns out that the solution of (1.3) is unique up to translations according to [16, 49, 55] and can be explicitly computed if $d=1$: see [45], or [23] and references therein for additional related results.

In order to state a Keller-Lieb-Thirring inequality for the Dirac operator, we need some definitions and preliminary properties. Let us start with the free Dirac operator on $\mathbb{R}^d$. We refer to [66] for a comprehensive list of results and properties. For simplicity, we choose units in which $c=1$, except in Appendix C in which we consider the non-relativistic limit as $c \to +\infty$. Let $d \ge 1$ and set $N := 2^{\lfloor(d+1)/2\rfloor}$ where $\lfloor x \rfloor = \max\{n \in \mathbb{Z} : n \le x\}$ denotes the integer part of $x$. Let $\alpha_1, \dots, \alpha_d$ and $\beta$ be $N \times N$ Hermitian matrices satisfying the following anti-commutation rules

$$ (1.4) \quad \forall j, k = 1, \dots, d, \quad \begin{cases} \alpha_j \alpha_k + \alpha_k \alpha_j = 2 \delta_{jk} \mathbb{I}_N \\ \alpha_j \beta + \beta \alpha_j = 0 \\ \beta^2 = \mathbb{I}_N \end{cases} $$

where $\delta_{jk}$ denotes the Kronecker symbol and $\mathbb{I}_N$ is the $N \times N$ identity matrix. See, e.g., [41] for an existence result for such matrices. The free Dirac operator in dimension $d$
is defined by

$$ \mathcal{D}_m := \sum_{j=1}^{d} \alpha_j (-i \partial_j) + m \beta = \alpha \cdot (-i \nabla) + m \beta $$

where we consider Cartesian coordinates $(x_1, \dots, x_d)$, $\partial_j := \partial/\partial x_j$ and $\alpha = (\alpha_k)_{k=1,\dots,d}$. With the Pauli matrices

$$ \sigma_1 := \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \quad \sigma_2 := \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} \quad \text{and} \quad \sigma_3 := \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}, $$

explicit expressions of $\mathcal{D}_m$ are given

(i) in dimension $d=1$, by $\alpha = \sigma_2$ and $\beta = \sigma_3$ so that

$$ \mathcal{D}_m := \sigma_2 (-i \partial_1) + m \sigma_3, $$

(ii) in dimension $d=2$, by $\alpha = (\sigma_j)_{j=1,2}$ and $\beta = \sigma_3$ so that

$$ \mathcal{D}_m := \sum_{j=1}^{2} \sigma_j (-i \partial_j) + m \sigma_3, $$

(iii) in dimension $d=3$, by $\alpha = (\alpha_k)_{k=1,2,3}$ and $\beta$ such that

$$ \alpha_k := \begin{pmatrix} 0 & \sigma_k \\ \sigma_k & 0 \end{pmatrix} \quad \text{and} \quad \beta := \begin{pmatrix} \mathbb{I}_2 & 0 \\ 0 & -\mathbb{I}_2 \end{pmatrix}. $$

The free Dirac operator satisfies $\mathcal{D}_m^2 = -\Delta + m^2$. It is self-adjoint on $L^2(\mathbb{R}^d, \mathbb{C}^N)$, with domain

$$ \mathrm{Dom}(\mathcal{D}_m) = \mathrm{H}^1(\mathbb{R}^d, \mathbb{C}^N) $$

and spectrum

$$ \sigma(\mathcal{D}_m) = \sigma_{\mathrm{ess}}(\mathcal{D}_m) = (-\infty, -m] \cup [m, +\infty). $$

Next we consider Dirac operators $\mathcal{D}_m - V$ with potentials $V \in L^p(\mathbb{R}^d, \mathbb{R}^+)$ where the notation $\mathcal{D}_m - V$ denotes $\mathcal{D}_m - V \mathbb{I}_N$. When switching on a potential $V$, we expect that some eigenvalues of $\mathcal{D}_m - V$ emerge from the upper essential spectrum $[m, +\infty)$. We shall prove in Section 2 that $\mathcal{D}_m - V$ can be defined as a self-adjoint operator with essential spectrum $\sigma_{\mathrm{ess}}(\mathcal{D}_m - V) = \sigma_{\mathrm{ess}}(\mathcal{D}_m)$. This allows us to define the ground state $\lambda_D(V)$ as the lowest eigenvalue in the gap $(-m, m)$.

Our first result states that the ground state is bounded by a function of $\|V\|_p$. Let

$$ (1.5) \qquad \Lambda_D(\alpha, p) := \inf \left\{ \lambda_D(V) : V \in L^p(\mathbb{R}^d, \mathbb{R}^+) \text{ and } \|V\|_p = \alpha \right\}. $$

**Theorem 1.1.** Assume that $p \ge d \ge 1$. There exists $\alpha_*(p) > 0$ such that the map $\alpha \mapsto \Lambda_D(\alpha, p)$ defined on $[0, \alpha_*(p))$ is continuous, strictly decreasing, takes values in $(-m, m]$, and such that

$$ \lim_{\alpha \to 0_+} \Lambda_D(\alpha, p) = m \quad \text{and} \quad \lim_{\alpha \to \alpha_*(p)} \Lambda_D(\alpha, p) = -m. $$
Moreover, if $(p, d) \neq (1, 1)$, the infimum (1.5) is attained on $(0, \alpha_*(p))$ and

$$ \forall \alpha \in (0, \alpha_*(p)), \quad \Lambda_D(\alpha, p) = \lambda_D(V_{\alpha, p}) $$

where $V_{\alpha,p} = |\Psi|^{2/(p-1)}$, and $\Psi \in L^2(\mathbb{R}^d, \mathbb{C}^N)$ solves the nonlinear Dirac equation

$$ (1.6) \qquad \mathcal{D}_m \Psi - |\Psi|^{\frac{2}{p-1}} \Psi = \Lambda_D(\alpha, p) \Psi $$

and satisfies the constraint $\int_{\mathbb{R}^d} |\Psi|^{2p/(p-1)} dx = \|V_{\alpha,p}\|_p^p = \alpha^p$.

The proof of Theorem 1.1 is given in Section 3 and relies on the properties of the inverse map of $\alpha \mapsto \Lambda_D(\alpha, p)$ defined by

$$ (1.7) \qquad \alpha_D(\lambda, p) := \inf \left\{ \|V\|_p : V \in L^p(\mathbb{R}^d, \mathbb{R}^+) \text{ and } \lambda_D(V) = \lambda \right\}. $$

The critical value is $\alpha_*(p) = \lim_{\lambda \to (-m)_+} \alpha_D(\lambda, p)$. It is such that

$$ \lim_{\alpha \to \alpha_*(p)_-} \lambda_D(V_{\alpha, p}) = -m $$

and this limit is the upper bound of the lower essential spectrum $(-\infty, -m]$ or, equivalently, the lower end of the gap. For sake of simplicity, we adopt the convention that $\alpha_*(p) = \alpha_D(-m, p)$. In the subcritical range of potentials, a simple consequence of Theorem 1.1 is the following *Keller-Lieb-Thirring estimate* for the Dirac operator $\mathcal{D}_m - V$.

**Corollary 1.2.** Assume $p \ge d \ge 1$. For all $V \in L^p(\mathbb{R}^d, \mathbb{R}^+)$ with $\|V\|_p < \alpha_*(p)$, we have the optimal bound

$$ (1.8) \qquad -m \le \Lambda_D (\|V\|_p, p) \le \lambda_D(V) \le m. $$

If $(p, d) \neq (1, 1)$, then $V_{\alpha,p}$ as in Theorem 1.1 realizes the equality case, i.e., $\lambda_D(V_{\alpha,p}) = \Lambda_D(\alpha, p)$.

Some plots of $\alpha \mapsto \Lambda_D(\alpha, p)$ are displayed in Fig. 1 (Right).

The nonlinear Dirac equation (1.6) plays for the Dirac operator $\mathcal{D}_m - V$ the same role as (1.3) for the Schrödinger operator $-\Delta + V$. However, $\Lambda_D(\alpha, V)$ is not obtained as the infimum but as a critical point of a Rayleigh quotient with infinitely many negative directions corresponding to a min-max principle (see [28]) and for this reason there is no simple interpolation inequality such as (1.1) in the case the Dirac operator. A more involved functional inequality holds: see Appendix C.

Nonlinear Dirac equations have been introduced to model extended fermions, as effective operators for nonlinear effects in graphene-like materials or Bose-Einstein condensates: see [32, Section 1.6] and [5, Introduction] for an introduction to the literature. Since the spinors in the Dirac equation have at least two components, many types of nonlinearities can be considered (see, e.g., [59] and references therein) and give rise to various phenomena. For instance, localized solutions to a nonlinear equation of the form

$$ \mathcal{D}_m \Psi - G(\Psi) = \lambda \Psi $$
for some function $G: \mathbb{C}^N \to \mathbb{C}^N$ correspond to solitary wave solutions to the time-dependent nonlinear Dirac equation and have attracted considerable attention: see, e.g., [4, 15, 38, 56].

It is a common assumption to consider a nonlinearity that preserves Lorentz, or particle-hole, symmetry. Such a non-linearity takes the form

$$ (1.9) \qquad \mathcal{D}_m \Psi - F(\langle \Psi, \beta \Psi \rangle_{\mathbb{C}^N}) \Psi = \lambda \Psi $$

and is called the *Soler-type nonlinearity*. The Soler nonlinearity formally appears when minimizing the first positive eigenvalue of $\mathcal{D}_m - \beta V$ but will not be studied in this paper. In contrast, the nonlinearity that appears in (1.6) is of the form $F(\langle \Psi, \Psi \rangle_{\mathbb{C}^N}) \Psi$, which is sometimes called a *Kerr-type nonlinearity* as in [5], apparently by extension of the cubic nonlinearity used in optics. Existence of localized solutions for (1.6) is studied in [5] in the critical exponent case $p = d = 2$, and in [8] in the critical exponent case $p = d$ for all dimensions $d \in \mathbb{N}$ with $m = 0$. Our results give an independent proof of the existence of a localized solution.

In [4, 38], the authors proved that equations of the form (1.9) have many solutions if $d \ge 2$ by looking for solutions of (1.9) in subspaces of fixed angular momentum. It seems that similar techniques could also be applied to (1.6). While it is reasonable to expect that the optimal potential is radially symmetric and the corresponding *ground state* $\Psi$ is the solution with lowest positive angular momentum and smallest number of oscillations, this is so far an open question: see Appendix A. In Appendix B, we also give numerical results that point in this direction.

We now focus on the dimension $d = 1$. It turns out that one can completely solve (1.6) using special functions. Explicit formulae are given below, where $B$ and ${}_2F_1$ respectively denote the Euler *Beta* function and the hypergeometric function.

**Theorem 1.3.** Let $d = 1$ and $p \in (1, +\infty)$. For all $\lambda \in [-m, m]$, the equation

$$ \mathcal{D}_m \Psi - |\Psi|^{\frac{2}{p-1}} \Psi = \lambda \Psi \quad \text{with} \quad \mathcal{D}_m := \begin{pmatrix} m & \partial_x \\ -\partial_x & -m \end{pmatrix} $$

has a unique solution $\Psi \in L^2(\mathbb{R}, \mathbb{C}^2) \setminus \{0\}$, up to a phase factor and a translation. Up to a translation, $V = |\Psi|^{2/(p-1)}$ is even, decreasing on $\mathbb{R}^+$ and such that $\alpha_D(\lambda, p) = \|V\|_p$.

* Subcritical regime $\lambda > -m$. With $A := \frac{p}{p-1}(m^2 - \lambda^2)$, $B := \frac{2}{p-1}\sqrt{m^2 - \lambda^2}$ and $z_0 := \frac{m-\lambda}{m+\lambda}$, we have

$$ (1.10) \qquad \forall x \in \mathbb{R}, \quad V(x) = \frac{A}{m \cosh(\mathbf{B} x) + \lambda} $$

and $(\alpha_D(\lambda, p))^p = p^p (\frac{m+\lambda}{p-1})^{p-1} z_0^{p-1/2} B(\frac{1}{2}, p) {}_2F_1(\frac{1}{2}, p; p + \frac{1}{2}; -z_0)$.

* Critical case $\lambda = -m$. With $\zeta = 2m/(p-1)$, we have

$$ (1.11) \qquad \forall x \in \mathbb{R}, \quad V(x) := \frac{\zeta p}{1 + \zeta^2 x^2} $$

and $(\alpha_*(p))^p = p^p (\frac{2m}{p-1})^{p-1} B(\frac{1}{2}, p - \frac{1}{2})$.
If $p = d = 1$, then $\alpha_D(\lambda, 1) = \arccos(\lambda/m)$ and

$$ \lim_{\lambda \to (-m)_+} \alpha_D(\lambda, 1) = \lim_{p \to 1_+} \alpha_*(p) = \pi. $$

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Let $d = 1$ and $m = 1$. (Left) The function $p \mapsto \alpha_*(p)$, with maximum at $p \approx 1.32$, satisfies $\lim_{p \to 1_+} \alpha_*(p) = \pi$ and $\lim_{p \to +\infty} \alpha_*(p) = 2$. (Right) For various values of $p$, the maps $\alpha \mapsto \Lambda_D(\alpha, p)$ take value $-1$ at $\alpha = \alpha_*(p)$. Upper (resp. lower) right plots correspond to $p < 1.32$ (resp. $p > 1.32$).</figcaption></figure>

See Fig. 1. With the notations of Theorem 1.1 and $\alpha = \alpha_D(\lambda, p)$, up to translations, we know that $V = V_{\alpha,p}$ in (1.10) and (1.11). For the proof of Theorem 1.3 and some additional details, see Section 5.1. Formally as $p \to 1_+$, the potential given by (1.10) converges to a $\delta$-Dirac distribution at $x = 0$ of mass $\arccos(\lambda/m)$ (see [67] for the study of self-adjoint extensions of $\mathcal{D}_m - \alpha \delta_0$). A remarkable consequence of the estimate in the case $p = d = 1$ is the Keller-Lieb-Thirring inequality

$$ (1.12) \qquad m \cos (\|V\|_1) \le \lambda_D(V) \le m $$

for any nonnegative potential $V \in L^1(\mathbb{R})$ with $\|V\|_1 \le \pi$. See Appendix D for a result on optimality cases in the case $p = d = 1$, with a proof.

The case of Theorem 1.3 presents some similarities with the results of [35]: in the case $p = 1$, it is expected that optimality is achieved only by singular measures. Our goals differ from those of [35] as we adopt the point of view of functional interpolation inequalities with Keller-type estimates as a subproduct, while [35] is concerned with the issue of the optimal charge distribution for a Dirac-Coulomb equation. In terms of methods, there are many similarities since we use Birman-Schwinger reformulations as well as classical tools of the concentration-compactness method. However, there are also significant differences because requesting that the potential is in $L^p(\mathbb{R}^d)$ means that the optimal $V$ is obtained through a non-linear Dirac equation which is not measure-valued as soon as $p > 1$.
Our results are not limited to estimates for the ground state and we also have a Lieb-Thirring inequality for the sum of eigenvalues in the gap $(-m, m)$ of Dirac operators of the form $\mathcal{D}_m - V$ with $V \in L^p(\mathbb{R}^d, \mathbb{R}^+)$. We denote by $-m < \lambda_1 \le \lambda_2 \le \cdots < m$ the possibly infinite sequence of eigenvalues in the gap $(-m, m)$, and write

$$e_k = e_k(m, V) := (m - \lambda_k) > 0,$$

so that $2m > e_1 \ge e_2 \ge \cdots > 0$. The quantity $e_k$ is the distance between the eigenvalue $\lambda_k$ and the bottom of the upper essential spectrum $+m$.

**Theorem 1.4.** For all $\gamma > d/2$ and $p \in (d, \gamma + d/2]$, there is a constant $L_{\gamma,d,p} > 0$ so that, for all $V \in L^p(\mathbb{R}^d, \mathbb{R}^+)$, and all $m > 0$, we have

$$(1.13) \quad \sum_{k \ge 1} e_k^\gamma(m, V) \le L_{\gamma,d,p} m^{\frac{d}{2}} \int_{\mathbb{R}^d} V_m^{\gamma+\frac{d}{2}-p} V^p dx \quad \text{with} \quad V_m := \min\{m, V\}.$$

If $V \in L^p(\mathbb{R}^d, \mathbb{R}^+) \cap L^{\gamma+d/2}(\mathbb{R}^d, \mathbb{R}^+)$, using the inequalities $V_m \le m$ and $V_m \le V$ gives respectively

$$\sum_{k \ge 1} e_k^\gamma \le L_{\gamma,d,p} m^{\gamma+d-p} \int_{\mathbb{R}^d} V^p dx \quad \text{and} \quad \sum_{k \ge 1} e_k^\gamma \le L_{\gamma,d,p} m^{\frac{d}{2}} \int_{\mathbb{R}^d} V^{\gamma+\frac{d}{2}} dx.$$

The inequality (1.13) is, in some sense, an interpolation between these two critical cases. In the proof, we use rough estimates: the method is constructive but there is a lot of space for improving on the constant $L_{\gamma,d,p}$.

**Structure of the paper.** This paper is organized as follows. In Section 2, we establish some properties of the operator $\mathcal{D}_m - V$ with $V \in L^p(\mathbb{R}^d)$: domain, associated Birman–Schwinger operator and self-adjointness. Section 3 is devoted to the variational problem associated with (1.5), after reformulation in the Birman–Schwinger framework. Theorem 3.1 is devoted to the existence of an optimal potential $V$ by concentration-compactness methods (Section 3.2). The regularity of the optimizers is studied in Section 3.3. Section 4 is devoted to the proof of Theorem 1.4. Explicit and numerical computations are performed in Section 5 in dimensions $d=1, 2$ and $3$. Open questions, numerical observations, remarks on the non-relativistic limit and Gagliardo-Nirenberg-Sobolev inequalities, and a result in the case $p=d=1$ are collected in Appendices A, B, C and D respectively.

## 2. Properties of Dirac operators

### 2.1. A self-adjoint realization

We assume that $V \in L^p(\mathbb{R}^d, \mathbb{R}^+)$ is positive valued and deal with the self-adjoint extensions of $\mathcal{D}_m - V$.

**Proposition 2.1.** Let $p \ge d \ge 1$ and $V \in L^p(\mathbb{R}^d, \mathbb{R}^+)$. Then the operator $\mathcal{D}_m - V$ is self-adjoint with domain:

$$\mathrm{Dom}(\mathcal{D}_m - V) := \{ \psi \in L^2(\mathbb{R}^d, \mathbb{C}^N) : \sqrt{V}\psi, (\mathcal{D}_m - V)\psi \in L^2(\mathbb{R}^d, \mathbb{C}^N) \}.$$
This is the unique self-adjoint realisation verifying:

$$H^1(\mathbb{R}^d, \mathbb{C}^N) \subseteq \mathrm{Dom}(\mathcal{D}_m - V) \subseteq H^{\frac{1}{2}}(\mathbb{R}^d, \mathbb{C}^N).$$

Moreover, we have the following properties.

(i) If $p$ satisfies

$$ (2.1) \qquad \left\{ \begin{array}{ll} p \ge 2 & \text{if } d=1, \\ p > 2 & \text{if } d=2, \\ p \ge d & \text{if } d \ge 3, \end{array} \right. $$

then $\mathrm{Dom}(\mathcal{D}_m - V) = H^1(\mathbb{R}^d, \mathbb{C}^N)$.

(ii) If $1 < p \le 2$ and $d=1$, then $\mathrm{Dom}(\mathcal{D}_m - V)$ is also included in $H^{\frac{3}{2}-\frac{1}{p}}(\mathbb{R}, \mathbb{C}^2)$, hence in $L^\infty(\mathbb{R}, \mathbb{C}^2)$.

We call the extension of Proposition 2.1 the distinguished extension, as it is the unique one whose domain is included in the formal form domain $\mathrm{Dom}(|\mathcal{D}_m|^{1/2}) = H^{1/2}(\mathbb{R}^d, \mathbb{C}^N)$. We will consider only this extension in what follows, so that the operator $\mathcal{D}_m - V$ is self-adjoint under the condition $p \ge d \ge 1$. The proof of the first part of Proposition 2.1 follows from [57]. For completeness, we provide a short proof using the associated Birman-Schwinger operator. Under Condition (2.1), the point (i) comes from the usual Kato-Rellich theorem [43,62]. The result in (ii) is derived by bootstrapping the Sobolev embedding theorem. See Section 2.3 for the proof of Proposition 2.1.

**Remark 2.2.** For comparison, it is interesting to consider limit cases. The Coulomb potential $V(x) = 1/|x|$ in dimension $d=3$ is in the weak Sobolev space $L^3_w(\mathbb{R}^d)$. The operator $\mathcal{D}_m - \kappa V$ is essentially self-adjoint if $0 \le \kappa \le \sqrt{3}/2$, it has a distinguished extension if $\sqrt{3}/2 < \kappa \le 1$, and no distinguished self-adjoint extension if $\kappa > 1$: see [13] and references therein. Also see Remark 2.5.

## 2.2. The Birman-Schwinger operator

The Birman-Schwinger operator is a powerful tool for analysing the spectral properties of $\mathcal{D}_m - V$ when $V$ belongs to a large class of perturbations. In the relativistic case, Klaus in [46] used it extensively to characterize and study the first eigenvalue of Dirac operators when proving the existence of a distinguished self-adjoint extension. For non-Hermitian potentials $V$, it can be employed to locate the eigenvalues of $\mathcal{D}_m - V$, as shown for example by Cuenin, Laptev and Tretter in [18], and by Fanelli and Krejčířík in [39]. Furthermore, it can be applied to discuss properties of the ground state of $\mathcal{D}_m - V$ when $V$ is a generalised Coulomb-type potential, see, e.g., [14,34,35,46]. Throughout this paper, following the approach by Kato [44] and by Konno and Kuroda [48], the Birman-Schwinger operator is used to define the self-adjoint extension of the operator $\mathcal{D}_m - V$. Then, with this rigorous definition at hand, we prove the existence of the optimization problem which defines the ground state by applying variational methods directly on the Birman-Schwinger reformulation of the problem.

For $z \notin \sigma(\mathcal{D}_m)$, let

$$ (2.2) \qquad R_0(z) := (\mathcal{D}_m - z)^{-1} $$
denote the resolvent operator. Recall that we assume $V \ge 0$. We introduce the *Birman-Schwinger operator*

$$ (2.3) \qquad K_V(z) := \sqrt{V} R_0(z) \sqrt{V} = \sqrt{V} \frac{1}{\mathcal{D}_m - z} \sqrt{V}. $$

**Lemma 2.3.** For all $p \ge d \ge 1$, all $V \in L^p(\mathbb{R}^d, \mathbb{R}^+)$ and all $z \notin \sigma(\mathcal{D}_m)$, the operator $K_V(z)$ is compact (hence bounded). In addition,

$$ \lim_{s \to \pm\infty} \|K_V(is)\|_{op} = 0. $$

This result follows from the *Kato-Seiler-Simon inequality*: see proof in Section 2.3. A consequence of Lemma 2.3 is the following result (also see Section 2.3 for the proof).

**Proposition 2.4.** Let $\mathcal{D}_m - V$ be the distinguished self-adjoint extension defined as in Proposition 2.1. Then

$$ \sigma_{\text{ess}}(\mathcal{D}_m - V) = \sigma_{\text{ess}}(\mathcal{D}_m) = (-\infty, -m] \cup [m, +\infty). $$

Moreover the Birman-Schwinger principle holds: for all $\lambda \in (-m, m)$, $\lambda$ is an eigenvalue of $\mathcal{D}_m - V$ if and only if 1 is an eigenvalue of $K_V(\lambda)$.

Let us point out some differences with Birman-Schwinger operators associated with Schrödinger operators (see Figure 2). In the Schrödinger case, the Birman-Schwinger operator is of the form

$$ \tilde{K}_V(\lambda) = \sqrt{V} \frac{1}{-\Delta - \lambda} \sqrt{V}. $$

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Birman-Schwinger principle. (Left). The spectrum of $\lambda \mapsto K_V(\lambda)$ (Dirac case) for $\lambda \in (-1, 1)$, and $V(x) = 2 \exp(-|x|^2/4)$ in dimension $d=2$. We only plotted the 10 largest (blue) and 10 lowest (red) eigenvalues. An energy $\lambda$ is an eigenvalue of $\mathcal{D}_m - V$ if one eigenvalue of $K_V(\lambda)$ crosses the black line 1. (Right) Same for $\lambda \mapsto \tilde{K}_V(\lambda)$ (Schrödinger case) with $\lambda \in (-2, 0)$.</figcaption></figure>

For any $\lambda < 0$, the operator $\tilde{K}_V(\lambda)$ is a positive compact operator and the map $\lambda \mapsto \tilde{K}_V(\lambda)$ is operator increasing on $\mathbb{R}^-$. In particular, if $\tilde{\mu}_1(\lambda) > \tilde{\mu}_2(\lambda) \ge \dots \ge 0$ denote
the eigenvalues of $\tilde{K}_V(\lambda)$, ranked in decreasing order and counted with multiplicities, all functions $\lambda \mapsto \tilde{\mu}_j(\lambda)$ are increasing on $\mathbb{R}^-$. In addition, the first eigenvalue $\tilde{\mu}_1$ is simple because the kernel $\tilde{K}_V(x, y)$ is pointwise positive, together with Krein-Rutman theorem: see [11, Theorem 6.13] for a statement and also [61, Section XIII.12].

In the Dirac case, the operator $K_V(\lambda)$ with $\lambda \in \mathbb{R}$ is defined only in the gap $(-m, m)$ of the essential spectrum. It is compact by Lemma 2.3 and symmetric because $\lambda$ is real, but it is not a positive operator. Its eigenvalues are real valued, and can be ranked as $\mu_1(\lambda) \ge \mu_2(\lambda) \ge \cdots \ge 0$ for the positive eigenvalues, and $\nu_1(\lambda) \le \nu_2(\lambda) \le \cdots \le 0$ for the negative ones. As the map $\lambda \mapsto (\mathcal{D}_m - \lambda)^{-1}$ is operator increasing on $(-m, m) \ni \lambda$, all maps $\lambda \mapsto \mu_j(\lambda)$ and $\lambda \mapsto \nu_j(\lambda)$ are increasing. This explains in particular why we expect eigenvalues to emerge from the upper essential spectrum in this setting. We do not know whether $\mu_1(\lambda)$ is always a simple eigenvalue or not (see Appendix A for more details on open questions).

For $\lambda \in (-m, m)$, $p \ge d \ge 1$, and $V \in L^p(\mathbb{R}^d, \mathbb{R}^+)$, let $\mu_1(K_V(\lambda))$ denote the largest (positive) eigenvalue of $K_V(\lambda)$. We rephrase the optimization problem (1.7) as

$$ (2.4) \qquad \alpha_D(\lambda, p) := \inf \left\{ \|V\|_p : V \in L^p(\mathbb{R}^d, \mathbb{R}^+) \text{ and } \mu_1(K_V(\lambda)) = 1 \right\}. $$

## 2.3. Proofs of Proposition 2.1, Lemma 2.3 and Proposition 2.4

We start by establishing that $K_V$ defined by (2.3) is a compact operator (Lemma 2.3) before proving Propositions 2.1 and 2.4.

*Proof of Lemma 2.3.* Assume first that $p > d$. We claim that, for $z \notin \sigma(\mathcal{D}_m)$, the operator $K_V(z)$ is compact. We have $R_0(z) = g_z(-i \nabla)$ with

$$ g_z(k) := \frac{1}{|k|^2 + m^2 - z^2} \left( \sum_{j=1}^{d} \alpha_j k_j + m \beta + z \mathbb{I}_N \right) =: \sum_{j=1}^{d} g_z^j(k) + g_z^m(k) + g_z^z(k), $$

with obvious notation. Let us focus on the $g_z^1(k)$ term. We write $g_z^1(k) = g_z^A(k) g_z^B(k)$ with

$$ g_z^A(k) := \frac{\sqrt{|k_1| \operatorname{sgn}(k_1)}}{\sqrt{|k|^2 + m^2 - z^2}} \alpha_1 \quad \text{and} \quad g_z^B(k) := \frac{\sqrt{|k_1|}}{\sqrt{|k|^2 + m^2 - z^2}} \mathbb{I}_d. $$

All components of the functions $g_z^A$ and $g_z^B$ are in $L^q(\mathbb{R}^d)$ for all $q > 2d$. Since $\sqrt{V}$ is in $L^{2p}(\mathbb{R}^d)$ with $2p > 2d$, we can use the Kato-Seiler-Simon inequality [65, Chapter 4, Theorem 4.1] and conclude that the operator

$$ K_V^1(z) := \sqrt{V} g_z^1(-i \nabla) \sqrt{V} $$

is in the Schatten class $\mathcal{S}_p(L^2(\mathbb{R}^d, \mathbb{C}^N))$ with

$$ \|K_V^1(z)\|_{\mathcal{S}_p} \le \left\| \sqrt{V(x)} g_z^A(-i \nabla) \right\|_{\mathcal{S}_{2p}} \left\| g_z^B(-i \nabla) \sqrt{V(x)} \right\|_{\mathcal{S}_{2p}} \le C \|V\|_p \|g_z^1\|_p. $$
In addition, we have $\|g_{z=is}^1\|_p \to 0$ as $s \to \pm\infty$. Similar computation for the other terms shows that $K_V$ is in the Schatten class $\mathfrak{S}_p$, and that $\lim_{s \to \pm\infty} \|K_V(i s)\|_{\text{op}} = 0$.

In the case $p = d$ with $d \ge 2$, we use that all components of the functions $g_z^A$ and $g_z^B$ are in the weak-Sobolev space $L_w^{2d}(\mathbb{R}^d)$. According to [65, Chapter 4], $g_z^A(-i \nabla) \sqrt{V}$ and $\sqrt{V} g_z^B(-i \nabla)$ are in the weak Schatten class $\mathfrak{S}_{2d,w}$. In particular, they are compact operators. This already proves that $K_V$ is compact as well. Note that $\|g_{z=i s}^1\|_{d,w}$ does not converge to 0 as $s \to \pm\infty$.

For any $R > 0$, $\sqrt{V_R} := \min(\sqrt{V}, R)$ belongs to $L^d(\mathbb{R}^d) \cap L^\infty(\mathbb{R}^d)$. We have

$$ \left\| g_z^A(-i \nabla) \sqrt{V} \right\|_{\text{op}} \le \left\| g_z^A(-i \nabla) \sqrt{V_R} \right\|_{\text{op}} + \left\| g_z^A(-i \nabla) (\sqrt{V} - \sqrt{V_R}) \right\|_{\text{op}}. $$

For $R$ large enough, the second term is small in the Schatten space $\mathfrak{S}_{2d,w}$, and for $z = i s$ with $|s|$ large, the first term is small in $\mathfrak{S}_q$ with $q > p$. Hence $\lim_{s \to \pm\infty} \|K_V(i s)\|_{\text{op}} = 0$.

Let us finally assume $p = d = 1$. In this case, with explicit computations, the kernel of the Birman-Schwinger operator $K_V(z)$ is given by

$$ \sqrt{V(x)} \cdot \frac{1}{2} \left( \begin{array}{cc} \frac{z+m}{k} & \operatorname{sign}(x-y) \\ \operatorname{sign}(x-y) & \frac{z-m}{k} \end{array} \right) e^{-k|x-y|} \cdot \sqrt{V(y)} \in L^2(\mathbb{R} \times \mathbb{R}, \mathbb{C}^2) $$

where $k = \sqrt{m^2 - z^2}$ is chosen with a positive real part. Thus, $K_V$ is a Hilbert-Schmidt operator (hence it is compact) and by the dominated convergence theorem we can conclude that $\lim_{s \to \pm\infty} \|K_V(i s)\|_{\text{op}} = 0$. ■

**Proof of Proposition 2.1.** We divide the proof in several steps.

**Distinguished self-adjoint extension.** We define the domain of self-adjointness for the operator $\mathcal{D}_m - V$ as a perturbation of $\mathcal{D}_m$ by applying the method of G. Nenciu in [57]. Using similar techniques as in the proof of Lemma 2.3, one can show that the operators $R_0(z) \sqrt{V}$ and $\sqrt{V} R_0(z)$ can be extended into bounded linear operators on $L^2(\mathbb{R}^d, \mathbb{C}^N)$. These operators are compact operators, in the Schatten class $\mathfrak{S}_{2p}$. We are now in the setting of [57]. Let $\Omega := \{z \in \mathbb{C} : 1 \notin \sigma(K_V(z))\}$ where $K_V$ is defined by (2.3). The set $\Omega$ is non-empty by Lemma 2.3. For $z \in \Omega$, define

$$ R(z) := R_0(z) + R_0(z) \sqrt{V} (1 - K_V(z))^{-1} \sqrt{V} R_0(z). $$

According to [57], the operator $\mathcal{D}_m - V$ has a unique self-adjoint extension whose resolvent is the operator $R(z)$ defined in (2.2). Its domain is $\text{Dom}(\mathcal{D}_m - V) := \text{Ran} R(z)$, which is independent of $z \in \Omega$. This is the unique extension which is included in the formal form domain $\text{Dom}(|\mathcal{D}_m|^{1/2}) = \mathbb{H}^{1/2}(\mathbb{R}^d, \mathbb{C}^N)$.

**Domain of the distinguished extension.** Define the maximal domain as

$$ \text{Dom}_{\max}(\mathcal{D}_m - V) := \{ \psi \in L^2(\mathbb{R}^d, \mathbb{C}^N) : (\mathcal{D}_m - V) \psi \in L^2(\mathbb{R}^d, \mathbb{C}^N) \}. $$

Then, the set

$$ \{ \psi \in L^2(\mathbb{R}^d, \mathbb{C}^N) : \sqrt{V} \psi, (\mathcal{D}_m - V) \psi \in L^2(\mathbb{R}^d, \mathbb{C}^N) \} $$
is also $\text{Dom}_{\max}(\mathcal{D}_m - V) \cap \text{Dom}(\sqrt{V})$. We write $\psi \in \text{Dom}(\mathcal{D}_m - V)$ as $R(z) f$ for some $f \in L^2(\mathbb{R}^d, \mathbb{C}^N)$. Then

$$ \sqrt{V} \psi = \sqrt{V} R_0(z) f + K_V (1 - K_V(z))^{-1} \sqrt{V} R_0(z) f \in L^2(\mathbb{R}^d, \mathbb{C}^N). $$

This proves that $\text{Dom}(\mathcal{D}_m - V) \subset \text{Dom}_{\max}(\mathcal{D}_m - V) \cap \text{Dom}(\sqrt{V})$. For the opposite inclusion, consider $\psi \in \text{Dom}_{\max}(\mathcal{D}_m - V) \cap \text{Dom}(\sqrt{V})$. We set $f := (\mathcal{D}_m + V - z) \psi \in L^2(\mathbb{R}^d, \mathbb{C}^N)$ and $\psi_0 = R(z) f \in \text{Dom}(\mathcal{D}_m - V)$. Note that $\sqrt{V} \psi_0 \in L^2(\mathbb{R}^d, \mathbb{C}^N)$, since

$$ \sqrt{V} R(z) = \sqrt{V} R_0(z) + K_V(z) (1 - K_V(z))^{-1} \sqrt{V} R_0(z) $$

is a bounded operator on $L^2(\mathbb{R}^d, \mathbb{C}^N)$. So $\phi := \psi - \psi_0$ is such that

$$ (\mathcal{D}_m - V - z) \phi = 0 \quad \text{and} \quad \phi \in \text{Dom}(\sqrt{V}). $$

From the relation

$$ (1 - R_0(z))^{-1} \sqrt{V} R_0(z) (\mathcal{D}_m - V - z) = \sqrt{V}, $$

we obtain $\sqrt{V} \phi = 0$, and from the relation

$$ R_0(z) (\mathcal{D}_m - V - z) = 1 - R_0(z) V = 1 - (R_0(z) \sqrt{V}) \sqrt{V}, $$

we finally get $\phi = 0$, hence $\psi = \psi_0 \in \text{Dom}(\mathcal{D}_m - V)$.

Finally, if $p \ge d \ge 1$ by the Hölder inequality and the Sobolev embedding theorem we get that $H^1(\mathbb{R}^d, \mathbb{C}^N) \subset \text{Dom}(\mathcal{D}_m - V)$ and this concludes the first part of the proof.

**Self-adjointness on** $H^1(\mathbb{R}^d, \mathbb{C}^N)$. Let us prove (i). Assume that $p$ satisfies (2.1). Thanks to (1.4) we have

$$ \forall \psi \in H^2(\mathbb{R}^d, \mathbb{C}^N), \quad \|\mathcal{D}_m \psi\|_2^2 = \|\nabla \psi\|_2^2 + m^2 \|\psi\|_2^2. $$

This shows that the graph norm of $\mathcal{D}_m$ is equivalent to the usual $H^1(\mathbb{R}^d, \mathbb{C}^N)$ norm. Set

$$ q := \begin{cases} \frac{2p}{p-2} & \text{if } p > 2 \\ +\infty & \text{if } d = 1 \text{ and } p = 2 \end{cases} \quad \text{so that} \quad \frac{1}{p} + \frac{1}{q} = \frac{1}{2}. $$

We write $V = V_1 + V_2$ with $V_1 := V 1_{V \ge R}$ and $V_2 := V 1_{V \le R}$. We have

$$ \begin{aligned} \|V \psi\|_2 &\le \|V_1 \psi\|_2 + \|V_2 \psi\|_2 \le \|V_1\|_p \|\psi\|_q + \|V_2\|_{\infty} \|\psi\|_2 \\ &\le C_S \|V_1\|_p \|\psi\|_{H^1} + \|V_2\|_{\infty} \|\psi\|_2 \end{aligned} $$

where, in the last inequality, we used Sobolev's embedding $H^1(\mathbb{R}^d) \hookrightarrow L^q(\mathbb{R}^d)$ and, according to (2.1), $q$ satisfies

$$ \begin{cases} 2 \le q \le +\infty & \text{if } d=1, \\ 2 \le q < +\infty & \text{if } d=2, \\ 2 \le q \le \frac{2d}{d-2} & \text{if } d \ge 3. \end{cases} $$
We choose $R$ large enough so that $C_S \|V_1\|_p < 1$ and conclude with the *Kato-Rellich theorem* (see [60, Theorem X.12]) that $\mathcal{D}_m - V$ is self-adjoint with domain $H^1(\mathbb{R}^d, \mathbb{C}^N)$. Since any self-adjoint operator only admits trivial self-adjoint extensions, we can conclude that $H^1(\mathbb{R}^d, \mathbb{C}^N) = \text{Dom}(\mathcal{D}_m - V)$.

**Regularity for** $d=1$. Let us focus on (ii) and assume that $d=1$ and $1<p\le 2$. Let us prove that $\text{Dom}(\mathcal{D}_m - V)$ is also included in $H^{\frac{3}{2}-\frac{1}{p}}(\mathbb{R}, \mathbb{C}^2)$. For any $\psi \in \text{Dom}(\mathcal{D}_m - V)$, we have

$$ (\mathcal{D}_m - V) \psi := f \in L^2(\mathbb{R}, \mathbb{C}^2), $$

hence $\mathcal{D}_m \psi = f + V\psi$. We recall the following negative Sobolev embeddings: for all $1<r\le 2$, we have $L^r(\mathbb{R}) \hookrightarrow H^{-s}(\mathbb{R})$ for all $s \ge \frac{2-r}{2r}$ and $L^2(\mathbb{R}) \hookrightarrow H^{-s}(\mathbb{R})$ for all $s \ge 0$, while

$$ \forall s \ge \frac{1}{2p}, \quad V \psi = \sqrt{V} \underbrace{(\sqrt{V} \psi)}_{\in L^2(\mathbb{R}, \mathbb{C}^2)} \in L^{\frac{2p}{p+1}}(\mathbb{R}, \mathbb{C}^2) \hookrightarrow H^{-s}(\mathbb{R}, \mathbb{C}^2). $$

We deduce that $\mathcal{D}_m \psi \in H^{-\frac{1}{2p}}(\mathbb{R}, \mathbb{C}^2)$, hence that $\psi \in H^{1-\frac{1}{2p}}(\mathbb{R}, \mathbb{C}^2)$. We now bootstrap the argument. For $p>1$, we have $1-\frac{1}{2p} > \frac{1}{2}$, so $\psi \in L^\infty(\mathbb{R}, \mathbb{C}^2)$ by Sobolev's embedding. This gives $V \psi \in L^p(\mathbb{R}, \mathbb{C}^2) \hookrightarrow H^{\frac{p-2}{2p}}(\mathbb{R}, \mathbb{C}^2)$. So $\mathcal{D}_m \psi = f - V \psi \in H^{\frac{p-2}{2p}}(\mathbb{R}, \mathbb{C}^2)$ as well, and we obtain $\psi \in H^{1+\frac{p-2}{2p}}(\mathbb{R}, \mathbb{C}^2)$ with $1+\frac{p-2}{2p} = \frac{3}{2}-\frac{1}{p}$, as wanted. ■

*Proof of Proposition 2.4.* Since $R_0(z) \sqrt{V}$ is compact, then $R(z)$ is a compact perturbation of the free resolvent $R_0(z)$. The result on $\sigma_{ess}(\mathcal{D}_m - V)$ follows from [66, Theorem 4.5] (also see [61, Theorem XIII.14 and Corollary 1]). Such a result is known in the literature as *Weyl's theorem*.

Moreover, by construction, the Birman-Schwinger principle holds for the distinguished self-adjoint extension defined as in Proposition 2.1: $\lambda \in (-m, m)$ is an eigenvalue of $\mathcal{D}_m - V$ if and only if 1 is an eigenvalue of $K_V(\lambda)$. See [2, Theorem 1.3] for a similar application of the Birman-Schwinger principle in a non-relativistic setting. ■

**Remark 2.5.** The self-adjointness of Dirac operators involving potentials with one Coulomb singularity or several Coulomb singularities has been intensively studied in respectively [3, 44, 57, 64, 69–71] (with additional references therein) and [46, 58]. In the alternative strategy of [36, 37] based on [28], a distinguished self-adjoint extension is built using the underlying Hardy inequality, which was related with the other constructions for Dirac-Coulomb operators in [33, 34]. Also see [31, 63] for further considerations on min-max principles, Hardy inequalities and self-adjointness issues. Optimal Hardy inequalities have been repeatedly used to establish optimal conditions for the existence of a ground state. For instance, in presence of a magnetic field as in [20, 25, 26], a critical magnetic field is obtained as the ground state energy approaches $-mc^2$, which determines the optimal constant of the corresponding Hardy inequality. In the approach of [34, 35] as well as in our paper, the Birman-Schwinger formula is essential as it was in [46, 47, 57]. Notice that we do not rely on Nenciu's method [57, Corollary 2.1], but instead use the method of Konno and Kuroda [48] and Kato's approach [44].
### 3. The variational problem

In this section, we consider the minimization problem (2.4) and prove Theorem 1.1 in a reformulation which relies on the Birman-Schwinger operator associated to $\mathcal{D}_m - V$, as introduced in Section 2.2. The proof of Theorem 1.1 is given below, right after the statement of Corollary 3.2, as a simple consequence of previous results in the Birman-Schwinger framework.

#### 3.1. An auxiliary maximization problem

First, we notice that, for all $t > 0$, we have $K_{tV}(\lambda) = t K_V(\lambda)$, hence $\mu_1(K_{tV}(\lambda)) = t \mu_1(K_V(\lambda))$. So, introducing the auxiliary problem

$$ (3.1) \qquad \mathcal{N}(\lambda, p) := \sup \left\{ \mu_1(K_W(\lambda)) : W \in L^p(\mathbb{R}^d, \mathbb{R}^+) , \|W\|_p = 1 \right\}, $$

we deduce that

$$ (3.2) \qquad \alpha_D(\lambda, p) = \frac{1}{\mathcal{N}(\lambda, p)}. $$

If $W$ is a maximizer for $\mathcal{N}(\lambda, p)$, then $V = W/\mathcal{N}(\lambda, p)$ is a minimizer for $\alpha_D(\lambda, p)$. In what follows, we study the maximization problem (3.1). We perform several changes of variables to study this problem. First, the min-max principle shows that $\mathcal{N}(\lambda, p)$ equals

$$ \mathcal{N}(\lambda, p) = \sup_{\substack{W \in L^p(\mathbb{R}^d, \mathbb{R}^+) \\ \|W\|_p = 1}} \sup_{\substack{\phi \in L^2(\mathbb{R}^d, \mathbb{C}^N) \\ \|\phi\|_2 = 1}} \left\langle \phi, \sqrt{W} R_0(\lambda) \sqrt{W} \phi \right\rangle. $$

We make the change of variable

$$ w := \sqrt{W} \phi $$

so that, by Hölder's inequality, $w \in L^{\frac{2p}{p+1}}(\mathbb{R}^d, \mathbb{C}^N)$, and, with the convention that $\|w\|_r = \|w|_{\mathbb{C}^N}\|_r$,

$$ \|w\|_{\frac{2p}{p+1}} \le \|W\|_{\frac{1}{p}} \|\phi\|_2 = 1. $$

In addition, there is equality if and only if $W^p$ is proportional to $|\phi|^2$, both proportional to $|w|^{\frac{2p}{p+1}}$. With

$$ q := \frac{2p}{p+1} \in (1, 2), $$

this shows that $\mathcal{N}(\lambda, p)$ is also solution to the optimization problem

$$ (3.3) \qquad \mathcal{N}(\lambda, p) = \sup \left\{ \langle w, R_0(\lambda) w \rangle : w \in L^q(\mathbb{R}^d, \mathbb{C}^N), \|w\|_q = 1 \right\}. $$

In addition, if $w \in L^q(\mathbb{R}^d, \mathbb{R}^+)$ is an optimizer of (3.3), then the corresponding optimal $W$ and $\phi$ are given by

$$ W = |w|^{\frac{q}{p}} = |w|^{\frac{2}{p+1}} \quad \text{and} \quad \phi = |w|^{\frac{q}{2}-1} w = |w|^{-\frac{1}{p+1}} w. $$
Thus, by showing the existence of an optimizer for (3.3), we solve problem (3.1), and by definition of the Birman-Schwinger operator, find an optimal potential and eigenfunction for our original problem (1.7).

Since $\alpha \mapsto \Lambda_D(\alpha, p)$ is the inverse map of $\lambda \mapsto \alpha_D(\lambda, p)$ according to (3.2), and since $\alpha_D(\lambda, p) = 1/\mathcal{N}(\lambda, p)$, it is enough to focus on the properties of $\mathcal{N}(\cdot, p)$.

**Theorem 3.1.** Let us consider $\mathcal{N}$ defined by (3.3). For all $\lambda \in (-m, m)$ and all $p > d$, we have $\mathcal{N}(\lambda, p) > 0$. All maximizing sequences for (3.3) are precompact up to translations, hence (3.3) has maximizers. If $w$ is such an optimizer, then $w$ satisfies the Euler-Lagrange equation

$$ (3.4) \qquad R_0(\lambda) w = \tau |w|^{-\frac{2}{p+1}} w \quad \text{with} \quad \tau = \mathcal{N}(\lambda, p). $$

Finally, the map $\lambda \mapsto \mathcal{N}(\lambda, p)$ is continuous, strictly increasing, and satisfies

$$ \lim_{\lambda \to -m} \mathcal{N}(\lambda, p) := \mathcal{N}_c(p) > 0 \quad \text{and} \quad \lim_{\lambda \to +m} \mathcal{N}(\lambda, p) = \infty. $$

The proof of the first part relies on the profile decomposition method (concentration-compactness) used by Lions [54], and is given in the next section. Theorem 3.1 implies the existence of an optimal potential and an optimal spinor.

**Corollary 3.2.** Under the assumptions of Theorem 1.1, the infimum (2.4) is attained for any $\lambda \in (-m, m)$ by a potential $V = |\Psi|^{2/(p-1)}$, where $\Psi \in L^2(\mathbb{R}^d, \mathbb{R}^N)$ solves the nonlinear Dirac equation

$$ (3.5) \qquad \mathcal{D}_m \Psi - |\Psi|^{\frac{2}{p-1}} \Psi = \lambda \Psi, $$

such that $\lambda_D(V) = \lambda$ and $(\int_{\mathbb{R}^d} |\Psi|^{2p/(p-1)} dx)^{1/p} = \|V\|_p = \alpha_D(\lambda, p) = 1/\mathcal{N}(\lambda, p)$.

*Proof of Theorem 1.1.* Our main Theorem 1.1 is a direct consequence of Theorem 3.1 and Corollary 3.2. Since $\mathcal{N}_c(p) > 0$, we have indeed $\alpha_c(p) := 1/\mathcal{N}_c(p) < \infty$. ■

*Proof of Corollary 3.2.* First, we translate the Euler-Lagrange equation for $w$ into an equation for the potential $V$ and an eigenfunction (not normalized) $\Psi$. We set

$$ \Psi = \tau^{\frac{1-p}{2}} |w|^{-\frac{2}{p+1}} w \quad \text{so that} \quad w = \tau^{\frac{p+1}{2}} |\Psi|^{\frac{2}{p-1}} \Psi. $$

Applying $\mathcal{D}_m - \lambda$ to (3.4) shows that $\Psi$ satisfies the nonlinear Dirac equation (3.5). The optimal potential $W$ for the $\mathcal{N}(\lambda, p)$ problem in (3.1) is $W = |w|^{\frac{2}{p+1}} = \tau |\Psi|^{\frac{2}{p-1}}$, and finally, the optimal potential $V$ for the $\alpha_D(\lambda, p)$ problem is, as wanted,

$$ V = \frac{W}{\mathcal{N}(\lambda, p)} = |\Psi|^{\frac{2}{p-1}}. $$

We recover the value of $\mathcal{N}(\lambda, p)$ and $\alpha_D(\lambda, p)$ from the solution $\Psi$ because

$$ \int_{\mathbb{R}^d} |\Psi|^{\frac{2p}{p-1}} dx = \tau^{-p} \int_{\mathbb{R}^d} |w|^{\frac{2p}{p+1}} dx = \tau^{-p} = \mathcal{N}(\lambda, p)^{-p} = \alpha_D(\lambda, p)^p. $$

Among all solutions of (3.5), $\Psi$ is the one with the smallest $L^{\frac{2p}{p-1}}(\mathbb{R}^d, \mathbb{C}^N)$ norm so that $\lambda = \Lambda_D(\alpha, p)$ and $\Psi$ actually solves (1.6). ■
## 3.2. Proof of Theorem 3.1

We now prove Theorem 3.1. We consider a more general case, and study a general optimization problem. In what follows, we use the notation

$$ \langle w, K * w \rangle := \iint_{\mathbb{R}^d \times \mathbb{R}^d} \langle w(x), K(x - y) w(y) \rangle_{\mathbb{C}^N} dx dy $$

and define for any $s > 0$ the maximization problem

$$ (3.6) \quad J(s) := \sup \left\{ \langle w, K * w \rangle : w \in L^q(\mathbb{R}^d, \mathbb{C}^N), \int_{\mathbb{R}^d} |w|^q dx = s \right\}. $$

Here, $K$ is a convolution operator, or equivalently a multiplication operator in Fourier space. In our case, $K(x - y) = R_0(\lambda)(x - y)$ is the kernel of the Dirac resolvent, but we state a more general result.

**Lemma 3.3.** Let $q \in (1, 2)$, set $q' := q/(q-1) \in (2, +\infty)$ and $r := q'/2 \in (1, +\infty)$. Let $K: \mathbb{R}^d \to \mathcal{M}_N(\mathbb{C})$ be a matrix-valued function satisfying $K(x) = K(-x)^*$, and such that one of the two properties holds:

(i) either $K \in L^r(\mathbb{R}^d, \mathcal{M}_N(\mathbb{C}))$,

(ii) or $K = R_0(\lambda)$ is a Dirac resolvent for some $\lambda \in (-m, m)$.

Then the map $w \mapsto \langle w, K * w \rangle$ is well-defined on $L^q(\mathbb{R}^d, \mathbb{C}^N)$ and real valued. Moreover, if $J(1) > 0$, then (3.6) admits maximizers.

Before proving this result, we make several remarks.

**Remark 3.4.** Lemma 3.3 fails at the endpoint $q = 2$. Indeed, by applying the Fourier transform we have

$$ \langle w, K * w \rangle = \int_{\mathbb{R}^d} \langle \hat{w}(k), \hat{K}(k) \hat{w}(k) \rangle_{\mathbb{C}^N} dk. $$

This means that all optimizing sequences must concentrate on Dirac masses in Fourier space at locations where $k \mapsto \sup \text{spec} (\hat{K}(k))$ has maxima. Since the Fourier transform is an isometry on $L^2(\mathbb{R}^d)$, we deduce that the maximization problem has no maximum in general. The same argument shows that the existence of optimizers is closely related to the fact that the Fourier transform is not a bijection between $L^q(\mathbb{R}^d)$ and $L^{q'}(\mathbb{R}^d)$ if $1 < q < 2$.

**Remark 3.5.** In the case of the Dirac operator, one has an explicit expression for $K = R_0(\lambda)$, the fundamental solution of $\mathcal{D}_m - \lambda$. Using that

$$ (\mathcal{D}_m - \lambda)^{-1} = (\mathcal{D}_m + \lambda) \frac{1}{-\Delta + m^2 - \lambda^2}, $$

we first deduce that $R_0(\lambda)(\cdot)$ is the Fourier transform of

$$ g_{\lambda}(k) = \left( \sum_{j=1}^{d} \alpha_j k_j + m \beta + \lambda \mathbb{I}_N \right) \frac{1}{k^2 + m^2 - \lambda^2}. $$
The function $k \mapsto g_\lambda(k)$ is analytic on $\mathbb{R}^d$ because there is no singularity in the denominator since $|\lambda| < m$, so its Fourier transform is exponentially decaying in $x$. Actually, we have

$$R_0(\lambda)(x) = \frac{c_{d,\lambda}}{|x|^{d/2-1}} \left( i \sum_{j=1}^{d} \alpha_j \frac{x_j}{|x|} \sqrt{m^2 - \lambda^2} K_{\frac{d}{2}} \left( \sqrt{m^2 - \lambda^2} |x| \right) + (m\beta + \lambda \mathbb{I}_N) K_{\frac{d}{2}-1} \left( \sqrt{m^2 - \lambda^2} |x| \right) \right)$$

where $c_{d,\lambda} = \frac{1}{2\pi} (\frac{\sqrt{m^2 - \lambda^2}}{2\pi})^{d/2-1}$ and $K_\nu$ is the modified Bessel function of the second kind. In particular, there is $C \ge 0$ so that

$$|R_0(\lambda)(x)| \le \begin{cases} C |x|^{1-d} & \text{as } |x| \to 0, \\ C e^{-\sqrt{m^2-\lambda^2}} |x| & \text{as } |x| \to +\infty. \end{cases}$$

So, in the Dirac case, we have $R_0(\lambda) \in L^r(\mathbb{R}^d)$ for all $r < \frac{d}{d-1}$ and $R_0(\lambda) \in L_w^{\frac{d}{d-1}}(\mathbb{R}^d)$. In particular, the case (ii) is not covered by (i) only in the case where $r = \frac{d}{d-1}$, which corresponds to the critical exponent case $p = d \ge 2$, that is, $q = \frac{2d}{d+2}$ in (3.3).

**Remark 3.6.** Let us consider the case $s = 1$ in (3.6). In order to see that $J(1) > 0$ in the Dirac case with $\lambda \in (-m, m)$, let $f \in L^q(\mathbb{R}^d, \mathbb{C})$ be a normalized function and let $\phi_+ \in \mathbb{C}^N$ be a normalized vector such that $\beta \phi_+ = \phi_+$. We find that

$$(\mathcal{D}_m + \lambda) f \phi_+ = (m + \lambda) f \phi_+ + (-i \nabla f) \cdot \alpha \phi_+.$$

Moreover, by (1.4), we have that $\langle \phi_+, \alpha_j \phi_+ \rangle_{\mathbb{C}^N} = 0$. Thus:

$$(3.7) \qquad J(1) \ge \langle f \phi_+, R_0(\lambda) f \phi_+ \rangle = \langle f \phi_+, (-\Delta + m^2 - \lambda^2)^{-1} (\mathcal{D}_m + \lambda) f \phi_+ \rangle \\ = (m + \lambda) \langle f, (-\Delta + m^2 - \lambda^2)^{-1} f \rangle_{L^2(\mathbb{R}^d, \mathbb{C})} > 0.$$

*Proof of Lemma 3.3.* First, we note that the condition $K(x) = K(-x)^*$ reads $\hat{K}(k) = \hat{K}(k)^*$, so the operator $K$ is symmetric.

In the first part of the proof, we cover both cases (i) and (ii) by assuming

$$(3.8) \quad K \in L_w^r(\mathbb{R}^d, \mathcal{M}_N(\mathbb{C})) \cap L^r(\mathcal{B}_1^c, \mathcal{M}_N(\mathbb{C})) \quad \text{with} \quad \mathcal{B}_R := \{x \in \mathbb{R}^d : |x| < R\}$$

with $\frac{2}{q} + \frac{1}{r} = 2$. From the Hardy-Littlewood-Sobolev inequality, and since $K \in L_w^r(\mathbb{R}^d)$, we have

$$(3.9) \quad \forall w_1, w_2 \in L^q(\mathbb{R}^d), \quad |\langle w_1, K * w_2 \rangle| \le C \|w_1\|_q \|w_2\|_q \|K\|_{r,w}.$$

In particular, $w \mapsto \langle w, K * w \rangle$ is well-defined and real valued on $L^q(\mathbb{R}^d)$.

Using the scaling $w_s = s^{1/q} w_1$, we obtain that

$$(3.10) \qquad J(s) = s^{2/q} J(1).$$
Since $J(1) > 0$, we deduce first that $J(s)$ is increasing. Also, since $2/q > 1$, $J(s)$ is convex and so we have the strong binding inequality

$$ (3.11) \qquad \forall s, s' > 0, \quad J(s + s') > J(s) + J(s'). $$

Let $(w_n)_{n \in \mathbb{N}}$ be a maximizing sequence for $J(1)$. Our argument relies on the concentration-compactness method for the sequence $(w_n)_{n \in \mathbb{N}}$, following the approach of Lions [54] and using Levy's functional. It differs from the concentration-compactness method used in [34], as we work directly with the Birman-Schwinger operator instead of the min-max quadratic form. We set

$$ Q(\rho) := \liminf_{n \to +\infty} Q_n(\rho) \quad \text{with} \quad Q_n(\rho) := \sup_{x \in \mathbb{R}^d} \int_{\mathcal{B}(x, \rho)} |w_n|^q dx. $$

It is clear from the definition that $\rho \mapsto Q(\rho)$ is non-decreasing, and that $Q(\rho) \le 1$ for all $\rho > 0$. We set

$$ \mu := \lim_{\rho \to +\infty} Q(\rho) \in [0, 1]. $$

We divide the proof in the classical steps of the concentration-compactness method and start by discarding the cases $\mu = 0$ (vanishing) and $\mu < 1$ (dichotomy).

• **Vanishing.** Fix $\varepsilon := J(1)/4 > 0$. Since $K \in L^r(\mathcal{B}_1^c)$, there is $R > 1$ large enough so that

$$ \|K\|_{L^r(\mathcal{B}_1^c)} \le \varepsilon. $$

By Young's inequality, since $\frac{2}{q} + \frac{1}{r} = 2$, we get that for all $w \in L^q(\mathbb{R}^d)$ with $\|w\|_q = 1$,

$$ \langle w, (1_{\mathcal{B}_R^c} K) * w \rangle \le \varepsilon \|w\|_q^2 \le \varepsilon. $$

We now estimate the contribution of $1_{\mathcal{B}_R} K$. For $z \in \mathbb{Z}^d$, let $C_z$ be the cube $z + [0, 1]^d$, so that $\{C_z\}_{z \in \mathbb{Z}^d}$ covers $\mathbb{R}^d$. For a function $w : \mathbb{R}^d \to \mathbb{C}^N$, we have

$$ \begin{aligned} \langle w, (1_{\mathcal{B}_R} K) * w \rangle &= \sum_{z, z' \in \mathbb{Z}^d} \iint_{C_z \times C_{z'}} \langle w(x), (1_{\mathcal{B}_R} K)(x-y) w(y) \rangle_{\mathbb{C}^N} dx dy \\ &\le \|K\|_{L_w^r(\mathbb{R}^d)} \sum_{z, z' \in \mathbb{Z}^d} \|w\|_{L^q(C_z)} \|w\|_{L^q(C_{z'})} \mathbb{1}_{\{|z-z'| \le R+2 \sqrt{d}\}} \end{aligned} $$

using again the Hardy-Littlewood-Sobolev inequality. The double sum can be seen as a discrete convolution, and we apply Young's inequality with $z \mapsto \|w\|_{L^q(C_z)} \in \ell^2(\mathbb{Z}^d)$ and $z \mapsto \mathbb{1}_{\{|z| \le R+2 \sqrt{d}\}} \in \ell^1(\mathbb{Z}^d)$ to bound

$$ \langle w, (1_{\mathcal{B}_R} K) * w \rangle \le C_R \sum_{z \in \mathbb{Z}^d} \|w\|_{L^q(C_z)}^2 \le C_R \sup_{z \in \mathbb{Z}^d} \|w\|_{L^q(C_z)}^{2-q} \|w\|_q^q $$

where $C_R$ is a positive constant which is independent of $w$: for all $w \in L^q(\mathbb{R}^d)$ with $\|w\|_q = 1$, we have

$$ \langle w, K * w \rangle \le \varepsilon + C_R \sup_{z \in \mathbb{Z}^d} \|w\|_{L^q(C_z)}^{2-q}. $$
Applying this estimate to a maximizing sequence $(w_n)_{n \in \mathbb{N}}$ for $J(1) = 4\varepsilon$, we obtain that, up to a subsequence,

$$ \frac{1}{2} J(1) \leq \langle w_n, K * w_n \rangle \leq \frac{1}{4} J(1) + C_R \sup_{z \in \mathbb{Z}^d} \|w_n\|_{L^q(C_z)}^{2-q} . $$

This implies

$$ Q_n(\sqrt{d}) \geq \sup_{z \in \mathbb{Z}^d} \|w_n\|_{L^q(C_z)}^q \geq \frac{J(1)}{(4C_R)^{\frac{q}{2-q}}} > 0 $$

and finally $\mu > 0$, which discards the vanishing case of the concentration-compactness method.

• Dichotomy. By definition of $Q_n$, there are sequences of centers $x_n \in \mathbb{R}^d$ and radii $\rho_n > 0$ going to infinity so that

$$ \lim_{n \to +\infty} \int_{B(x_n, \rho_n)} |w_n|^q dx = \mu. $$

Without loss of generality, by translating the functions $w_n$, we may assume $x_n = 0$. In addition, up to a non-displayed subsequence, we have that for all $\varepsilon > 0$, there is $n_0$ large enough so that, for all $n \ge n_0$, we have

$$ \int_{\rho_n < |x| < 2 \rho_n} |w_n|^q dx < \varepsilon \quad \text{and} \quad \left| 1 - \mu - \int_{|x| > 2 \rho_n} |w_n|^q dx \right| < \varepsilon. $$

We set

$$ \left\{ \begin{array}{ll} w_n^{(1)} & := w_n \mathbb{1}_{\{|x| \le \rho_n\}}, \\ w_n^{(2)} & := w_n \mathbb{1}_{\{\rho_n \le |x| \le 2 \rho_n\}}, \\ w_n^{(3)} & := w_n \mathbb{1}_{\{|x| > 2 \rho_n\}}. \end{array} \right. $$

Introducing $E(w_1, w_2) := \langle w_1, K * w_2 \rangle$ and $\mathcal{E}(w) := E(w, w)$, we have

$$ \begin{aligned} \mathcal{E}(w_n) &= \mathcal{E}(w_n^{(1)}) + \mathcal{E}(w_n^{(2)}) + \mathcal{E}(w_n^{(3)}) \\ &\quad + 2 \operatorname{Re} (E(w_n^{(1)}, w_n^{(2)}) + E(w_n^{(1)}, w_n^{(3)}) + E(w_n^{(2)}, w_n^{(3)})) \\ &\le J(\mu) + J(\varepsilon) + J(1 - \mu + \varepsilon) \\ &\quad + 2 \operatorname{Re} (E(w_n^{(1)}, w_n^{(2)}) + E(w_n^{(1)}, w_n^{(3)}) + E(w_n^{(2)}, w_n^{(3)})) . \end{aligned} $$

From (3.9), and the fact that $\|w_n^{(2)}\|_q \le \varepsilon^{1/q}$, we get that

$$ E(w_n^{(1)}, w_n^{(2)}) \le C \mu \varepsilon^{1/q} \quad \text{and} \quad E(w_n^{(2)}, w_n^{(3)}) \le C (1 - \mu) \varepsilon^{1/q}. $$

Finally, we have

$$ \begin{aligned} |E(w_n^{(1)}, w_n^{(3)})| &\le \int_{|x| \le \rho_n} \int_{|y| \ge 2\rho_n} |\langle w_n(x), K(x-y)w_n(y) \rangle_{\mathbb{C}^N}| dx dy \\ &\le \iint_{\mathbb{R}^d \times \mathbb{R}^d} |w_n(x)| |w_n(y)| (K \mathbb{1}_{\mathcal{B}_{\rho_n}^c})(x-y) dx dy \\ &\le \|K \mathbb{1}_{\mathcal{B}_{\rho_n}^c}\|_r C \varepsilon \end{aligned} $$
for $n$ large enough, where in the last line we used Young's inequality, and the fact that $\rho_n \to +\infty$. Thanks to these facts, we can conclude that

$$J(1) \le J(\mu) + J(1 - \mu + \varepsilon) + J(\varepsilon) + C \varepsilon^{1/q}.$$

In the limit as $\varepsilon \to 0$, we obtain $J(1) \le J(\mu) + J(1 - \mu)$, which contradicts (3.11) if $\mu \ne 1$. So $\mu = 1$, which discards the *dichotomy* case of the concentration-compactness method.

• **Convergence for tight sequences.** At this point, we proved that for all $\varepsilon > 0$ there is $\rho > 0$ and $n_0$ large enough so that, for all $n > n_0$, and after appropriate translations and subsequences,

$$(3.12) \qquad \|1_{B_{\rho}^c} w_n\|_q \le \varepsilon.$$

In other words, the sequence $(w_n)_{n \in \mathbb{N}}$ is tight in $L^q(\mathbb{R}^d, \mathbb{C}^N)$. The sequence $(w_n)_{n \in \mathbb{N}}$ is bounded in the reflexive Banach space $L^q(\mathbb{R}^d, \mathbb{C}^N)$. Hence, up to a non-displayed sub-sequence, $(w_n)_{n \in \mathbb{N}}$ converges weakly to some $w \in L^q(\mathbb{R}^d, \mathbb{C}^N)$, and we have $\|w\|_q \le 1$.

Let us prove that $\mathcal{E}(w) = J(1)$. Let $\varepsilon > 0$, and let $\rho > 0$ be large enough so that (3.12) holds. In particular, by Hardy-Littlewood-Sobolev, we have

$$\left| \langle w_n \mathbb{1}_{B_{\rho}^c}, K * w_n \rangle \right| \le C \|w_n\|_q \|w_n \mathbb{1}_{B_{\rho}^c}\|_q \le C \varepsilon,$$

and we have a similar inequality with $w$ instead of $w_n$. On the other hand, we have

$$\langle w_n \mathbb{1}_{B_{\rho}}, K * w_n \rangle = \langle w_n, T w_n \rangle_{L^q, L^{q'}},$$

where $T$ is the operator from $L^q(\mathbb{R}^d)$ to $L^{q'}(\mathbb{R}^d)$ with kernel $T(x, y) = \mathbb{1}_{B_{\rho}}(x) K(x - y)$. The operator $T: L^q(\mathbb{R}^d, \mathbb{C}^N) \to L^{q'}(\mathbb{R}^d, \mathbb{C}^N)$ is bounded. We claim that $T$ is a compact operator. In the Dirac case (ii), this comes from the fact that $K * w_n \in W^{1,q}$ with $\|K * w_n\|_{W^{1,q}(\mathbb{R}^d, \mathbb{C}^N)} \le C \|w_n\|_q$ together with the Rellich-Kondrachov compact embedding theorem. In the case (i), where $K \in L^r(\mathbb{R}^d)$, setting $\tau_h f(x) := f(x - h)$, we have

$$\|\tau_h(K * w) - K * w\|_{L^{q'}(B_{\rho})} = \|(\tau_h K - K) * w\|_{q'} \le \|\tau_h K - K\|_r \|w\|_q.$$

Since $K \in L^r(\mathbb{R}^d)$, we have $\|\tau_h K - K\|_r \to 0$ as $h \to 0$, and we conclude with the Kolmogorov-Riesz-Fréchet theorem (see for instance [11, Theorem 4.26]).

As a consequence, $(Tw_n)_{n \in \mathbb{N}}$ converges strongly to $Tw$ in $L^{q'}(\mathbb{R}^d, \mathbb{C}^N)$. In particular, we obtain that

$$\lim_{n \to +\infty} \langle w_n \mathbb{1}_{B_{\rho}}, K * w_n \rangle = \langle w \mathbb{1}_{B_{\rho}}, K * w \rangle.$$

Gathering the two inequalities gives

$$\begin{aligned} & |\langle w_n, K * w_n \rangle - \langle w, K * w \rangle| \\ & \le |\langle w_n \mathbb{1}_{B_{\rho}}, K * w_n \rangle - \langle w \mathbb{1}_{B_{\rho}}, K * w \rangle| + |\langle w_n \mathbb{1}_{B_{\rho}^c}, K * w_n \rangle| + |\langle w \mathbb{1}_{B_{\rho}^c}, K * w \rangle| \\ & \le |\langle w_n \mathbb{1}_{B_{\rho}}, K * w_n \rangle - \langle w \mathbb{1}_{B_{\rho}}, K * w \rangle| + 2C\varepsilon. \end{aligned}$$

Sending first $n$ to $+\infty$, and then $\varepsilon$ to 0 shows that $\langle w, K * w \rangle = J(1)$. Finally, since $\|w\|_q \le 1$, by (3.10) we deduce that $\|w\|_q = 1$. This proves that $(w_n)_{n \in \mathbb{N}}$ converges strongly to $w$ in $L^q(\mathbb{R}^d, \mathbb{C}^N)$ and that $w$ is an optimizer. ■
It is an open question to decide whether $T$ is compact or not under the condition (3.8).

**Proof of Theorem 3.1.** In the setting of Theorem 3.1, we take $K = R_0(\lambda)$. We have $J(1) > 0$ in this case, as noticed in Remark 3.6, so by Lemma 3.3, the problem (3.3) admits maximizers. By standard arguments, optimizers satisfy the Euler-Lagrange equation (3.4).

The fact that $\lambda \mapsto \mathcal{N}(\lambda, p)$ is strictly increasing comes from the fact that $\lambda \mapsto R_0(\lambda)$ is operator strictly increasing: for instance, we have $\partial_\lambda R_0(\lambda) = (R_0(\lambda))^2 > 0$. Let us prove the continuity. Let $-m < \lambda' < \lambda < m$, and let $w_\lambda$ be the optimizer for $\mathcal{N}(\lambda, p)$. Using that $\mathcal{N}(\cdot)$ is strictly increasing and the resolvent identity

$$R_0(\lambda') = R_0(\lambda) - (\lambda - \lambda') R_0(\lambda') R_0(\lambda),$$

we obtain

$$0 < \mathcal{N}(\lambda, p) - \mathcal{N}(\lambda', p) \le (\lambda - \lambda') \langle w_\lambda, R_0(\lambda') R_0(\lambda) w_\lambda \rangle.$$

Using that $R_0$ is a bounded operator from $L^q(\mathbb{R}^d)$ to $L^{q'}(\mathbb{R}^d)$, and from $L^{q'}(\mathbb{R}^d)$ into itself, with uniform bounds in a neighborhood of $\lambda$, we deduce that there is $C > 0$ so that

$$|\langle w_\lambda, R_0(\lambda') R_0(\lambda) w_\lambda \rangle| \le C \|w_\lambda\|_q^2 = C,$$

This proves that $\mathcal{N}(\cdot, p)$ is locally Lipschitz, hence continuous.

We now prove the bounds on $\lim_{\lambda \to \pm m} \mathcal{N}(\lambda, p)$. To prove that $\lim_{\lambda \to m} \mathcal{N}(\lambda, p) = +\infty$, we go back to (3.7) and take a function $f = L^{-d/q} g(\cdot/L)$, where $g$ is an arbitrary test function that is normalized in $L^q(\mathbb{R}^d)$. This gives

$$\mathcal{N}(\lambda, p) \ge L^{-2d/q} (m+\lambda) \left\langle g(\cdot/L), (-\Delta + m^2 - \lambda^2)^{-1} g(\cdot/L) \right\rangle.$$

We bound the resolvent as

$$(-\Delta + m^2 - \lambda^2)^{-1} \ge (m^2 - \lambda^2)^{-1} (1 + (m^2 - \lambda^2)^{-1} \Delta)$$

and change variables to obtain

$$\mathcal{N}(\lambda, p) \ge L^{d(1-\frac{2}{q})} (m-\lambda)^{-1} \left( \|g\|_2^2 - L^{-2} (m^2 - \lambda^2)^{-1} \|\nabla g\|_2^2 \right).$$

Since $1 - \frac{2}{q} = p$, we may take $L = (m - \lambda)^{-\alpha}$ for any $\alpha \in (1/2, p/d)$ and conclude that $\lim_{\lambda \to m} \mathcal{N}(\lambda, p) = +\infty$.

Finally, to prove that $\lim_{\lambda \to -m} \mathcal{N}(\lambda, p) > 0$, we claim that

(3.13) there exists a function $w \in L^2 \cap L^q(\mathbb{R}^d, \mathbb{C}^N)$ such that $\|w\|_q = 1$ and $Pw = w$,

where $P := \mathbb{1}_{m<\mathcal{P}_m<2m}$ is the spectral projection of the free Dirac operator onto $(m, 2m)$. This would give

$$\mathcal{N}(\lambda, p) \ge \langle w, R_0(\lambda) w \rangle = \left\langle w, \frac{P}{\mathcal{P}_m - \lambda} w \right\rangle + \left\langle w, \frac{P^{\perp}}{\mathcal{P}_m - \lambda} w \right\rangle.$$
The second term is null since $P^{\perp} w = 0$. For the first term, we have $m < \mathcal{D}_m < 2m$ on the range of $P$, and in particular $P(\mathcal{D}_m - \lambda)^{-1} P \ge P(2m - \lambda)^{-1} P$, hence $\mathcal{N}(\lambda, p) \ge (2m - \lambda)^{-1} \|w\|_2$. Taking $\lambda \to -m$ shows that $\lim_{\lambda \to -m} \mathcal{N}(\lambda, p) \ge (3m)^{-1} \|w\|_2 > 0$.

It remains to prove (3.13). Recall that $\mathcal{D}_m = \mathcal{F}M(k)\mathcal{F}^*$, where $\mathcal{F}$ denotes the Fourier transform and $M(k)$ is the $d \times d$ matrix $M(k) := \alpha \cdot k + m\beta$, which satisfies $M(k) = M(k)^*$, $M(k)^2 = (|k|^2 + m^2) \mathbb{I}_d$, and $\sigma(M(k)) = \{\pm (|k|^2 + m^2)^{1/2}\}$. Let $v \mapsto v(k)$ be a smooth family of spinors from some open ball $B(k=0, \varepsilon)$ to $\mathbb{C}^d$, with $0 < \varepsilon < m$, so that $M(k)v(k) = (|k|^2 + m^2)^{1/2} v(k)$. To construct such a local family of spinors, one can consider $v_0$ a normalized eigenfunction of $M(k=0)$, and set,

$$v(k) := \frac{P(k)v_0}{\|P(k)v_0\|_2}, \quad P(k) := 1(M(k) > 0).$$

Since $k \mapsto P(k)$ is smooth locally around 0 ($P(k)$ can be written as a Cauchy integral $P(k) = (2i\pi)^{-1} \oint_{\mathcal{C}} (z - M(k))^{-1} dz$ with a contour enclosing $m$), so is $k \mapsto v(k)$. Let also $\chi(k) : \mathbb{R}^d \to \mathbb{R}^+$ be a non null smooth compactly supported function, with $\chi(k) = 0$ for $|k| > \varepsilon$. We consider the function

$$w := \frac{\tilde{w}}{\|\tilde{w}\|_q} \quad \text{with} \quad \tilde{w} := \mathcal{F}(\chi(k)v(k)).$$

By construction, we have $\tilde{w} \ne 0$, and since $\tilde{w}$ has a Fourier transform which is smooth and compactly supported, it belongs to the Schwartz class $S(\mathbb{R}^d, \mathbb{C}^N)$. Finally, since on the support of $\chi$, we have $M(k)v(k) = (|k|^2 + m^2)^{1/2} v(k)$ with $m < (|k|^2 + m^2)^{1/2} < \sqrt{2}m^2$, we deduce that

$$P \tilde{w} = \mathcal{F}(1_{m<M(k)<2m} \chi(k)v(k)) = \mathcal{F}(\chi(k)v(k)) = \tilde{w},$$

which concludes the proof of (3.13). ■

### 3.3. Regularity of the solutions of the non-linear Dirac equation

Under Condition (2.1), solutions of (3.5) with $\Psi \in L^2(\mathbb{R}^d, \mathbb{C}^N)$ are in $\text{Dom}(\mathcal{D}_m - V) = \text{H}^1(\mathbb{R}^d, \mathbb{C}^N)$. Let us consider the other cases of Proposition 2.1. If $d=1$ and $1<p \le 2$, any optimal function for (3.3) obtained in Theorem 3.1 gives rise to a solution $\Psi \in \text{W}^{1,q}(\mathbb{R}, \mathbb{C}^2)$ of (3.5) with $q=2p/(p+1)$. We conclude that $\Psi$ is continuous. If $p=d=2$ and $q=4/3$, the corresponding solution $\Psi$ of (3.5) is in $\text{W}^{1,q}(\mathbb{R}, \mathbb{C}^2) \hookrightarrow \text{H}^{1/2}(\mathbb{R}, \mathbb{C}^2)$, hence $V|\Psi|^2 = |\Psi|^{2p/(p-1)}$ is integrable and $\Psi \in \text{Dom}(\mathcal{D}_m - V)$ of the distinguished extension of Proposition 2.1 but we do not know whether $\Psi \in \text{H}^1(\mathbb{R}^2, \mathbb{C}^2)$ or not.

In dimension $d=1$, an explicit expression of the solutions of (3.5) such that

$$\lim_{x \to \pm\infty} \Psi(x) = (0, 0)^{\top}$$

is given in Theorem 1.3. In the case $p=d=2$, it is unclear how to obtain $\Psi \in \text{H}^1(\mathbb{R}^d, \mathbb{C}^N)$ by general arguments, as pointed out in [6]. However, any solution to (3.5) (and not only the ones found in Theorem 3.1) have additional regularity properties under Condition (2.1).
**Proposition 3.7.** Let $\lambda \in [-m, m)$ and either $p \ge d$ if $d \ge 3$, or $p > d$ in dimension $d = 1$ and $2$. If $\Psi \in H^1(\mathbb{R}^d, \mathbb{C}^N)$ solves (3.5), then $\Psi \in C^\infty(\mathbb{R}^d, \mathbb{C}^N)$.

*Proof.* Let us first prove that $\Psi \in L^\infty(\mathbb{R}^d, \mathbb{C}^N)$ with a usual bootstrap argument. If $\Psi \in L^q(\mathbb{R}^d, \mathbb{C}^N)$, then $|\Psi|^{\frac{2}{p-1}} \Psi \in L^{\frac{p-1}{p+1}} q(\mathbb{R}^d, \mathbb{C}^N)$. Also, if $q > 2 \frac{p+1}{p-1}$, then $2 < \frac{p-1}{p+1} q < q$, so if $\Psi \in L^2(\mathbb{R}^d, \mathbb{C}^N) \cap L^q(\mathbb{R}^d, \mathbb{C}^N)$, then $\lambda \Psi + |\Psi|^{\frac{2}{p-1}} \Psi \in L^{\frac{p-1}{p+1}} q(\mathbb{R}^d, \mathbb{C}^N)$. In particular, $\mathcal{D}_m \Psi \in L^{\frac{p-1}{p+1}} q(\mathbb{R}^d, \mathbb{C}^N)$, hence $\Psi \in W^{1, \frac{p-1}{p+1}} q(\mathbb{R}^d, \mathbb{C}^N) \hookrightarrow L^{\tilde{q}}(\mathbb{R}^d, \mathbb{C}^N)$, with $\tilde{q} = +\infty$ if $\frac{p-1}{p+1} q > d$ and

$$ \frac{1}{\tilde{q}} = \frac{p+1}{p-1} \frac{1}{q} - \frac{1}{d} $$

otherwise. As a first step of an iteration scheme, we proved that if $\Psi \in L^q(\mathbb{R}^d, \mathbb{C}^N)$, then $\Psi \in L^{\tilde{q}}(\mathbb{R}^d, \mathbb{C}^N)$ as well. For the initialization, we note that $H^1(\mathbb{R}^d, \mathbb{C}^N) \hookrightarrow L^q(\mathbb{R}^d, \mathbb{C}^N)$ for all $q$ such that $2 \le q \le \frac{2d}{d-2} =: 2^*$ if $d \ge 3$ and $2 \le q < +\infty =: 2^*$ if $d = 2$. Hence with $2 \frac{p+1}{p-1} < 2^*$, there is $q_0 > 2 \frac{p+1}{p-1}$ so that $\Psi \in L^2(\mathbb{R}^d, \mathbb{C}^N) \cap L^{q_0}(\mathbb{R}^d, \mathbb{C}^N)$. The map $F: x \mapsto \frac{p+1}{p-1} x - \frac{1}{d}$ satisfies $F(x) < x$ for $x \in [0, x^*]$ with $x^* = \frac{p-1}{2d} < \frac{1}{2} \frac{p-1}{p+1}$. We easily deduce that there is $n \in \mathbb{N}$ so that $F^{(n)}(\frac{1}{q_0}) < 0$, which proves $\Psi \in L^\infty(\mathbb{R}^d, \mathbb{C}^N)$ as wanted.

Since $\mathcal{D}_m \Psi \in L^\infty(\mathbb{R}^d, \mathbb{C}^N)$, we have $\Psi \in W^{1, \infty}(\mathbb{R}^d, \mathbb{C}^N) \hookrightarrow C^{0, \alpha}(\mathbb{R}^d, \mathbb{C}^N)$ for all $0 \le \alpha < 1$, by bootstrapping again, we obtain $\Psi \in C^\infty(\mathbb{R}^d, \mathbb{C}^N)$. ■

## 4. Lieb-Thirring inequality

This section contains the proof of Theorem 1.4. We closely follow the original proof by Lieb and Thirring [51, 53] (see also [52]). This is possible since we are assuming $V \ge 0$. In the general case where $V$ has no sign, some results can be found in the works of Cuenin [17], and Frank-Simon [40], where the authors control the Riesz-mean

$$ \sum_k \mathrm{dist}(\lambda_k, \sigma(\mathcal{D}_m - V))^\gamma, $$

that is, the distance to the whole spectrum. Actually, without assuming a sign on $V$, one cannot expect to control the sums in (1.13), since, for $V \le 0$ small, the eigenvalues of $\mathcal{D}_m - V$ emerge from the bottom essential spectrum (hence have a distance of order $2m > 0$ to the upper essential spectrum). Here, since $V$ is nonnegative, the eigenvalues emerge from the upper essential spectrum as the strength of the potential increases.

*Proof of Theorem 1.4.* It is sufficient to prove the result for $V$ bounded and compactly supported. By the *Birman-Schwinger principle* introduced in Section 2.2, we know that $\lambda$ is an eigenvalue for $\mathcal{D}_m - V$ acting on $\mathbb{C}^N$ valued spinors if and only if 1 is an eigenvalue of $K_V(\lambda)$ defined by (2.3): see Proposition 2.4. We also proved that $\lambda \mapsto K_V(\lambda)$ is operator increasing. In particular, if we set

$$ N_e(V) := \text{number of eigenvalues of } \mathcal{D}_m - V \text{ in } [-m, m-e] $$

and

$$ B_e(V) := \text{number of eigenvalues of } K_V(m-e) \text{ which are greater or equal than 1,} $$
then we have $N_e(V) \le B_e(V)$. We have equality if the highest eigenvalues of $K_V(\lambda)$ gets strictly smaller than 1 as $\lambda \to -m$. This happens for instance if $\|V\|_p \le \alpha_*(p)$.

With $R_0$ defined by (2.2), using the operator inequality

$$R_0(\lambda) \le 1_{\mathbb{C}^N} (\sqrt{-\Delta + m^2} - \lambda)^{-1},$$

we can estimate $B_e(V)$ by $N B_e^{\text{pr}}(V)$, where $B_e^{\text{pr}}(V)$ is the number of eigenvalues above 1 of the pseudo-relativistic Birman-Schwinger operator

$$K_V^{\text{pr}}(m-e) := \sqrt{V} (\sqrt{-\Delta + m^2} - m + e)^{-1} \sqrt{V}.$$

In addition, with the definition

$$N_e^{\text{pr}}(V) := \text{number of eigenvalues of } (\sqrt{-\Delta + m^2} - m) - V \text{ less or equal than } -e,$$

the usual Birman-Schwinger principle shows that $B_e^{\text{pr}}(V) = N_e^{\text{pr}}(V)$. To sum up, we have

$$(4.1) \qquad N_e(V) \le B_e(V) \le N B_e^{\text{pr}}(V) = N N_e^{\text{pr}}(V).$$

The operator $\sqrt{-\Delta + m^2} - m$ is sometimes called the Chandrasekhar (or pseudo-relativistic) kinetic energy operator. It is a positive operator, $\sqrt{-\Delta + m^2} - m - V$ is bounded from below, and the min-max formula applies. We can now repeat the usual arguments of Lieb and Thirring for the pseudo-relativistic operator.

First, for $\gamma > 0$, the cake-layer representation gives

$$(4.2) \qquad \sum_{k \ge 1} e_k^\gamma = \gamma \int_0^{2m} e^{\gamma-1} N_e(V) \, de \le \gamma N \int_0^{2m} e^{\gamma-1} B_e^{\text{pr}}(V) \, de.$$

Note that for the pseudo-relativistic model, if $-e_1^{\text{pr}} \le -e_2^{\text{pr}} \le \dots < 0$ are the negative eigenvalues of $(\sqrt{-\Delta + m^2} - m) - V$, we have

$$\sum_{k \ge 1} (e_k^{\text{pr}})^\gamma = \gamma \int_0^\infty e^{\gamma-1} N_e^{\text{pr}}(V) \, de = \gamma \int_0^\infty e^{\gamma-1} B_e^{\text{pr}}(V) \, de,$$

and the integral runs over $e \in \mathbb{R}^+$ instead of $e \in (0, 2m)$. Actually, the previous two inequalities together with (4.1) show that

$$\sum_{k \ge 1} e_k^\gamma \le N \sum_{k \ge 1} (e_k^{\text{pr}})^\gamma.$$

In other words, the Riesz-mean of the eigenvalues increases when one replaces the Dirac operator by the pseudo-relativistic one (up to the $N$ factor). Lieb-Thirring inequalities for the last sum have been derived by Daubechies in [19] (and used, e.g., in [50]). In what follows, we derive another inequality specifically for the Dirac operator. We use in particular the fact that the integral in (4.2) only runs for $e$ in the bounded interval $(0, 2m)$ instead of $\mathbb{R}^+$.
• Bound for $B_e^{\text{pr}}(V)$. Assume $V \in L^p(\mathbb{R}^d)$ with $d < p$. The number of eigenvalues above 1 of $K_V^{\text{pr}}(m-e)$ is bounded from above by $\|K_V^{\text{pr}}(m-e)\|_{\mathfrak{S}_p}^p$. We estimate this norm using the Kato-Simon-Seiler inequality (see [65, Theorem 4.2]). Using a decomposition similar to the one in the proof of Lemma 2.3, we obtain

$$B_e^{\text{pr}}(V) \le \|K_V^{\text{pr}}(m-e)\|_{\mathfrak{S}_p}^p \le C_p \|g_{m,e}\|_p^p \|V\|_p^p,$$

where we introduced the function

$$g_{m,e}(k) := \left(\sqrt{k^2 + m^2} - (m - e)\right)^{-1}.$$

Note that $g_{m,e} \in L^p(\mathbb{R}^d)$ since $p > d$, and

$$\|g_{m,e}\|_p^p = \int_{\mathbb{R}^d} \frac{\mathrm{d}k}{(\sqrt{k^2 + m^2} - m + e)^p} = |\mathbb{S}^{d-1}| \int_0^\infty \frac{r^{d-1} \mathrm{d}r}{(\sqrt{r^2 + m^2} - m + e)^p}.$$

To estimate this norm, we make the change of variable $X = \frac{1}{e}(\sqrt{r^2 + m^2} - m)$, so that $r = \sqrt{(eX + m)^2 - m^2} = \sqrt{eX(eX + 2m)}$. We obtain

$$\|g_{m,e}\|_p^p = \frac{|\mathbb{S}^{d-1}|}{e^{p-\frac{d}{2}}} \int_0^\infty \frac{(X(eX + 2m))^{\frac{d}{2}-1} (eX + m) \mathrm{d}X}{(X+1)^p}.$$

The last integral is an increasing function of $e$ (and has a finite value as $e \to 0$ by the monotone convergence theorem). Since $e \in (0, 2m)$, we can bound this integral by its value at $e = 2m$. We deduce that there is a constant $C_{p,d}$ such that

$$(4.3) \qquad B_e^{\text{pr}}(V) \le C_{p,d} \|V\|_p^p \frac{m^{d/2}}{e^{p-d/2}}.$$

• Proof of the Lieb-Thirring estimate. We now follow [51–53]. The min-max principle for the pseudo-relativistic operator shows that its eigenvalues are decreasing when $V$ increases. Since $V \le [V - e/2]_+ + e/2$, we may bound

$$\begin{aligned} B_e^{\text{pr}}(V) &= N_e^{\text{pr}}(V) \le N_e^{\text{pr}}([V - e/2]_+ + e/2) \\ &= N_{e/2}^{\text{pr}}([V - e/2]_+) = B_{e/2}^{\text{pr}}([V - e/2]_+). \end{aligned}$$

For any $p > d$, we can apply the bound in (4.3) to estimate $B_{e/2}^{\text{pr}}([V - e_2]_+)$. Inserting this estimate into (4.2), we get

$$\begin{aligned} \sum_{k \ge 1} e_k^\gamma \le N C_{p,d} \gamma m^{\frac{d}{2}} \int_0^{2m} (e/2)^{\gamma-1+\frac{d}{2}-p} \|[V-e/2]_+\|_p^p \mathrm{d}e \\ = C_{\gamma,d,p} m^{\frac{d}{2}} \int_{\mathbb{R}^d} \int_0^{2m} e^{\gamma-1+\frac{d}{2}-p} [V(x)-e/2]_+^p \mathrm{d}x \mathrm{d}y \\ \le C_{\gamma,d,p} m^{\frac{d}{2}} \int_{\mathbb{R}^d} V^{\gamma+\frac{d}{2}}(x) \int_0^{s^*(x)} s^{\gamma-1+\frac{d}{2}-p} (1-s)_+^p \mathrm{d}s \mathrm{d}x, \end{aligned}$$
where $s^*(x) := \min\{m/V(x), 1\}$, with the convention that $s^*(x) = 1$ if $V(x) = 0$. The second integral converges whenever $p < \gamma + d/2$. We can simply use the bound $(1-s)^p \le 1$ in the last integral, and finally obtain

$$ (4.4) \qquad \sum_{k \ge 1} e_k^\gamma \le L_{\gamma, d, p} m^{\frac{d}{2}} \int_{\mathbb{R}^d} V_m^{\gamma + \frac{d}{2} - p} V^p dx \quad \text{with} \quad V_m := \min\{m, V\}. $$

This inequality is valid for all $d < p < \gamma + d/2$. Note that $C_{\gamma,d,p}$ stays bounded in the limit as $p \to \gamma + d/2$, so a similar inequality also holds if $p = \gamma + d/2$. ■

**Remark 4.1.** The result of Theorem 1.4 can be extended to the case of a potential $V \in L^p(\mathbb{R}^d, \mathbb{R}^+) + L^{\gamma+d/2}(\mathbb{R}^d, \mathbb{R}^+)$ by noticing that the right-hand side of (4.4) is continuous for $V$ in this space.

# 5. Explicit computations

## 5.1. The case $d = 1$: proof of Theorem 1.3

In this section, we prove the uniqueness and the symmetry up to translations of the solution of the nonlinear Dirac equation (3.5). We also compute the map $\alpha_D(\lambda, p)$.

*Proof of Theorem 1.3.* In the one-dimensional case, Equation (3.5) can be rewritten for the components of $\Psi := (\varphi, \chi)^\top$ as

$$ (5.1) \qquad \left\{ \begin{array}{l} \varphi' = - \left( \lambda + m + (|\chi|^2 + |\varphi|^2)^{\frac{1}{p-1}} \right) \chi, \\ \chi' = \left( \lambda - m + (|\chi|^2 + |\varphi|^2)^{\frac{1}{p-1}} \right) \varphi. \end{array} \right. $$

The corresponding potential is $V = (|\chi|^2 + |\varphi|^2)^{\frac{1}{p-1}}$. This system conserves

$$ \begin{aligned} H(\varphi, \chi) &:= m (|\chi|^2 - |\varphi|^2) + \lambda (|\chi|^2 + |\varphi|^2) + \frac{p-1}{p} (|\chi|^2 + |\varphi|^2)^{\frac{p}{p-1}}, \\ G(\varphi, \chi) &:= \bar{\chi} \varphi - \bar{\varphi} \chi. \end{aligned} $$

Since we are looking for solutions vanishing at $\pm\infty$, they satisfy $H(\varphi(x), \chi(x)) = 0$, $G(\varphi(x), \chi(x)) = 0$ for all $x \in \mathbb{R}$. This second condition shows that solutions can be chosen real valued. For real valued variables in the $(\varphi, \chi)$-plane, the level set $H(\varphi, \chi) = 0$ has the shape of an infinity sign. Among real valued functions, uniqueness up to translations follows from the phase plane analysis. We can choose the unique solution with $\chi(0) = 0$, $\varphi(0) > 0$, given. For this solution, $\varphi$ is even and $\chi$ is odd and positive on $\mathbb{R}^+$. Hence symmetry and uniqueness, up to translations and multiplication by a phase, are granted by elementary considerations. Next, we have

$$ V' = \frac{1}{p-1} \frac{(\chi^2 + \varphi^2)'}{(\chi^2 + \varphi^2)^{\frac{p}{p-1}}} \quad \text{with} \quad (\chi^2 + \varphi^2)' = -4m\chi\varphi, $$
which proves that $V$ is increasing in the quadrant $\{\chi < 0, \varphi > 0\}$ and decreasing in the quadrant $\{\chi > 0, \varphi > 0\}$. Hence $V$ is even and decreasing on $\mathbb{R}^+$, while on $\mathbb{R}^+$ both $\chi$ and $\varphi$ are positive valued.

Now let us compute $\|V\|_p$. It is enough to do the computation on $\mathbb{R}^+$. First, the equation $H(\varphi, \chi) = 0$ can be rewritten as

$$2 m \varphi^2 = (m + \lambda) V^{p-1} + \frac{p-1}{p} V^p,$$

and so

$$\varphi = \sqrt{\frac{1}{2m} V^{p-1} \left(m + \lambda + \frac{p-1}{p} V\right)}.$$

Next, from the equation $V^{p-1} = \chi^2 + \varphi^2$, we deduce that

$$\chi = \sqrt{V^{p-1} - \varphi^2} = \sqrt{\frac{1}{2m} V^{p-1} \left(m - \lambda - \frac{p-1}{p} V\right)}.$$

Finally, we have

$$(p-1) V^{p-2} V' = (V^{p-1})' = (\chi^2 + \varphi^2)' = -4 m \chi \varphi.$$

Collecting the three last equalities shows that $V$ solves the autonomous differential equation

$$V' = - \frac{2}{p-1} V \sqrt{\left(m - \lambda - \frac{p-1}{p} V\right) \left(m + \lambda + \frac{p-1}{p} V\right)}.$$

At $x = 0$, we have $V'(0) = 0$, which implies

$$V(0) = \frac{p}{p-1} (m - \lambda).$$

• Subcritical regime $\lambda > -m$. The function

$$Z(x) := \frac{p-1}{p(m+\lambda)} V \left( \frac{p-1}{2(m+\lambda)} x \right)$$

satisfies

$$(5.2) \qquad Z' = -Z \sqrt{(z_0 - Z)(1+Z)}, \quad Z(0) = z_0 = \frac{m-\lambda}{m+\lambda}.$$

One can directly check that the solution of (5.2) is

$$Z(x) = \frac{2 z_0}{(1 + z_0) \cosh \left( \sqrt{z_0} x \right) + 1 - z_0}.$$

This gives (1.10). The $L^p(\mathbb{R})$ norm of $V$ is computed as

$$\|V\|_p^p = \frac{p^p (m + \lambda)^{p-1}}{2 (p-1)^{p-1}} \|Z\|_p^p.$$
Using that $Z$ is even, monotone decreasing on $\mathbb{R}^+$, with the change of variable $z = Z(x)$ and $t = z/z_0$, we obtain, using (5.2),

$$ \begin{aligned} \|Z\|_p^p &= 2 \int_0^{+\infty} Z^p(x) \, dx \\ &= 2 \int_0^{z_0} \frac{z^{p-1}}{\sqrt{(z_0-z)(1+z)}} \, dz = 2 z_0^{p-\frac{1}{2}} \int_0^1 \frac{t^{p-1}}{\sqrt{(1-t)(1-(-z_0)t)}} \, dt \\ &= 2 z_0^{p-\frac{1}{2}} B\left(\frac{1}{2}, p\right) {}_2F_1\left(\frac{1}{2}, p; p+\frac{1}{2}; -z_0\right). \end{aligned} $$

See [1, 15.3.1 p. 558] for the last equality. This completes the computation of $\alpha_D(\lambda, p)$. By taking the limit as $p \to 1_+$, we obtain $\alpha_D(\lambda, 1) = \arccos(\lambda/m)$.

• Critical case $\lambda = -m$. The function

$$ Z(x) := \frac{p-1}{2mp} V\left(\frac{p-1}{2m} x\right) $$

solves

$$ Z' = -2 Z^{3/2} \sqrt{1-Z}, \quad Z(0) = 1 $$

on $\mathbb{R}^+$. The solution is

$$ \forall x \in \mathbb{R}, \quad Z(x) = \frac{1}{1+x^2}. $$

This gives (1.11), and the expression of $\alpha_*(p)$ follows from

$$ \|V\|_p^p = \frac{p^p (2m)^{p-1}}{(p-1)^{p-1}} \|Z\|_p^p $$

with

$$ \|Z\|_p^p = \int_{\mathbb{R}} \frac{dx}{(1+x^2)^p} = B\left(\frac{1}{2}, p-\frac{1}{2}\right) $$

according to [42, 8.380.3 p. 917]. This concludes the proof of Theorem 1.3.

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

Notice that $\lim_{z_0 \to +\infty} \sqrt{z_0} B\left(\frac{1}{2}, p\right) {}_2F_1\left(\frac{1}{2}, p; p+\frac{1}{2}; -z_0\right) = B\left(\frac{1}{2}, p-\frac{1}{2}\right)$, so that $\lim_{\lambda \to (-1)_+} \alpha_D(\lambda, p) = \alpha_*(p)$.

## 5.2. The radial case in dimension $d=2$

We now provide some numerical simulations to obtain upper bounds for the maps $\Lambda_D(\alpha, p)$.

First, we restrict the minimization problem (1.5) to radial potentials, that is, we compute

$$ \Lambda_D^{\text{rad}}(\alpha, p) := \inf \left\{ \lambda_D(V) : V \in L^p(\mathbb{R}^d, \mathbb{R}^+), V \text{ radial and } \|V\|_p = \alpha \right\}. $$

Below in Appendix B, we provide some numerical evidences that the optimal potentials are radial. We abusively write $V(x) = V(r)$ with $r = |x|$, $x \in \mathbb{R}^2$, use polar coordinates $(x, y) = (r \cos \theta, r \sin \theta)$, and write

$$ \partial_x = \cos \theta \partial_r - \frac{1}{r} \sin \theta \partial_\theta \quad \text{and} \quad \partial_y = \sin \theta \partial_r + \frac{1}{r} \cos \theta \partial_\theta. $$
In these coordinates, the Dirac operator becomes

$$ \mathcal{D}_m = \begin{pmatrix} m & -i \partial_x - \partial_y \\ -i \partial_x + \partial_y & -m \end{pmatrix} = \begin{pmatrix} m & e^{-i\theta} \left(-i\partial_r - \frac{1}{r}\partial_\theta\right) \\ e^{i\theta} \left(-i\partial_r + \frac{1}{r}\partial_\theta\right) & -m \end{pmatrix}. $$

This suggests to decompose a spinor $\Psi$ in Fourier modes with the convention

$$ \Psi(r, \theta) = \sum_{n \in \mathbb{Z}} \begin{pmatrix} \varphi_n(r) e^{i n \theta} \\ i \chi_n(r) e^{i (n+1) \theta} \end{pmatrix}. $$

If $\Phi := (\mathcal{D}_m - V) \Psi$ with corresponding Fourier modes $((\tilde{\varphi}_n, \tilde{\chi}_n)^\top)_{n \in \mathbb{Z}}$, then we have

$$ \begin{pmatrix} \tilde{\varphi}_n \\ \tilde{\chi}_n \end{pmatrix} = (\mathcal{D}_m^{(n)} - V) \begin{pmatrix} \varphi_n \\ \chi_n \end{pmatrix} \quad \text{with} \quad \mathcal{D}_m^{(n)} = \begin{pmatrix} m & \partial_r + \frac{n+1}{r} \\ -\partial_r + \frac{n}{r} & -m \end{pmatrix}. $$

The operator $\mathcal{D}_m^{(n)}$ is self-adjoint in the Hilbert space $L^2(\mathbb{R}^+ \times (0, 2\pi), r dr d\theta)$ because $(\partial_r)^* = -\partial_r - \frac{1}{r}$. Let $\lambda_D^{(n)}(V)$ denote the lowest eigenvalue of $\mathcal{D}_m^{(n)} - V$ in the gap $(-m, m)$, and let

$$ \Lambda_D^{\text{rad},(n)}(\alpha, p) := \inf \left\{ \lambda_D^{(n)}(V) : V \in L^p(\mathbb{R}^2, \mathbb{R}^+), V \text{ radial and } \|V\|_p = \alpha \right\} $$

and

$$ \Lambda_D^{\text{rad}}(\alpha, p) := \inf_{n \in \mathbb{Z}} \Lambda_D^{\text{rad},(n)}(\alpha, p). $$

We have the estimates

$$ (5.3) \qquad \Lambda_D(\alpha, p) \le \Lambda_D^{\text{rad}}(\alpha, p) \le \Lambda_D^{\text{rad},(0)}(\alpha, p). $$

A wavefunction $\Psi(r, \theta) = (\varphi(r) e^{i n \theta}, i \chi(r) e^{i (n+1) \theta})^\top$ solves the non-linear Dirac equation (3.5) if and only if

$$ (5.4) \qquad \begin{cases} \varphi' - \frac{n}{r} \varphi = - \left( \lambda + m + (|\chi|^2 + |\varphi|^2)^{\frac{1}{p-1}} \right) \chi, \\ \chi' + \frac{n+1}{r} \chi = \left( \lambda - m + (|\chi|^2 + |\varphi|^2)^{\frac{1}{p-1}} \right) \varphi. \end{cases} $$

This system with $n=0$ is studied by W. Borrelli in [5]. It is an open question to decide whether $\Lambda_D^{\text{rad}}(\alpha, p)$ is attained by $\Lambda_D^{\text{rad},(n)}(\alpha, p)$ with $n=0$ or not, and if equality holds in (5.3) so that $\Lambda_D(\alpha, p) = \Lambda_D^{\text{rad},(0)}(\alpha, p)$. See Fig. 3 for some numerical results.

### 5.3. The radial case in dimension $d=3$

As in the two dimensional case, we restrict the minimization problem (1.5) to radially symmetric decreasing potentials. The corresponding Dirac operator decomposes as a direct sum in eigenspaces of the spin-orbit operator

$$ K = \beta (2 S \cdot L + 1) = \beta (J^2 - L^2 + 1/4), \quad \text{spec}(K) = \pm 1, \pm 2, \dots $$
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Radial case with $d=2$ and $m=1$. (Left) The function $p \mapsto \alpha_{\star}^{\text{rad},(n=0)}(p)$ is an upper bound for $\alpha_{\star}(p)$ and reaches its maximum for $p \approx 2.66$. (Right) The maps $\alpha \mapsto \Lambda_{D}^{\text{rad},(n=0)}(\alpha, p)$ for values of $p$ corresponding either to $p < 2.66$ (top) or $p > 2.66$ (bottom). Numerically the case $n = -1$ gives worse estimates.</figcaption></figure>

and the *total angular momentum* in the z-direction $J_3$, with $\text{spec}(J_3) = \frac{1}{2}\{1, 2, 3, \dots\}$. See [66, Section 4.6.4] for details. For any $\kappa \in \text{spec}(K)$, we introduce the operator

$$ \mathcal{D}_m^{(\kappa)} := \begin{pmatrix} m - V & \partial_r + \frac{\kappa+1}{r} \\ -\partial_r + \frac{\kappa-1}{r} & -m - V \end{pmatrix} $$

as a self-adjoint operator acting on $L^2(\mathbb{R}^+, r^2 dr)$. If $\lambda_D^{(\kappa)}(V)$ denotes the lowest eigenvalue of $\mathcal{D}_m^{(\kappa)} - V$ in the gap $(-m, m)$, let us define

$$ \Lambda_D^{\text{rad},(\kappa)}(\alpha, p) := \inf \left\{ \lambda_D^{(\kappa)}(V) : V \in L^p(\mathbb{R}^3, \mathbb{R}^+) \text{, } V \text{ radial and } \|V\|_p = \alpha \right\}. $$

We have $\Lambda_D^{\text{rad}}(\alpha, p) = \inf_{\kappa \in \mathbb{Z}\setminus\{0\}} \Lambda_D^{\text{rad},(\kappa)}(\alpha, p)$ and

$$ \Lambda_D(\alpha, p) \le \Lambda_D^{\text{rad}}(\alpha, p) \le \Lambda_D^{\text{rad},(\kappa=1)}(\alpha, p). $$

It is an open question to decide whether the above inequalities are in fact equalities or not. If $\kappa = 1$, we look for an eigenstate of $\mathcal{D}_m - V$ in the Wakano form of [68], that is,

$$ \Psi(r, \theta, \phi) = \begin{pmatrix} \varphi(r) \\ 0 \\ i\chi(r) \cos\theta \\ i e^{i\phi} \chi(r) \sin\theta \end{pmatrix} $$
so that the nonlinear equation becomes

$$ (5.5) \qquad \left\{ \begin{array}{l} \varphi' = - \left( (\lambda + m) + (|\chi|^2 + |\varphi|^2)^{\frac{1}{p-1}} \right) \chi, \\ \chi' + \frac{2}{r} \chi = \left( (\lambda - m) + (|\chi|^2 + |\varphi|^2)^{\frac{1}{p-1}} \right) \varphi. \end{array} \right. $$

System (5.5) provides us with numerical upper estimates of $\Lambda_D(\alpha, p)$: see Fig. 4.

<figure><img src="image_7.png" /><figcaption>Figure 4. Radial case with $d = 3$ and $m = 1$. (Left) The function $p \mapsto \alpha_*^{\text{rad},(\kappa=1)}(p)$ reaches its maximum at $p \approx 3.86$. (Right) The maps $\alpha \mapsto \Lambda_D^{\text{rad},(\kappa=1)}(\alpha, p)$ for values of $p$ corresponding either to $p < 3.86$ (top) or $p > 3.86$ (bottom).</figcaption></figure>

## 5.4. An explicit bound in the radial case in dimensions $d = 2$ or $d = 3$

Let us assume that $m = 1$ and consider at $\lambda = -1$ (lower end of the gap) the system

$$ (5.6) \qquad \varphi' = -W\chi, \quad \chi' + \frac{\delta}{r}\chi = (W-2)\varphi, \quad W^{p-1} = |\varphi|^2 + |\chi|^2. $$

According to the previous section, the radial case $d = 2$ corresponds to $\delta = 1$ (that is $n = 0$ in (5.4)), and the radial case $d = 3$ to $\delta = 2$ (that is $\kappa = 1$ in (5.5)). Writing $\chi(r) = f(r)\varphi(r)$, the equation becomes

$$ \varphi' = -W f \varphi, \quad f' = W(f^2 + 1) - \frac{\delta}{r} f - 2, \quad W^{p-1} = |\varphi|^2 (1 + |f|^2). $$

We now notice that this system admits a solution with $f(r) = r/\mu$ (so that all functions in the middle equality are constant functions). Explicitly, assuming $\delta < p - 1$, with $\mu := \frac{1}{2}(p - 1 - \delta)$, we find a solution of (5.6) given by

$$ (5.7) \qquad \varphi_p(r) = \frac{(p \mu)^{\frac{p-1}{2}} \mu}{(\mu^2 + r^2)^{p/2}}, \quad \chi_p(r) = \frac{(p \mu)^{\frac{p-1}{2}} r}{(\mu^2 + r^2)^{p/2}} \quad \text{and} \quad W_p(r) = \frac{p \mu}{\mu^2 + r^2}. $$
This solution is reminiscent of the solution of [9, Corollary 1.4]. Up to a slight abuse of notations, we can consider $W_p$ as a function of $x \in \mathbb{R}^d$ with $r = |x|$.

**Lemma 5.1.** For all $p \ge d \ge 2$ or $p > 1$ if $d = 1$, and all $\delta < p - 1$, the potential $W_p$ in (5.7), seen as a radial function in $L^p(\mathbb{R}^d)$, satisfies

$$ (5.8) \qquad \|W_p\|_p^p = p^p \pi^{\frac{d}{2}} \left(\frac{2}{p-1-\delta}\right)^{p-d} \frac{\Gamma(p-\frac{d}{2})}{\Gamma(p)} $$

so that in particular $\lim_{p \to d_+} \|W_p\|_p = d \sqrt{\pi} \left(\frac{\Gamma(\frac{d}{2})}{\Gamma(d)}\right)^{1/d}$.

Applied either with $d = 2$ and $\delta = 1$, or $d = 3$ and $\delta = 2$, the expression (5.8) gives an upper bound for $\alpha_*(p)$ in dimension $d = 2$ and $d = 3$. We find that

$$ \alpha_{*}(p)^{p} \leq \frac{p^{p}}{p-1} \left(\frac{2}{p-2}\right)^{p-2} \pi \quad \text{if } d=2, $$

$$ \alpha_{*}(p)^{p} \leq p^{p} \left(\frac{2}{p-3}\right)^{p-3} \frac{\Gamma\left(p-\frac{3}{2}\right)}{\Gamma(p)} \pi^{\frac{3}{2}} \quad \text{if } d=3. $$

In particular,

$$ \alpha_{*}^{\text{rad},(n=0)}(2) \leq 2\sqrt{\pi} \approx 3.54491 \quad \text{if } d=2, $$

$$ \alpha_{*}^{\text{rad},(\kappa=1)}(3) \leq 3 \left(\frac{\pi}{2}\right)^{2/3} \approx 4.05385 \quad \text{if } d=3. $$

The upper bound given by this expression for $d = 1$ (with $\delta = 0$) coincides with the expression found in Theorem 1.3, and we conjecture that we actually have equality in $d = 2$ and $d = 3$ as well. Numerically, the curve $p \mapsto \|W_p\|_p$ coincides with the numerical solution $p \mapsto \alpha_{*}^{\text{rad},(n=0)}(p)$ if $d = 2$ and $p \mapsto \alpha_{*}^{\text{rad},(\kappa=1)}(p)$ if $d = 3$ of Figs. 3 and 4. It is however an open question to decide whether $\varphi_p, \chi_p$ and $W_p$ is the unique solution of (5.6) and if it is optimal among radial optimal functions, and also among non-radial optimal functions (see Appendix B).

# Appendices

## A. Open questions

In this article, we study the ground state defined the lowest eigenvalue in the gap $\lambda_D(V)$ of a general Dirac operator $\not{D}_m - V$ with $V \in L^p(\mathbb{R}^d, \mathbb{R}^+)$ using Birman-Schwinger techniques, and prove that this quantity always makes sense if the $L^p(\mathbb{R}^d)$ norm of $V$ is small enough. To our knowledge, there are several open questions concerning this lowest eigenvalue, which we recall here.

* Is the map $V \mapsto \lambda_D(V)$ concave?

* Is $\lambda_D(V)$ always a simple eigenvalue, or equivalently, is $\mu_1(K_V)$ always simple?

Assuming that the answer of the last question is positive, we denote by $\Psi$ the corresponding eigenfunction for the Dirac operator. We decompose it as $\Psi = \Psi_+ + \Psi_-$ with $\beta \Psi_+ = \Psi_+$ (upper component) and $\beta \Psi_- = -\Psi_-$ (lower component).
• If $V$ is radial (decreasing), is $\Psi_+$ also radial (decreasing)?

Concerning the variational problem associated with (1.5), we recall two questions that were already raised earlier:

• Is the optimal potential $V$ radial (decreasing) if $d \ge 2$?

• If so, is the corresponding ground state $\Psi$ the solution with lowest angular momentum and smallest number of oscillations, as it is suggested in Sections 5.2 and 5.3?

## B. Is the optimal potential radial? A numerical answer

In dimension $d = 2$, we investigate numerically whether the optimal potential $V$ for (1.5) is radial, or equivalently whether the optimal potential $W$ for (3.1) is radial. In order to do so, we run the following self-consistent algorithm<fnref n="1" />. Recall that $K_W := \sqrt{W} R_0(\lambda) \sqrt{W}$ where $R_0$ denotes the resolvent of the free Dirac operator. For $p > d = 2$ and $\lambda \in [-m, m)$, we choose an initial potential $W_0$ at random, and set

$$ \left\{ \begin{array}{ll} \phi_k & := \text{normalized eigenvector corresponding to } \mu_1(K_{W_k}), \\ W_{k+1} & := |\phi_k|^{2/p}. \end{array} \right. $$

In practice, the potential $W_{k+1}$ is also translated so that its maximum is at the origin. We can check that the quantity $\mu_1(K_{W_k})$ is increasing, and that the sequence $(W_k)_{k \in \mathbb{N}}$ converges to some limit potential $W_*$ in $L^p(\mathbb{R}^2)$. A typical run of the algorithm is displayed in Fig 5. In order to check whether $W_*$ is radial or not, we compute the $L^p(\mathbb{R}^2)$ norm of its angular derivative. For $\lambda \in [-0.9, 0.9]$, $m = 1$ and $p \in (2, 8)$, this norm is always much smaller than 1 and usually of the order of $10^{-2}$ or $10^{-3}$, after less than 100 iterations, depending on the parameters we chose. These numerical results suggest that the optimal potentials might be radial, up to translations.

## C. A nonlinear interpolation inequality for the Dirac operator

### C.1. Non-relativistic limit and Keller-Lieb-Thirring inequalities

In order to consider the non-relativistic limit $c \to +\infty$, it is interesting to reintroduce the parameters $\hbar$, $m$ and $c$. The eigenvalue problem

$$ (\mathcal{D}_m^{\hbar,c} - W) \psi = \mu \psi \quad \text{where} \quad \mathcal{D}_m^{\hbar,c} := -i \hbar c \alpha \cdot \nabla + m c^2 \beta $$

is reduced to the eigenvalue problem corresponding to $\hbar = c = m = 1$ by the change of variables

$$ \psi(x) = \Psi\left(\frac{mc}{\hbar}x\right), \quad W(x) = mc^2 V\left(\frac{mc}{\hbar}x\right), \quad \mu = mc^2 \lambda. $$

<footnote n="1">¹The code is available upon request to the authors.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Contour lines of the potential $W_k$ during the iterations, for $p = 3$ and $\lambda = 1/2$, for some $W_0$ chosen at random. The quantities $W_k$ and $\phi_k$ are computed on a square $[-a, a]^2$ with $a = 6$, $L = 100$ discretization points per direction and periodic boundary conditions. The Dirac operator and its inverse are computed in Fourier space and the $L^p(\mathbb{R}^2)$ integrals in direct space.</figcaption></figure>

As a consequence, the ground state $\lambda_D^{\hbar,c}(W)$ of $\mathcal{D}_m^{\hbar,c} - W$ defined as its lowest eigenvalue in the gap $(-mc^2, mc^2)$ and estimated by $\lambda_D(V) \ge \Lambda_D^{(m=1)} (\|V\|_p, p)$ according to the Keller-Lieb-Thirring inequality for the Dirac operator (1.8) becomes

$$ (C.1) \qquad \lambda_D^{\hbar,c}(W) \ge m c^2 \Lambda_D^{(m=1)} \left( \hbar^{-\frac{d}{p}} m^{\frac{d}{p}-1} c^{\frac{d}{p}-2} \|W\|_p, p \right) $$

using the above change of variables. Here $\Lambda_D^{(m=1)}$ stands for $\Lambda_D$ when we assume $m=1$ in notations of Theorem 1.1.

**Proposition C.1.** Let either $d \ge 1$ and $p > 1$ if $d=1$, or $p \ge d$ if $d \ge 2$. With $\Lambda_D^{(m=1)}(\alpha, p)$ defined by (1.5), $\eta = 2p/(2p-d)$ and $K_p$ as in (1.2), we have

$$ 1 - \Lambda_D^{(m=1)}(\alpha, p) = 2^{\frac{d}{2p-d}} K_p \alpha^\eta (1 + o(1)) \quad \text{as} \quad \alpha \to 0_+. $$

If $d=1$, we obtain that

$$ K_p = (p^p (p-1)^{-(p-1)} B(\frac{1}{2}, p))^{-\frac{2}{2p-1}} $$

by expanding the expression of $\alpha_D(\lambda, p)$ given in Theorem 1.3 as $\lambda \to 1_-$ . This is consistent with $K_p = C_q^{-\eta}$ and the expression of the explicit, optimal value of the constant $C_q$ in (1.1) if the dimension is $d=1$: we refer to [23] and references therein for details.
*Proof.* Let us consider the general case $d \ge 1$. The non-relativistic limit of the ground state $\lambda_D^{h,c}(W)$ of the Dirac operator $\mathcal{D}_m^{h,c} - W$ is, up to the mass energy $mc^2$, given by the ground state of the Schrödinger operator

$$-\frac{\hbar^2}{2m} \Delta - W$$

by standard results: see for instance [32, Section 2.4]. Hence

$$\lim_{c \to +\infty} (mc^2 - \lambda_D^{h,c}(W)) = \lambda_S^-(W_\mu) \quad \text{where} \quad W_\mu(x) := W(\mu x) \quad \text{and} \quad \mu = \frac{\hbar}{\sqrt{2m}}.$$

Here $-\lambda_S^-(W_\mu)$ denotes, if it exists, the negative ground state of the Schrödinger operator $-\Delta - W_\mu$. The factor $\mu = \hbar/\sqrt{2m}$ arises from a scaling argument. By definition (1.5), we obtain

$$\lim_{c \to +\infty} mc^2 \left( 1 - \Lambda_D^{(m=1)} \left( \hbar^{-\frac{d}{p}} m^{\frac{d}{p}-1} c^{\frac{d}{p}-2} \|W\|_p, p \right) \right) \le K_p \|W_\mu\|_p^\eta = K_p \mu^{-\frac{d\eta}{p}} \|W\|_p^\eta$$

but there is in fact equality if we use as test function an optimal function $W$ for (1.2). Taking $\alpha = \hbar^{-\frac{d}{p}} m^{\frac{d}{p}-1} c^{\frac{d}{p}-2} \|W\|_p$ in the limit as $c \to +\infty$ concludes the proof. ■

Proposition C.1 is in fact equivalent to

$$(C.2) \qquad \lim_{c \to +\infty} (mc^2 - \lambda_D^{h,c}(W)) \le K_p (\frac{2m}{\hbar^2})^{\frac{d}{2p-d}} \|W\|_p^\eta$$

written with the physical constants. In other words, we recover a standard *Keller-Lieb-Thirring inequality for the Schrödinger operator* (1.2) in the non-relativistic limit. In dimension $d=1$, a tedious but elementary computation directly shows that the constant obtained by taking the non-relativistic limit in the Keller-Lieb-Thirring inequality for the Dirac operator written with optimal constant is the optimal constant in the Keller-Lieb-Thirring inequality for the Schrödinger operator, as it can be deduced for instance from [23, 45].

The definition (1.7) can be generalized to the case $(\hbar, c) \ne (1, 1)$ using the monotonicity of $\alpha \mapsto \Lambda_D^{(m=1)}(\alpha, p)$ stated in Theorem 1.1 and (C.1). If $\alpha_D^{(m=1)}$ denotes the inverse of $\alpha \mapsto \Lambda_D^{(m=1)}(\alpha, p)$, the condition

$$(C.3) \qquad \|W\|_p \le \alpha_D^{h,c}(\lambda, p) := \hbar^{\frac{d}{p}} m^{1-\frac{d}{p}} c^{2-\frac{d}{p}} \alpha_D^{(m=1)}(\frac{\lambda}{mc^2}, p)$$

guarantees that $\lambda_D^{h,c}(W) \ge \lambda$. Notice that $p \ge d$ implies that

$$\lim_{c \to \infty} \|W\|_p \le \alpha_D^{h,c}(\lambda, p) = \infty.$$

## C.2. An interpolation inequality for the Dirac operator

Using a min-max principle as in [28], it is possible to write an optimal interpolation inequality of Gagliardo-Nirenberg-Sobolev type which plays for the free Dirac operator
the same role as (1.1). The inequality is somewhat involved, but Inequality (1.1) is recovered in the non-relativistic limit as $c \to +\infty$. For sake of simplicity, we consider only the case $d=1$.

Let us start by a short and formal summary of the min-max principle applied to the determination of the ground state of the Dirac operator. If $(\varphi, \chi)^{\top}$ is an eigenspinor of the operator $\mathcal{D}_m^{\frac{\hbar,c}{2}} - V$ with eigenvalue $\lambda \in (-mc^2, mc^2)$, then, as in (5.1) we have

$$ \begin{cases} \hbar c \varphi' = -(\lambda + mc^2 + V) \chi, \\ \hbar c \chi' = (\lambda - mc^2 + V) \varphi. \end{cases} $$

The first line gives

$$ \chi = -\hbar c \frac{\varphi'}{\lambda + mc^2 + V} $$

so that the problem amounts to solving

$$ -(\hbar c)^2 \left( \frac{\varphi'}{\lambda + mc^2 + V} \right)' + (mc^2 - \lambda - V) \varphi = 0. $$

Multiplying by $\varphi$ and integrating suggests to introduce the functional

$$ \mathcal{E}[\mu, V, \phi] := (\hbar c)^2 \int_{\mathbb{R}} \frac{|\phi'|^2}{\mu + mc^2 + V} dx + \int_{\mathbb{R}} (mc^2 - \mu - V) |\phi|^2 dx. $$

Clearly, we have $\mathcal{E}[\lambda, V, \phi] = 0$. In addition, for all fixed $V$ and $\phi$, the map $\mu \mapsto \mathcal{E}[\mu, V, \phi]$ is decreasing. It is proved in [63, Lemma 2.4], that for all $-m < \mu < \lambda_D^{\frac{\hbar,c}{2}}(V)$, the quadratic map $\phi \mapsto \mathcal{E}[\mu, V, \phi]$ is positive definite, and that, for $\mu = \lambda_D(V)$, we have $\mathcal{E}[\mu, V, \phi] = 0$ if and only if $\phi = \varphi$, up to a multiplicative constant. In particular, we have

$$ \forall \phi \in C_0^\infty(\mathbb{R}), \quad \forall V \in L^p(\mathbb{R}), \quad \|V\|_p \le \alpha_D^{\frac{\hbar,c}{2}}(\lambda, p) \implies \mathcal{E}[\lambda, V, \phi] \ge 0. $$

Minimizing $\mathcal{E}[\lambda, W, \phi]$ in $W$ such that $\|W\|_p = \alpha \le \alpha_D^{\frac{\hbar,c}{2}}(\lambda, p)$ shows that the optimal $W$ solves the Euler-Lagrange equation of the implicit form

$$ (C.4) \qquad \nu W^{p-1} = |\phi|^2 + \frac{(\hbar c)^2 |\phi'|^2}{(\lambda + mc^2 + W)^2}, $$

where $\nu \ge 0$ is now the Lagrange multiplier for the constraint $\|W\|_p = \alpha$. Note that for all fixed $a, b, c \ge 0$, the equation

$$ \nu X^{p-1} = a + \frac{b}{(c+X)^2} $$

has a unique solution in $X_\nu \ge 0$, as the left-hand side is an increasing function of $X$, while the right-hand side is decreasing, and that $\nu \mapsto X_\nu$ is increasing. So for fixed $\nu \ge 0$, there is a unique $W = V_\nu[\phi]$ satisfying (C.4) and the map $\nu \mapsto V_\nu$ is pointwise decreasing, hence so is the map $\nu \mapsto \|V_\nu\|_p$. With $\alpha_D^{\frac{\hbar,c}{2}}(\lambda, p)$ given by (C.3), we define

$$ \nu_*(\lambda, p, \phi) := \inf \left\{ \nu > 0 : \|V_\nu[\phi]\|_p \le \alpha_D^{\frac{\hbar,c}{2}}(\lambda, p) \right\}. $$
Summarizing, we proved that for all $\phi \in C_0^\infty(\mathbb{R})$ and all $\nu \ge \nu_*(\lambda, p, \phi)$,

$$ (C.5) \qquad (\hbar c)^2 \int_{\mathbb{R}} \frac{|\phi'|^2}{\lambda + m c^2 + V_\nu[\phi]} dx + \int_{\mathbb{R}} (m c^2 - \lambda - V_\nu[\phi]) |\phi|^2 dx \ge 0, $$

which can be interpreted as a *Gagliardo–Nirenberg type inequality* for $\phi$ alone. Such an inequality is known for a fixed, given potential $V$ from [21, 27, 28] and it is then of Hardy-type, as for instance the new Hardy inequality in [35], but the novelty in this paper is that we take $V = V_\nu[\phi]$ thus making it a non-linear interpolation inequality. While the form (C.5) is non-explicit, it allows to recover the usual Gagliardo–Nirenberg inequality in the non-relativistic limit as $c \to \infty$. By writing $\lambda = m c^2 + E$ for some $E < 0$, (C.5) becomes

$$ (\hbar c)^2 \int_{\mathbb{R}} \frac{|\phi'|^2}{2m c^2 + E + V_\nu[\phi]} dx - \int_{\mathbb{R}} (E + V_\nu[\phi]) |\phi|^2 dx \ge 0. $$

Let us choose $\nu = \|\phi\|_{2p/(p-1)}^2$. As $c \to \infty$, we get from (C.4) that $V_\nu[\phi]$ converges to $\|\phi\|_{2p/(p-1)}^{-2/(p-1)} |\phi|^{2/(p-1)}$. Together with (C.2), we get that $\nu \ge \nu_*(\lambda, p, \phi)$ in the limit $c \to \infty$ whenever $|E| \ge K_p (2m/\hbar^2)^{d/(2p-d)}$. We obtain

$$ \frac{\hbar^2}{2m} \int_{\mathbb{R}} |\phi'|^2 dx - \|\phi\|_{\frac{2p}{p-1}}^2 \ge E \int_{\mathbb{R}} |\phi|^2 dx. $$

This inequality is the Gagliardo–Nirenberg inequality (1.1) written in non-scale invariant form, for an appropriate choice of the parameter $\lambda$ in (1.1).

## D. The case $p = d = 1$

This appendix deals with the limit case $p=1$ of Theorem 1.3 devoted to the one-dimension- al Keller estimates. We give a computation of $\alpha_D(\lambda, 1)$ which is not based on the limit as $p \to 1_+$ of the nonlinear estimates and prove that any sequence of optimizing potentials concentrates into a Dirac $\delta$ distribution.

**Proposition D.1.** If $d=1$, then $\alpha_D(\lambda, 1) = \arccos(\lambda/m)$. More specifically, for all $\alpha \in (0, \pi)$, all $V \in L^1(\mathbb{R}, \mathbb{R}^+)$ with $\|V\|_1 = \alpha$, if $\lambda \in (-m, m)$ is an eigenvalue of $\mathcal{D}_m - V$, then we have the strict inequality

$$ m \cos \alpha < \lambda. $$

In addition, any sequence of nonnegative potentials $(V_n)_{n \in \mathbb{N}}$ with $\|V_n\|_1 = \alpha$ and eigen- values $\lambda_n$ approaching $m \cos \alpha$, converges as $n \to +\infty$ to a Dirac $\delta$ distribution.

According to [67], “the method of directly solving the Dirac equation with a $\delta$-function potential and the method of obtaining the solution by first solving the Dirac equation with a short-range potential and afterward taking the $\delta$-function limit, lead to different results” [concerning the spectrum]. This issue is known as Klein’s paradox. Although the Keller-Lieb-Thirring (1.12) makes sense for any nonnegative potential $V \in L^1(\mathbb{R}^d)$, it is a natural question to investigate by direct methods whether the bound is achieved in the larger set of bounded nonnegative measures and consider sequences of optimizing potentials.
*Proof.* We start with a calculation for a bounded and compactly supported potential $V$. In this case, the eigenvalue equation rewrites as

$$ \Psi' = (i\sigma_2(V + \lambda) + m\sigma_1)\Psi. $$

We decompose $\Psi$ on the (not-orthonormal) basis given by the eigenvectors $e_{\pm}$ of the matrix $i\lambda\sigma_2 - m\sigma_1$ defined by

$$ e_{\pm} := \begin{pmatrix} \sqrt{m^2 - \lambda^2} \\ \pm(m - \lambda) \end{pmatrix} \quad \text{such that} \quad \begin{pmatrix} 0 & m + \lambda \\ m - \lambda & 0 \end{pmatrix} e_{\pm} = \pm \sqrt{m^2 - \lambda^2} e_{\pm}. $$

Decomposing $\Psi(x) = a(x)e_+ + b(x)e_-$ and using the identities

$$ \begin{align*} \langle i\sigma_2 e_{\pm}, e_{\pm} \rangle &= 0, & \langle i\sigma_2 e_{\pm}, e_{\mp} \rangle &= \pm 2(m-\lambda)\sqrt{m^2-\lambda^2}, \\ \langle i\sigma_2 e_{\pm}, i\sigma_2 e_{\pm} \rangle &= 2m(m-\lambda), & \langle i\sigma_2 e_{\pm}, i\sigma_2 e_{\mp} \rangle &= 2\lambda(m-\lambda), \end{align*} $$

gives, with $W := V/\sqrt{m^2 - \lambda^2}$,

$$ \begin{align*} a' &= (\sqrt{m^2 - \lambda^2} - \lambda W)a - mWb, \\ b' &= -(\sqrt{m^2 - \lambda^2} - \lambda W)b + mWa. \end{align*} $$

Since $V$ (and $W$) are compactly supported, a square-integrable solution must have $b = 0$ in a neighborhood of $-\infty$ and $a = 0$ in a neighborhood of $+\infty$. Without loss of generality, we take a solution with $a(x) > 0$ for $x$ near $-\infty$. Since

$$ (ab)' = mW(a^2 - b^2) $$

is nonnegative if $|a| > |b|$, such a solution enters the first (a, b) quadrant and stays in the first quadrant until the first value of $x$ such that $a(x) = 0$. We denote this value by $x_1$, with $x_1 = +\infty$ if a does not change sign. In the interval $(-\infty, x_1)$, the ratio $t := b/a$ is well-defined and satisfies

$$ \begin{align*} t' &= \frac{1}{a^2} \left( -2ab\left(\sqrt{m^2 - \lambda^2} - \lambda W\right) + mW(a^2 + b^2) \right) \\ &= -2\sqrt{m^2 - \lambda^2}t + W(m + mt^2 + 2\lambda t). \end{align*} $$

We finally define the angle

$$ \theta_{\lambda}(t) := \arctan \left( \frac{mt + \lambda}{\sqrt{m^2 - \lambda^2}} \right), $$

such that $\lim_{x \to x_1} \theta_{\lambda}(t(x)) = \pi/2$ and

$$ (D.1) \qquad \frac{1}{\sqrt{m^2 - \lambda^2}} (\theta_{\lambda} \circ t)' = \frac{t'}{m + 2\lambda t + m t^2} = W - 2 \frac{\sqrt{m^2 - \lambda^2} t}{m + 2\lambda t + m t^2}. $$

Integrating for $x \in (-\infty, x_1)$, we obtain

$$ (D.2) \qquad \frac{\pi/2 - \theta_{\lambda}(0)}{\sqrt{m^2 - \lambda^2}} = \int_{-\infty}^{x_1} W(s) \, ds - 2 \int_{-\infty}^{x_1} \frac{\sqrt{m^2 - \lambda^2} t}{m + 2\lambda t + m t^2} \, ds < \frac{\alpha}{\sqrt{m^2 - \lambda^2}}. $$
Since $\theta_{\lambda}(0) = \arcsin(\lambda/m)$, we obtain

$$ \arccos(\lambda/m) < \alpha \quad \text{or} \quad \lambda > m \cos \alpha. $$

In order to approximate unbounded potentials, we need an estimate on the negative term in (D.2). Take any number $c > 1$. Since $t$ is continuous, there is an interval $I_V(c) \subset (-\infty, x_1]$ such that $t(x) \in (1/c, c)$ for all $x \in I_V(c)$. We have the bound (note that the integrand is symmetric under $t \mapsto 1/t$)

$$ \int_{I_V(c)} \frac{2(m^2 - \lambda^2) t(s)}{m + m t(s)^2 + 2 \lambda t(s)} ds \geq \frac{2(m^2 - \lambda^2) c}{m + m c^2 + 2 \lambda c} |I_V(c)| $$

and therefore

$$ (D.3) \qquad \arccos(\lambda/m) \le \alpha - \frac{2(m^2 - \lambda^2)c}{m + m c^2 + 2\lambda c} |I_V(c)|. $$

To prove that $|I_V(c)|$ cannot be arbitrarily small, we integrate (D.1) on $I_V(c)$, which gives

$$ (D.4) \qquad \theta_{\lambda}(c) - \theta_{\lambda}(1/c) \le \int_{I_V(c)} V(s) ds \le Q_V(|I_V(c)|) $$

where we have defined

$$ Q_V(r) := \sup_{x \in \mathbb{R}} \int_{x-r/2}^{x+r/2} V(s) ds. $$

Now, assume that $(V_n)_{n \in \mathbb{N}}$ is a sequence of potentials with $\|V_n\|_1 = \alpha$ and eigenvalues $\lambda_n$ converging to $\lambda := m \cos \alpha$. Without loss of generality, we may assume that each $V_n$ is bounded and compactly supported. By (D.3), in order to approach the equality case, we need that $|I_{V_n}(c)|$ tends to zero for each $c > 1$. We now use (D.4) to show that this implies the convergence (after suitable translations) to a Dirac $\delta$ distribution.

Fix $\epsilon > 0$ and $r > 0$. Fix $c > 1$ and $n_0$ such that for all $n \ge n_0$, we have

$$ \theta_{\lambda_n}(c) - \theta_{\lambda_n}(1/c) \ge \theta_{\lambda}(+\infty) - \theta_{\lambda}(0) - \epsilon = \alpha - \epsilon. $$

Upon increasing $n_0$, we can assume $|I_{V_n}(c)| \le r$ for all $n \ge n_0$. From (D.4), this gives

$$ Q_{V_n}(r) \ge Q_{V_n}(|I_{V_n}(c)|) \ge \alpha - \epsilon. $$

Since $r$ and $\epsilon$ are arbitrary, we have shown that $Q_{V_n}$ converges pointwise to $\alpha$. In the language of concentration-compactness, this excludes *vanishing* and *dichotomy* and implies that, after a sequence of translations, $V_n$ converges to a measure of total mass $\alpha$ supported at the origin, hence, to a Dirac $\delta$ distribution. ■

**Acknowledgements.** The authors thank the referees for their suggestions which have all been taken into account and led to a significant improvement of the paper.
**Funding.** This work was partially supported by the Project EFI (ANR-17-CE40-0030) of the French National Research Agency (ANR). HVDB received funding from the Center for Mathematical Modeling (Universidad de Chile & CNRS IRL 2807) through ANID/Basal projects #FB210005 and #ACE210010, as well as ANID/Fondecyt project #11220194, and MathAmSud project EEquaddii 20-MATH-04. This work was partially developed when FP was employed at CNRS & CEREMADE-Université Paris-Dauphine and supported by the project ANR-17-CE29-0004 molQED of the ANR and by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement MDFT No. 725528). He is also supported by the project PID2021-123034NB-100 funded by MCIN/AEI/10.13039/501100011033/ FEDER, UE.

© 2023 by the authors. Reproduction of this article by any means permitted for noncommercial purposes. CC-BY 4.0

References

[1] M. ABRAMOWITZ AND I. A. STEGUN, *Handbook of mathematical functions with formulas, graphs, and mathematical tables. 10th printing, with corrections.* National Bureau of Standards. A Wiley-Interscience Publication. New York etc.: John Wiley & Sons. xiv, 1046 pp. £42.70 (1972)., 1972.

[2] S. ALAMA, P. A. DEIFT, AND R. HEMPEL, *Eigenvalue branches of the Schrödinger operator $H - λW$ in a gap of $σ(H)$*, Comm. Math. Phys., 121 (1989), pp. 291–321.

[3] N. ARRIZABALAGA, J. DUOANDIKOETXEA, AND L. VEGA, *Self-adjoint extensions of Dirac operators with Coulomb type singularity*, J. Math. Phys., 54 (2013), pp. 041504, 20 [Erratum: J. Math. Phys., 59 (7): 079902, 2 (2018)].

[4] M. BALABANE, T. CAZENAVE, A. DOUADY, AND F. MERLE, *Existence of excited states for a nonlinear Dirac field*, Comm. Math. Phys., 119 (1988), pp. 153–176.

[5] W. BORRELLI, *Stationary solutions for the 2D critical Dirac equation with Kerr nonlinearity*, J. Differential Equations, 263 (2017), pp. 7941–7964.

[6] ——, *Multiple solutions for a self-consistent Dirac equation in two dimensions*, J. Math. Phys., 59 (2018), pp. 041503, 13.

[7] ——, *Weakly localized states for nonlinear Dirac equations*, Calc. Var. Partial Differential Equations, 57 (2018), pp. Paper No. 155, 21.

[8] W. BORRELLI AND R. L. FRANK, *Sharp decay estimates for critical Dirac equations*, Trans. Amer. Math. Soc., 373 (2020), pp. 2045–2070.

[9] W. BORRELL, A. MALCHIODI, AND R. WU, *Ground state Dirac bubbles and Killing spinors*, Comm. Math. Phys., 383 (2021), pp. 1151–1180.

[10] R. BOSI, J. DOLBEAULT, AND M. J. ESTEBAN, *Estimates for the optimal constants in multipolar Hardy inequalities for Schrödinger and Dirac operators*, Commun. Pure Appl. Anal., 7 (2008), pp. 533–562.

[11] H. BREZIS, *Functional Analysis, Sobolev Spaces and Partial Differential Equations*, Springer New York, 2011.

[12] E. A. CARLEN, R. L. FRANK, AND E. H. LIEB, *Stability estimates for the lowest eigenvalue of a Schrödinger operator*, Geom. Funct. Anal., 24 (2014), pp. 63–84.

[13] B. CASSANO AND F. PIZZICHELLO, *Self-adjoint extensions for the Dirac operator with Coulomb-type spherically symmetric potentials*, Lett. Math. Phys., 108 (2018), pp. 2635–2667.
[14] B. CASSANO, F. PIZZICHILLO, AND L. VEGA, *A Hardy-type inequality and some spectral characterizations for the Dirac-Coulomb operator*, Rev. Mat. Complut., 33 (2020), pp. 1-18.

[15] T. CAZENAVE AND L. VÁZQUEZ, *Existence of localized solutions for a classical nonlinear Dirac field*, Comm. Math. Phys., 105 (1986), pp. 35-47.

[16] C. V. COFFMAN, *Uniqueness of the ground state solution for $\Delta u - u + u^3 = 0$ and a variational characterization of other solutions*, Arch. Rational Mech. Anal., 46 (1972), pp. 81-95.

[17] J.-C. CUENIN, *Eigenvalue bounds for Dirac and fractional Schrödinger operators with complex potentials*, J. Funct. Anal., 272 (2017), pp. 2987-3018.

[18] J.-C. CUENIN, A. LAPTEV, AND C. TRETTER, *Eigenvalue estimates for non-selfadjoint Dirac operators on the real line*, Ann. Henri Poincaré, 15 (2014), pp. 707-736.

[19] I. DAUBECHIES, *An uncertainty principle for fermions with generalized kinetic energy*, Communications in Mathematical Physics, 90 (1983), pp. 511-520.

[20] J. DOLBEAULT, M. ESTEBAN, AND M. LOSS, *Characterization of the critical magnetic field in the Dirac-Coulomb equation*, Journal of Physics A: Mathematical and Theoretical, 41 (2008), p. 185303 (13pp).

[21] J. DOLBEAULT, M. J. ESTEBAN, J. DUOANDIKOETXA, AND L. VEGA, *Hardy-type estimates for Dirac operators*, Ann. Sci. Éc. Norm. Supér., 40 (2007), pp. 885-900.

[22] J. DOLBEAULT, M. J. ESTEBAN, AND A. LAPTEV, *Spectral estimates on the sphere*, Analysis & PDE, 7 (2014), pp. 435-460.

[23] J. DOLBEAULT, M. J. ESTEBAN, A. LAPTEV, AND M. LOSS, *One-dimensional Gagliardo-Nirenberg-Sobolev inequalities: remarks on duality and flows*, J. London Math. Soc., 90 (2014), pp. 525-550.

[24] ——, *Interpolation inequalities and spectral estimates for magnetic operators*, Ann. Henri Poincaré, 19 (2018), pp. 1439-1463.

[25] J. DOLBEAULT, M. J. ESTEBAN, AND M. LOSS, *Relativistic hydrogenic atoms in strong magnetic fields*, Ann. Henri Poincaré, 8 (2007), pp. 749-779.

[26] J. DOLBEAULT, M. J. ESTEBAN, AND M. LOSS, *Critical magnetic field for 2d magnetic Dirac-Coulomb operators and Hardy inequalities*, in *Partial Differential Equations, Spectral Theory, and Mathematical Physics, The Ari Laptev Anniversary Volume*, EMS Series of Congress Reports, EMS Press, jun 2021, pp. 41-63.

[27] J. DOLBEAULT, M. J. ESTEBAN, M. LOSS, AND L. VEGA, *An analytical proof of Hardy-like inequalities related to the Dirac operator*, J. Funct. Anal., 216 (2004), pp. 1-21.

[28] J. DOLBEAULT, M. J. ESTEBAN, AND E. SÉRÉ, *On the eigenvalues of operators with gaps. Application to Dirac operators*, J. Funct. Anal., 174 (2000), pp. 208-226.

[29] ——, *General results on the eigenvalues of operators with gaps, arising from both ends of the gaps. Application to Dirac operators*, J. Eur. Math. Soc. (JEMS), 8 (2006), pp. 243-251.

[30] ——, *Corrigendum to: "On the eigenvalues of operators with gaps. Application to Dirac operators"* [J. Funct. Anal. 174 (1) (2000) 208-226], J. Funct. Anal., 284 (2023), p. 109651.

[31] J. DOLBEAULT, M. J. ESTEBAN, AND E. SÉRÉ, *Distinguished self-adjoint extension and eigenvalues of operators with gaps. Application to Dirac-Coulomb operators*, arXiv: 2206.11679, to appear in J. Spectral Theory, (2023).

[32] M. J. ESTEBAN, M. LEWIN, AND E. SÉRÉ, *Variational methods in relativistic quantum mechanics*, Bull. Amer. Math. Soc., 45 (2008), pp. 535-535.
[33] ——, *Domains for Dirac-Coulomb min-max levels*, Rev. Mat. Iberoam., 35 (2019), pp. 877–924.

[34] ——, *Dirac-Coulomb operators with general charge distribution I. Distinguished extension and min-max formulas*, Annales Henri Lebesgue, 4 (2021), pp. 1421–1456.

[35] ——, *Dirac-Coulomb operators with general charge distribution II. The lowest eigenvalue*, Proc. Lond. Math. Soc. (3), 123 (2021), pp. 345–383.

[36] M. J. ESTEBAN AND M. LOSS, *Self-adjointness for Dirac operators via Hardy-Dirac inequalities*, J. Math. Phys., 48 (2007), pp. 112107, 8.

[37] ——, *Self-adjointness via partial Hardy-like inequalities*, in Mathematical results in quantum mechanics, World Sci. Publ., Hackensack, NJ, 2008, pp. 41–47.

[38] M. J. ESTEBAN AND E. SÉRÉ, *Stationary states of the nonlinear Dirac equation: A variational approach*, Comm. Math. Phys., 171 (1995), pp. 323–350.

[39] L. FANELLI AND D. KREJČÍŘÍK, *Location of eigenvalues of three-dimensional non-self-adjoint Dirac operators*, Lett. Math. Phys., 109 (2019), pp. 1473–1485.

[40] R. L. FRANK AND B. SIMON, *Critical Lieb-Thirring bounds in gaps and the generalized Nevai conjecture for finite gap Jacobi matrices*, Duke Math. J., 157 (2011), pp. 461–493.

[41] T. FRIEDRICH, *Dirac Operators in Riemannian Geometry*, American Mathematical Society, aug 2000.

[42] I. S. GRADSHTEYN AND I. M. RYZHIK, *Table of integrals, series, and products*, Elsevier / Academic Press, Amsterdam, eighth ed., 2015. Translated from the Russian, Translation edited and with a preface by Daniel Zwillinger and Victor Moll, Revised from the seventh edition.

[43] T. KATO, *Fundamental properties of Hamiltonian operators of Schrödinger type*, Trans. Amer. Math. Soc., 70 (1951), p. 195.

[44] ——, *Holomorphic families of Dirac operators*, Mathematische Zeitschrift, 183 (1983), pp. 399–406.

[45] J. B. KELLER, *Lower bounds and isoperimetric inequalities for eigenvalues of the Schrödinger equation*, J. Mathematical Phys., 2 (1961), pp. 262–266.

[46] M. KLAUS, *Dirac operators with several Coulomb singularities*, Helv. Phys. Acta, 53 (1980), pp. 463–482 (1981).

[47] M. KLAUS AND R. WÜST, *Spectral properties of Dirac operators with singular potentials*, J. Math. Anal. Appl., 72 (1979), pp. 206–214.

[48] R. KONNO AND S. T. KURODA, *On the finiteness of perturbed eigenvalues*, J. Fac. Sci. Univ. Tokyo Sect. I, 13 (1966), pp. 55–63 (1966).

[49] M. K. KWONG, *Uniqueness of positive solutions of $\Delta u - u + u^p = 0$ in $\mathbb{R}^n$*, Arch. Rational Mech. Anal., 105 (1989), pp. 243–266.

[50] E. LIEB, M. LOSS, AND H. SIEDENTOP, *Stability of relativistic matter via Thomas-Fermi theory*, Helvetica Physica Acta, 69 (1996), pp. 974–984.

[51] E. LIEB AND W. THIRRING, *Bound on kinetic energy of fermions which proves stability of matter*, Phys. Rev. Lett., 35 (1975), pp. 687–689.

[52] E. H. LIEB AND R. SEIRINGER, *The stability of matter in quantum mechanics*, Cambridge University Press, 2010.
[53] E. H. LIEB AND W. E. THIRRING, *Inequalities for the moments of the eigenvalues of the Schrödinger Hamiltonian and their relation to Sobolev inequalities*, Essays in Honor of Valentine Bargmann, E. Lieb, B. Simon, A. Wightman Eds. Princeton University Press, 1976, pp. 269–303.

[54] P.-L. LIONS, *The concentration-compactness principle in the calculus of variations. The locally compact case*, Part I, Ann. Inst. H. Poincaré Anal. Non Linéaire, 1 (1984), pp. 109–149.

[55] K. McLEOD, *Uniqueness of positive radial solutions of $\Delta u + f(u) = 0$ in $\mathbb{R}^n$. II*, Trans. Amer. Math. Soc., 339 (1993), pp. 495–505.

[56] F. MERLE, *Existence of stationary states for nonlinear Dirac equations*, J. Differential Equations, 74 (1988), pp. 50–68.

[57] G. NENCIU, *Self-adjointness and invariance of the essential spectrum for Dirac operators defined as quadratic forms*, Comm. Math. Phys., 48 (1976), pp. 235–247.

[58] G. NENCIU, *Distinguished self-adjoint extension for Dirac operator with potential dominated by multicenter Coulomb potentials*, Helv. Phys. Acta, 50 (1977), pp. 1–3.

[59] D. PELINOVSKY, *Survey on global existence in the nonlinear Dirac equations in one spatial dimension*, in Harmonic analysis and nonlinear partial differential equations, RIMS Kôkyûroku Bessatsu, B26, Res. Inst. Math. Sci. (RIMS), Kyoto, 2011, pp. 37–50.

[60] M. REED AND B. SIMON, *Methods of modern mathematical physics. II. Fourier analysis, self-adjointness*, Academic Press [Harcourt Brace Jovanovich, Publishers], New York-London, 1975.

[61] ——, *Methods of modern mathematical physics. IV. Analysis of operators*, Academic Press [Harcourt Brace Jovanovich, Publishers], New York-London, 1978.

[62] F. RELLICH, *Störungstheorie der Spektralzerlegung*, Math. Ann., 113 (1937), pp. 600–619.

[63] L. SCHIMMER, J. P. SOLOVEJ, AND S. TOKUS, *Friedrichs extension and min-max principle for operators with a gap*, Ann. Henri Poincaré, 21 (2020), pp. 327–357.

[64] U.-W. SCHMINCKE, *Distinguished selfadjoint extensions of Dirac operators*, Math. Z., 129 (1972), pp. 335–349.

[65] B. SIMON, *Trace ideals and their applications*, vol. 120 of Mathematical Surveys and Monographs, American Mathematical Society, Providence, RI, second ed., 2005.

[66] B. THALLER, *The Dirac equation*, Texts and Monographs in Physics, Springer-Verlag, Berlin, 1992.

[67] P. ŠEBA, *Klein's paradox and the relativistic point interaction*, Lett. Math. Phys., 18 (1989), pp. 77–86.

[68] M. WAKANO, *Intensely localized solutions of the classical Dirac-Maxwell field equations*, Progress of Theoretical Physics, 35 (1966), pp. 1117–1141.

[69] R. WÜST, *A convergence theorem for selfadjoint operators applicable to Dirac operators with cutoff potentials*, Math. Z., 131 (1973), pp. 339–349.

[70] ——, *Distinguished self-adjoint extensions of Dirac operators constructed by means of cutoff potentials*, Math. Z., 141 (1975), pp. 93–98.

[71] ——, *Dirac operations with strongly singular potentials. Distinguished self-adjoint extensions constructed with a spectral gap theorem and cut-off potentials*, Math. Z., 152 (1977), pp. 259–271.
**Jean Dolbeault**

CEREMADE (CNRS UMR n° 7534), PSL University, Université Paris-Dauphine, Place de Lattre de Tassigny, 75775 Paris 16, France;
dolbeaul@ceremade.dauphine.fr

**David Gontier**

CEREMADE (CNRS UMR n° 7534), PSL University, Université Paris-Dauphine, Place de Lattre de Tassigny, 75775 Paris 16, France & ENS/PSL University, DMA, F-75005, Paris, France;
gontier@ceremade.dauphine.fr

**Fabio Pizzichillo**

Departamento de Matemática e Informática aplicadas a las Ingenierías Civil y Naval, Universidad Politécnica de Madrid, Escuela Técnica Superior de Ingenieros de Caminos, Canales y Puertos, Calle Profesor Aranguren, 3. 28040, Madrid, Spain;
fabio.pizzichillo@upm.es

**Hanne Van Den Bosch**

Departamento de Ingeniería Matemática & Center for Mathematical Modeling, Facultad de Ciencias Físicas y Matemáticas, Universidad de Chile and CNRS IRL 2807, Beauchef 851, Piso 5, Santiago, Chile;
hvdbosch@dim.uchile.cl