# Some Remarks on the Regularized Hamiltonian for Three Bosons with Contact
Interactions

Daniele Ferretti  Department of Mathematics G. Castelnuovo, University of Rome
“La Sapienza”
Piazzale Aldo Moro, 5, 00185 Rome, Italy<EMAIL_ADDRESS>and
Alessandro Teta  Department of Mathematics G. Castelnuovo, University of Rome
“La Sapienza”,
Piazzale Aldo Moro, 5, 00185 Rome, Italy<EMAIL_ADDRESS>
###### Abstract.

We discuss some properties of a model Hamiltonian for a system of three bosons
interacting via zero-range forces in three dimensions. In order to avoid the
well known instability phenomenon, we consider a regularized Hamiltonian with
a repulsive three-body force. We review the main result recently obtained in
[1] where, starting from a suitable quadratic form $Q$, the self-adjoint and
bounded from below Hamiltonian $\mathcal{H}$ is constructed provided that the
strength $\gamma$ of the three-body repulsion is larger than a threshold
parameter $\gamma_{c}\mspace{1.5mu}$. We also show that the threshold value
$\gamma_{c}$ found in [1] is optimal, in the sense that the quadratic form $Q$
is unbounded from below if $\gamma<\gamma_{c}\mspace{1.5mu}$. Finally, we give
an alternative and much simpler proof of the result in [1] whenever
$\gamma>\gamma^{\prime}_{c}\mspace{1.5mu}$, with $\gamma^{\prime}_{c}$
strictly larger than $\gamma_{c}\mspace{1.5mu}$.

The authors acknowledge the support of the GNFM Gruppo Nazionale per la Fisica
Matematica - INdAM

_Keywords: Zero-range interactions; Three-body Hamiltonians; Schrödinger
operators._

_MSC 2020: 81Q10; 81Q15; 70F07; 46N50;_

##  1\. Introduction

In this note we discuss some properties of a model Hamiltonian describing the
dynamics of three identical bosons interacting via zero-range forces in
dimension three. Since the seminal papers by Minlos and Faddeev ([12, 13]), it
is known that a natural candidate for such Hamiltonian turns out to be
unbounded from below, giving rise to the so-called Thomas effect. Here natural
means that the boundary condition defining the Hamiltonian (known as TMS
boundary condition) is the direct generalization to the three-body case of the
boundary condition characterizing the Hamiltonian of the two-body problem.
Roughly speaking, the reason of such instability is due to the interaction
becoming too singular as all the three particles are close to each other. We
note that this pathology is absent in dimension one, where perturbation theory
of quadratic forms can be used, and in dimension two, where the renormalized
two-body boundary condition is sufficient to avoid the collapse (see e.g. [3,
5]).

It is worth to underline that the construction of a self-adjoint and bounded
from below Hamiltonian for three, or more, interacting bosons with zero-range
forces in dimension three is a challenging open problem in Mathematical
Physics. Following a suggestion contained in [12], it has been recently
studied ([7], [11], [1]) a regularized version of the Hamiltonian for a system
of three bosons (see also [6] for the case of $N$ bosons interacting with an
impurity). The main idea is to introduce a three-body repulsion that reduces
to zero the strength of the contact interaction between two particles if the
third particle approaches the common position of the first two. On the other
hand, when the third particle is far enough, the usual two-body point
interaction is restored. The result is that the regularized Hamiltonian is
self-adjoint and bounded from below if the strength $\gamma$ of the three-body
interaction is larger than a threshold value $\gamma_{c}\mspace{1.5mu}$.

The aim of this paper is to describe the construction of such regularized
Hamiltonian following the approach developed in [1] and also to prove two
further results. More precisely, in section 2 we introduce the notation and we
formulate the main result of [1], essentially based on the analysis of a
suitable quadratic form $Q$.

In section 3 we prove that the threshold value $\gamma_{c}$ obtained in [1] is
optimal, in the sense that for $\gamma<\gamma_{c}$ the quadratic form $Q$ is
unbounded from below.

In section 4 we give a different proof of the main result in [1] based on a
new approach in position space. The proof is surely less general since it is
valid only for $\gamma>\gamma^{\prime}_{c}\mspace{1.5mu}$, where
$\gamma^{\prime}_{c}>\gamma_{c}\mspace{1.5mu}$. On the other hand it has the
advantage to be much simpler and to show that the choice of the three-body
force is not arbitrary but it is dictated by the inherent singularity of the
problem.

##  2\. Regularized Hamiltonian

Let us consider a system composed of three identical spinless bosons of mass
$\frac{1}{2}$ in three dimensions and let us fix the center of mass reference
frame so that $\bm{x}_{1}$, $\bm{x}_{2}$ and
$\bm{x}_{3}=-\bm{x}_{1}-\bm{x}_{2}$ represent the Cartesian coordinates of the
three particles. We also introduce the Jacobi coordinates

$\begin{cases}\bm{r}_{k}=\frac{1}{2}\sum_{i,j=1}^{3}\epsilon_{ijk}(\bm{x}_{i}-\bm{x}_{j}),\\\
\bm{\rho}_{k}=\frac{3}{2}\,\bm{x}_{k}-\frac{1}{2}\sum_{\ell=1}^{3}\bm{x}_{\ell},\end{cases}\qquad
k\in\\{1,2,3\\}$ (2.1)

where $\epsilon_{ijk}$ is the Levi-Civita symbol, so that one has the
following identities

$\begin{cases}\bm{r}_{k\pm 1}=-\frac{1}{2}\bm{r}_{k}\mp\bm{\rho}_{k},\\\
\bm{\rho}_{k\pm
1}=\pm\frac{3}{4}\bm{r}_{k}-\frac{1}{2}\bm{\rho}_{k},\end{cases}\qquad
k\in\mathbb{Z}\mspace{1.5mu}\diagup_{\mspace{-4.5mu}\\{3\\}}.$ (2.2)

Denoting by $\bm{x}=\bm{r}_{1}=\bm{x}_{2}-\bm{x}_{3}$ and
$\bm{y}=\bm{\rho}_{1}=\bm{x}_{1}-\frac{\bm{x}_{2}\,+\,\bm{x}_{3}}{2}$, the
Hilbert space of the system is

$L^{2}_{\mathrm{sym}}(\mathbb{R}^{6}):=\left\\{\psi\in
L^{2}(\mathbb{R}^{6})\,\big{|}\;\psi(\bm{x},\bm{y})=\psi(-\bm{x},\bm{y})=\psi\\!\left(\tfrac{1}{2}\,\bm{x}+\bm{y},\tfrac{3}{4}\,\bm{x}-\tfrac{1}{2}\,\bm{y}\right)\right\\}\\!.$
(2.3)

Indeed, notice that the symmetry conditions in (2.3) corresponds to the
exchange of particles $2,3$ and $1,2$ that implies also the condition
$\psi(\bm{x},\bm{y})=\psi\\!\left(\tfrac{1}{2}\,\bm{x}-\bm{y},-\tfrac{3}{4}\,\bm{x}-\tfrac{1}{2}\,\bm{y}\right)$,
associated with the exchange of particles $3,1$. If the bosons interact via
zero-range forces, then the system is described, at least formally, by the
Hamiltonian

$\hat{\mathcal{H}}=\mathcal{H}_{0}+\mu\mspace{1.5mu}\delta(\bm{x})+\mu\mspace{1.5mu}\delta(\tfrac{1}{2}\,\bm{x}+\bm{y})+\mu\mspace{1.5mu}\delta(\tfrac{1}{2}\,\bm{x}-\bm{y})$
(2.4)

where $\mu\\!\in\\!\mathbb{R}$ is a coupling constant and $\mathcal{H}_{0}$ is
the free Hamiltonian of the system, i.e.

$\mathcal{H}_{0}=-\Delta_{\bm{x}}-\tfrac{3}{4}\,\Delta_{\bm{y}}.$ (2.5)

In order to define a rigorous counterpart of $\hat{\mathcal{H}}$, one needs to
build a perturbation of the free Hamiltonian supported on the coincidence
hyperplanes

$\pi_{k}:=\left\\{(\bm{r}_{k},\bm{\rho}_{k})\\!\in\mathbb{R}^{6}\,\big{|}\;\bm{r}_{k}\\!=\bm{0}\right\\},\qquad\pi:=\textstyle{\bigcup_{k=1}^{3}\pi_{k}}$
(2.6a) or, equivalently,
$\begin{split}\pi_{1}:=\\!\left\\{(\bm{x},\bm{y})\\!\in\mspace{-1.5mu}\mathbb{R}^{6}\right.&\left.\\!\big{|}\;\bm{x}=\bm{0}\right\\}\\!,\quad\pi_{2}:=\\!\left\\{(\bm{x},\bm{y})\\!\in\mspace{-1.5mu}\mathbb{R}^{6}\,\big{|}\;\bm{y}=-\tfrac{1}{2}\,\bm{x}\right\\}\\!,\\\
&\pi_{3}:=\\!\left\\{(\bm{x},\bm{y})\\!\in\mspace{-1.5mu}\mathbb{R}^{6}\,\big{|}\;\bm{y}=\tfrac{1}{2}\,\bm{x}\right\\}\\!.\end{split}$
(2.6b)

In other words, we look for a self-adjoint and bounded from below extension in
$L^{2}_{\mathrm{sym}}(\mathbb{R}^{6})$ of the following symmetric and densely
defined operator

$\dot{\mathcal{H}}_{0}:=\mathcal{H}_{0}\big{\rvert}_{\mathscr{D}(\dot{\mathcal{H}}_{0})},\qquad\mathscr{D}(\dot{\mathcal{H}}_{0}):=H^{2}_{0}(\mathbb{R}^{6}\setminus\pi)\cap
L^{2}_{\mathrm{sym}}(\mathbb{R}^{6})$ (2.7)

that is closed according to the graph norm of
$\mathcal{H}_{0},\mathscr{D}(\mathcal{H}_{0})=H^{2}(\mathbb{R}^{6})\cap
L^{2}_{\mathrm{sym}}(\mathbb{R}^{6})$. In particular, we are interested in the
family of self-adjoint extensions studied in [1] (see also [11]) which, at
least formally, are characterized by the boundary condition

$\psi(\bm{x},\bm{y})=\frac{\xi(\bm{y})}{\lvert\bm{x}\rvert}+\alpha(\bm{y})\mspace{1.5mu}\xi(\bm{y})+o(1),\qquad\lvert\bm{x}\rvert\rightarrow
0,$ (2.8)

where $\alpha$ is a position dependent parameter given by

$\alpha:\bm{y}\longmapsto-\frac{1}{\mathfrak{a}}+\frac{\gamma}{\lvert\bm{y}\rvert}\,\theta(\lvert\bm{y}\rvert)$
(2.9)

with $\mathfrak{a}$ the two-body scattering length, $\gamma$ a positive
parameter representing the strength of the regularization and $\theta$ a real
measurable function with compact support such that

$1-\frac{s}{b}\leq\theta(s)\leq 1+\frac{s}{b}\,,\qquad s\geq 0$ (2.10)

for some $b\\!>\mspace{-1.5mu}0$. Notice that assumption (2.10) forces the
function $\theta$ to be continuous at zero, with $\theta(0)=1$. Moreover, the
simplest choice for $\theta$ is the characteristic function of the ball of
radius $b$ centered in the origin. We also stress that, due to the symmetry
constraints of $L^{2}_{\mathrm{sym}}(\mathbb{R}^{6})$, the boundary condition
(2.8) implies

$\displaystyle\psi(\bm{x},\bm{y})=\frac{\xi(\bm{x})}{\left\lvert\bm{y}+\frac{1}{2}\,\bm{x}\right\rvert}+\alpha(\bm{x})\mspace{1.5mu}\xi(\bm{x})+o(1),\;\;\;\qquad\bm{y}\rightarrow-\tfrac{1}{2}\mspace{1.5mu}\bm{x},$
$\displaystyle\psi(\bm{x},\bm{y})=\frac{\xi(-\bm{x})}{\left\lvert\bm{y}-\frac{1}{2}\,\bm{x}\right\rvert}+\alpha(-\bm{x})\mspace{1.5mu}\xi(-\bm{x})+o(1),\qquad\bm{y}\rightarrow\tfrac{1}{2}\mspace{1.5mu}\bm{x}.$

Observe that for $\gamma=0$ equation (2.8) reduces to the standard TMS
boundary condition, which leads to the Thomas effect. Then, for $\gamma>0$ we
are introducing a three-body repulsion meant to regularize the ultraviolet
singularity occurring when the positions of all particles coincide. However,
since $\mathrm{supp}\,\theta$ is a compact, the usual two-body point
interaction is restored when the third particle is far enough.

The procedure adopted in [1] for the rigorous construction of the Hamiltonian
is the following: one first introduces the quadratic form $Q$ in
$L^{2}_{\mathrm{sym}}(\mathbb{R}^{6})$ describing, at least formally, the
expectation value of the energy of our three-body system. Then one defines a
suitable form domain $\mathscr{D}(Q)$ and proves that $Q,\mathscr{D}(Q)$ is
closed and bounded from below. Finally, the Hamiltonian is defined as the
unique self-adjoint and bounded from below operator associated to the
quadratic form.

In order to define the quadratic form $Q,\mathscr{D}(Q)$, we first introduce
an auxiliary hermitian quadratic form $\Phi^{\lambda}$ in
$L^{2}(\mathbb{R}^{3})$ given by [1, equation (3.1)], namely

$\Phi^{\lambda}:=\Phi^{\lambda}_{\mathrm{diag}}\\!+\Phi^{\lambda}_{\mathrm{off}}+\Phi_{\mathrm{reg}}\mspace{-1.5mu}+\Phi_{0},\qquad\mathscr{D}(\Phi^{\lambda})=H^{1/2}(\mathbb{R}^{3}),$
(2.11)

where

$\displaystyle\Phi^{\lambda}_{\mathrm{diag}}[\xi]$
$\displaystyle:=12\pi\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{p}\,\sqrt{\tfrac{3}{4}\mspace{1.5mu}p^{2}+\lambda\mspace{1.5mu}}\,\lvert\hat{\xi}(\bm{p})\rvert^{2},$
(2.12a) $\displaystyle\Phi^{\lambda}_{\mathrm{off}}[\xi]$
$\displaystyle:=-\frac{12}{\pi}\\!\int_{\mathbb{R}^{6}}\mspace{-7.5mu}d\bm{p}\mspace{1.5mu}d\bm{q}\;\frac{\overline{\hat{\xi}(\bm{p})}\,\hat{\xi}(\bm{q})}{p^{2}+q^{2}+\bm{p}\\!\cdot\mspace{-1.5mu}\bm{q}+\lambda},$
(2.12b) $\displaystyle\Phi_{\mathrm{reg}}[\xi]$
$\displaystyle:=\frac{6\gamma}{\pi}\\!\int_{\mathbb{R}^{6}}\mspace{-7.5mu}d\bm{p}\mspace{1.5mu}d\bm{q}\;\frac{\overline{\hat{\xi}(\bm{p})}\,\hat{\xi}(\bm{q})}{|\bm{p}-\bm{q}|^{2}},$
(2.12c) $\displaystyle\Phi_{0}[\xi]$
$\displaystyle:=12\pi\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{y}\>\beta(\bm{y})\lvert\xi(\bm{y})\rvert^{2},\qquad\beta:\bm{y}\longmapsto-\frac{1}{\mathfrak{a}}+\gamma\,\frac{\theta(y)-1}{y}.$
(2.12d)

By assumption (2.10), one has
$\beta\mspace{-1.5mu}\in\mspace{-1.5mu}L^{\infty}(\mathbb{R}^{3})$ and
therefore $\Phi_{0}$ is bounded. The proof of the fact that $\Phi^{\lambda}$
is well defined in $H^{1/2}(\mathbb{R}^{3})$ is relatively standard and it is
given in [1, proposition 3.1]. The more relevant point concerning
$\Phi^{\lambda}$ is that it is coercive for $\lambda$ large enough as long as
$\gamma>\gamma_{c}\mspace{1.5mu}$, with

$\gamma_{c}=\frac{4}{3}-\frac{\sqrt{3}}{\pi}\approx 0.782004.$ (2.13)

The proof is given in [1, proposition 3.6] and it requires a rather long and
non trivial analysis performed in the momentum representation. The conclusion
is that there exists $\lambda_{0}\mspace{-1.5mu}>\mspace{-1.5mu}0$ such that
$\Phi^{\lambda}$ is closed and bounded from below by a positive constant for
each $\lambda\\!>\\!\lambda_{0}$ and $\gamma>\gamma_{c}\mspace{1.5mu}$.
Therefore one can uniquely define a self-adjoint, positive and invertible
operator $\Gamma^{\lambda}$ in $L^{2}(\mathbb{R}^{3})$ such that

$\Phi^{\lambda}[\xi]=\langle\xi,\,\Gamma^{\lambda}\xi\rangle_{L^{2}(\mathbb{R}^{3})},\qquad\forall\,\xi\in
D$ (2.14)

with $D=\mathscr{D}(\Gamma^{\lambda})$ a dense subspace independent of
$\lambda$. Furthermore, defining the continuous111Here
$\mathscr{D}(\mathcal{H}_{0})$ must be intended as a Hilbert subspace of
$L^{2}_{\mathrm{sym}}(\mathbb{R}^{6})$ endowed with the graph norm of
$\mathcal{H}_{0}$. operator

$\begin{split}\tau:\mathscr{D}(\mathcal{H}_{0})\longrightarrow
L^{2}(\mathbb{R}^{3}),\\\ \varphi\longmapsto
12\pi\,\varphi\rvert_{\pi_{1}}\end{split}$ (2.15)

satisfying $\mathrm{ran}(\tau)=H^{1/2}(\mathbb{R}^{3})$ and
$\mathrm{ker}(\tau)=\mathscr{D}(\dot{\mathcal{H}}_{0})$, one can check that
the injective operator $G(z):=(\tau
R_{\mathcal{H}_{0}}(\bar{z}))^{\ast}\\!\in\mathscr{B}(L^{2}(\mathbb{R}^{3}),L^{2}_{\mathrm{sym}}(\mathbb{R}^{6}))$
is represented in the Fourier space by

$(\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[\widthof{G(z)\mspace{1.5mu}\xi}]{\kern
0.1pt\mathchar 866\relax\kern
0.1pt}{\rule{0.0pt}{505.89pt}}}{}}{2.4ex}}\stackon[-6.9pt]{G(z)\mspace{1.5mu}\xi}{\tmpbox})(\bm{k},\bm{p})=\sqrt{\frac{2}{\pi}}\;\frac{\hat{\xi}(\bm{p})+\hat{\xi}\\!\left(\bm{k}\\!-\frac{1}{2}\bm{p}\right)\mspace{-1.5mu}+\hat{\xi}\\!\left(-\bm{k}\\!-\frac{1}{2}\bm{p}\right)\\!}{k^{2}+\frac{3}{4}p^{2}-z},\qquad
z\in\rho(\mathcal{H}_{0}).$ (2.16)

We are now in position to introduce the quadratic form in
$L^{2}_{\mathrm{sym}}(\mathbb{R}^{6})$ ([1, definition 2.1])

$\displaystyle\mathscr{D}(Q)$
$\displaystyle:=\\!\left\\{\psi\mspace{-1.5mu}\in\mspace{-1.5mu}L^{2}_{\mathrm{sym}}(\mathbb{R}^{6})\,\big{|}\;\psi\mspace{-1.5mu}=\phi_{\lambda}\mspace{-1.5mu}+G(-\lambda)\mspace{1.5mu}\xi,\;\phi_{\lambda}\\!\in\mspace{-1.5mu}H^{1}(\mathbb{R}^{6}),\;\xi\mspace{-1.5mu}\in\mspace{-1.5mu}H^{1/2}(\mathbb{R}^{3}),\;\lambda\\!>\mspace{-1.5mu}0\right\\}\\!,$
$\displaystyle Q[\psi]$
$\displaystyle:=\lVert\mathcal{H}_{0}^{1/2}\phi_{\lambda}\rVert^{2}+\lambda\lVert\phi_{\lambda}\rVert^{2}-\lambda\lVert\psi\rVert^{2}+\Phi^{\lambda}[\xi].$
(2.17)

Using the properties of $\Phi^{\lambda}$ and $G(-\lambda)$, it is now easy to
show that the above quadratic form is closed and bounded from below if
$\gamma>\gamma_{c}\mspace{1.5mu}$. Hence it uniquely defines a self-adjoint
and lower semi-bounded operator $\mathcal{H}$ which, by definition, is the
Hamiltonian of the three bosons system.

Following an equivalent approach, one can consider the densely defined and
closed operator $\Gamma(z)$, given by

$\begin{split}&\Gamma(z):D\subset L^{2}(\mathbb{R}^{3})\longrightarrow
L^{2}(\mathbb{R}^{3})\\\
\Gamma(z):=\Gamma^{\lambda}-&(\lambda+z)\mspace{1.5mu}G(\bar{z})^{\ast}G(-\lambda),\qquad\lambda\\!>\\!\lambda_{0},\,z\in\rho(\mathcal{H}_{0})\end{split}$
(2.18)

which represents a sort of analytic continuation of $\Gamma^{\lambda}$, $D$.
Actually, one can prove that $\Gamma(z)$ fulfils

$\displaystyle\Gamma(z)^{\ast}=\Gamma(\bar{z}),\qquad\forall\,z\in\rho(\mathcal{H}_{0}),$
(2.19a)
$\displaystyle\Gamma(z)-\Gamma(w)=(z-w)\mspace{1.5mu}G(\bar{z})^{\ast}G(w),\qquad\forall\,w,z\in\rho(\mathcal{H}_{0}),$
(2.19b)
$\displaystyle\forall\,z\in\mathbb{C}:\;\Re(z)<-\lambda_{0}\>\lor\>\Im(z)>0,\quad
0\in\rho(\Gamma(z)).$ (2.19c)

These properties imply, according to e.g. [14], that for any $z\in\mathbb{C}$
such that $\Gamma(z)$ has a bounded inverse, the operator

$R(z)=(\mathcal{H}_{0}-z)^{-1}+G(z)\Gamma(z)^{-1}G(\bar{z})^{\ast}\,$ (2.20)

defines the resolvent of a self-adjoint and bounded from below operator which
coincides with the Hamiltonian $\mathcal{H}$ obtained with the approach based
on the quadratic form and

$\left\\{z\in\mathbb{C}\,\big{|}\;\Re(z)<-\lambda_{0}\>\lor\>\Im(z)>0\right\\}\\!\subseteq\mspace{-1.5mu}\rho(\mathcal{H}).$

Moreover, one can verify that $\mathcal{H}$ coincides with $\mathcal{H}_{0}$
on $\mathscr{D}(\dot{\mathcal{H}}_{0})$, satisfies the boundary condition
(2.8) in the $L^{2}$ sense (see [1, remark 4.1]) and it is characterized by

$\begin{split}\mathscr{D}(\mathcal{H})&=\\!\left\\{\psi\mspace{-1.5mu}\in\mathscr{D}(Q)\,\big{|}\;\phi_{z}\\!\in\mspace{-1.5mu}\mathscr{D}(\mathcal{H}_{0}),\;\xi\mspace{-1.5mu}\in\mspace{-1.5mu}D,\,\Gamma(z)\mspace{1.5mu}\xi=\tau\phi_{z}\right\\}\\!,\\\
\mathcal{H}\psi&=\mathcal{H}_{0}\mspace{1.5mu}\phi_{z}+z\mspace{1.5mu}G(z)\mspace{1.5mu}\xi.\end{split}$
(2.21)

##  3\. Optimality of $\gamma_{c}$

In this section we prove the optimality of the threshold parameter
$\gamma_{c}$ defined by (2.13). More precisely our goal is to prove the
following theorem.

######  Theorem 3.1.

Whenever $\gamma<\gamma_{c}$, the quadratic form $Q$ given by (2.17) is
unbounded from below.

In order to achieve the result, we shall adapt the ideas contained in [8,
section 5]. Denote for short $G^{\lambda}\\!:=\mspace{-1.5mu}G(-\lambda)$ for
any $\lambda\\!>\\!0$ and let
$\\{u_{n}\\}_{n\in\,\mathbb{N}}\\!\subset\\!\mathscr{D}(Q)$ be a sequence of
trial functions given by

$\displaystyle
u_{n}(\bm{x},\bm{y})=(G^{\lambda}\mspace{1.5mu}\eta_{n})(\bm{x},\bm{y}),$
(3.1)
$\displaystyle\eta_{n}(\bm{y})=n^{2}\mspace{1.5mu}f(n\mspace{1.5mu}\bm{y}),\qquad
f\in H^{1/2}(\mathbb{R}^{3}).$ (3.2)

We stress that, by an explicit estimate due to (2.16), one finds

$\inf_{n\in\mspace{1.5mu}\mathbb{N}}\,\lVert
G^{\lambda}\eta_{n}\rVert_{L^{2}_{\mathrm{sym}}(\mathrm{R}^{6})}\gneq 0.$
(3.3)

Indeed,

$\displaystyle\lVert G^{\lambda}\eta_{n}\rVert^{2}$
$\displaystyle=\frac{2}{\pi}\\!\int_{\mathbb{R}^{6}}\mspace{-7.5mu}d\bm{k}d\bm{p}\>\frac{1}{n^{2}}\frac{\left\lvert\hat{f}\\!\left(\tfrac{\bm{p}}{n}\right)+\hat{f}\mspace{-1.5mu}\left(\tfrac{\bm{k}}{n}-\tfrac{\bm{p}}{2n}\right)+\hat{f}\mspace{-1.5mu}\left(\tfrac{\bm{k}}{n}+\tfrac{\bm{p}}{2n}\right)\right\rvert^{2}}{\left(k^{2}+\frac{3}{4}p^{2}+\lambda\right)^{2}}$
$\displaystyle=\frac{2}{\pi}\\!\int_{\mathbb{R}^{6}}\mspace{-7.5mu}d\bm{\kappa}d\mspace{-1.5mu}\bm{q}\;\frac{\left\lvert\hat{f}(\bm{q})+\hat{f}\mspace{-1.5mu}\left(\bm{\kappa}-\tfrac{\bm{q}}{2}\right)+\hat{f}\mspace{-1.5mu}\left(\bm{\kappa}+\tfrac{\bm{q}}{2}\right)\right\rvert^{2}}{\left(\kappa^{2}+\frac{3}{4}q^{2}+\frac{\lambda}{n^{2}}\right)^{2}}>\lVert
G^{\lambda}f\rVert^{2}.$

Our goal is to show that whenever $\gamma$ is smaller than the threshold value
$\gamma_{c}$ given by (2.13), one has

$\lim_{n\rightarrow\,+\infty}Q[u_{n}]=-\infty,$ (3.4)

According to (2.17), we have

$Q[u_{n}]=-\lambda\,\lVert
G^{\lambda}\mspace{1.5mu}\eta_{n}\rVert_{L^{2}_{\mathrm{sym}}(\mathrm{R}^{6})}+\Phi^{\lambda}[\eta_{n}]\leq-\lambda\,\lVert
G^{\lambda}f\rVert_{L^{2}_{\mathrm{sym}}(\mathrm{R}^{6})}+\Phi^{\lambda}[\eta_{n}]$
(3.5)

and therefore the theorem is proven if we exhibit some
$f\\!\in\mspace{-1.5mu}H^{1/2}(\mathbb{R}^{3})$ such that

$\lim_{n\rightarrow\,+\infty}\Phi^{\lambda}[\eta_{n}]=-\infty.$ (3.6)

######  Lemma 3.2.

Let $\Phi^{\lambda}$ and $\eta_{n}$ be given by (2.11) and (3.2),
respectively. Then, one has

$\Phi^{\lambda}[\eta_{n}]=n^{2}(\Phi^{0}_{\mathrm{diag}}\\!+\Phi^{0}_{\mathrm{off}}\mspace{-1.5mu}+\Phi_{\mathrm{reg}})[f]+\mathcal{O}(n).$

###### Proof.

First of all, we can neglect the bounded component $\Phi_{0}\mspace{1.5mu}$,
since

$\displaystyle\Phi_{0}[\eta_{n}]$ $\displaystyle=12\pi
n^{4}\\!\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{y}\>\beta(y)\,\lvert
f(n\mspace{1.5mu}\bm{y})\rvert^{2}=12\pi\mspace{1.5mu}n\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{t}\>\beta\\!\left(\tfrac{t}{n}\right)\lvert
f(\bm{t})\rvert^{2}$ $\displaystyle\leq
12\pi\mspace{1.5mu}n\,\lVert\beta\rVert_{L^{\infty}(\mathbb{R}^{3})}\mspace{1.5mu}\lVert
f\rVert^{2}_{L^{2}(\mathbb{R}^{3})}\implies\Phi_{0}[\eta_{n}]=\mathcal{O}(n),\qquad
n\rightarrow+\infty.$

Next, rescaling properly the variables in computing
$\Phi^{\lambda}_{\mathrm{diag}}[\eta_{n}]$, one gets

$\displaystyle\Phi^{\lambda}_{\mathrm{diag}}$
$\displaystyle[\eta_{n}]=12\pi\mspace{1.5mu}n\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{\kappa}\>\sqrt{\tfrac{3}{4}\,n^{2}\kappa^{2}+\lambda\,}\,\lvert\hat{f}(\bm{\kappa})\rvert^{2}$
$\displaystyle=6\sqrt{3}\,\pi\mspace{1.5mu}n^{2}\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{\kappa}\>\lvert\bm{\kappa}\rvert\,\lvert\hat{f}(\bm{\kappa})\rvert^{2}+12\pi\mspace{1.5mu}n\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{\kappa}\mspace{1.5mu}\left(\mspace{-1.5mu}\sqrt{\tfrac{3}{4}\,n^{2}\kappa^{2}+\lambda\,}-\sqrt{\tfrac{3}{4}}\,n\kappa\right)\mspace{-1.5mu}\lvert\hat{f}(\bm{\kappa})\rvert^{2}$
$\displaystyle=n^{2}\Phi^{\mspace{1.5mu}0}_{\mathrm{diag}}[f]+o(n).$

Indeed, exploiting the elementary inequality $\sqrt{a^{2}+b^{2}}-\lvert
a\rvert\leq\lvert b\rvert$, we can use the dominated convergence theorem to
obtain

$\lim_{n\rightarrow+\infty}\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{\kappa}\left(\mspace{-1.5mu}\sqrt{\tfrac{3}{4}\,n^{2}\kappa^{2}+\lambda\,}-\sqrt{\tfrac{3}{4}}\,n\kappa\right)\mspace{-1.5mu}\lvert\hat{f}(\bm{\kappa})\rvert^{2}=0.$

Concerning the regularizing contribution, one simply has

$\displaystyle\Phi_{\mathrm{reg}}[\eta_{n}]$
$\displaystyle=\frac{6\gamma}{\pi}\mspace{1.5mu}n^{2}\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{p}\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{q}\;\frac{\overline{\hat{f}(\bm{p})}\,\hat{f}(\bm{q})}{\lvert\bm{p}-\bm{q}\rvert^{2}}$
$\displaystyle=n^{2}\Phi_{\mathrm{reg}}[f].$

Finally, we compute $\Phi^{\lambda}_{\mathrm{off}}[\eta_{n}]$

$\displaystyle\Phi^{\lambda}_{\mathrm{off}}$
$\displaystyle[\eta_{n}]=-\frac{12}{\pi}\mspace{1.5mu}n^{2}\\!\int_{\mathbb{R}^{6}}\mspace{-7.5mu}d\bm{p}\mspace{1.5mu}d\bm{q}\;\frac{\overline{\hat{f}(\bm{p})}\,\hat{f}(\bm{q})}{p^{2}+q^{2}+\bm{p}\\!\cdot\mspace{-1.5mu}\bm{q}+\frac{\lambda}{n^{2}}}$
$\displaystyle=$
$\displaystyle-\frac{12}{\pi}\mspace{1.5mu}n^{2}\\!\int_{\mathbb{R}^{6}}\mspace{-7.5mu}d\bm{p}\mspace{1.5mu}d\bm{q}\;\frac{\overline{\hat{f}(\bm{p})}\,\hat{f}(\bm{q})}{p^{2}+q^{2}+\bm{p}\\!\cdot\mspace{-1.5mu}\bm{q}}\,+$
$\displaystyle+\frac{12}{\pi}\mspace{1.5mu}\lambda\\!\int_{\mathbb{R}^{6}}\mspace{-7.5mu}d\bm{p}\mspace{1.5mu}d\bm{q}\;\frac{\overline{\hat{f}(\bm{p})}\,\hat{f}(\bm{q})}{\big{(}p^{2}+q^{2}+\bm{p}\\!\cdot\mspace{-1.5mu}\bm{q}+\frac{\lambda}{n^{2}}\big{)}\big{(}p^{2}+q^{2}+\bm{p}\\!\cdot\mspace{-1.5mu}\bm{q}\big{)}}.$

Defining the integral operator in $L^{2}(\mathbb{R}^{3})$ given by

$(P_{\\!n}\mspace{1.5mu}\hat{\varphi})(\bm{p}):=\frac{12}{\pi}\mspace{1.5mu}\lambda\\!\int_{\mathrm{R}^{3}}\mspace{-7.5mu}d\bm{q}\;\frac{\hat{\varphi}(\bm{q})}{\big{(}p^{2}+q^{2}+\bm{p}\\!\cdot\mspace{-1.5mu}\bm{q}+\frac{\lambda}{n^{2}}\big{)}\big{(}p^{2}+q^{2}+\bm{p}\\!\cdot\mspace{-1.5mu}\bm{q}\big{)}},$
(3.7)

we can write

$\Phi^{\lambda}_{\mathrm{off}}[\eta_{n}]=n^{2}\Phi^{\mspace{1.5mu}0}_{\mathrm{off}}[f]+\\!\int_{\mathrm{R}^{3}}\mspace{-7.5mu}d\bm{p}\>\overline{\hat{f}(\bm{p})}\,(P_{\\!n}\mspace{1.5mu}\hat{f})(\bm{p}).$

We notice that $P_{\\!n}$ is a Hilbert-Schmidt operator and

$\displaystyle\lVert
P_{\\!n}\rVert^{2}_{\mathscr{B}(L^{2}(\mathbb{R}^{3}))}\\!$
$\displaystyle\leq\frac{124\mspace{1.5mu}\lambda^{2}\\!}{\pi^{2}}\\!\int_{\mathbb{R}^{6}}\mspace{-7.5mu}d\bm{p}\mspace{1.5mu}d\bm{q}\;\frac{1}{\big{(}p^{2}+q^{2}+\bm{p}\\!\cdot\mspace{-1.5mu}\bm{q}+\frac{\lambda}{n^{2}}\big{)}^{2}\big{(}p^{2}+q^{2}+\bm{p}\\!\cdot\mspace{-1.5mu}\bm{q}\big{)}^{2}}$
$\displaystyle\mspace{-18.0mu}\leq\frac{124\mspace{1.5mu}\lambda^{2}\\!}{\pi^{2}}\\!\int_{\mathbb{R}^{6}}\mspace{-7.5mu}d\bm{p}\mspace{1.5mu}d\bm{q}\;\frac{4}{\big{(}\frac{p^{2}+\,q^{2}\\!}{2}+\frac{\lambda}{n^{2}}\big{)}^{2}\big{(}p^{2}+q^{2}\big{)}^{2}}=496\mspace{1.5mu}\pi\lambda^{2}\\!\int_{0}^{+\infty}\mspace{-24.0mu}dk\;\frac{k}{\big{(}\frac{k^{2}\\!}{2}+\frac{\lambda}{n^{2}}\big{)}^{2}\\!}$
$\displaystyle\mspace{-18.0mu}=496\mspace{1.5mu}\pi\lambda\mspace{1.5mu}n^{2}.$

Using the above estimate, we find

$\Phi^{\lambda}_{\mathrm{off}}[\eta_{n}]=n^{2}\Phi^{\mspace{1.5mu}0}_{\mathrm{off}}[f]+O(n)$

and the lemma is proven.

∎

In light of lemma 3.2, it is straightforward to see that (3.6) is achieved as
soon as we exhibit a function $f\\!\in\mspace{-1.5mu}H^{1/2}(\mathbb{R}^{3})$
such that, whenever $\gamma\\!<\mspace{-1.5mu}\gamma_{c}\mspace{1.5mu}$, holds

$\Phi^{\mspace{1.5mu}0}_{\mathrm{diag}}[f]+\Phi^{\mspace{1.5mu}0}_{\mathrm{off}}[f]+\Phi_{\mathrm{reg}}[f]<0.$

A relevant feature of the previous lemma is that the leading order of
$\Phi^{\lambda}[\eta_{n}]$ as $n$ goes to infinity does not depend on
$\lambda$ and, therefore, we have reduced the problem to the study of the
hermitian quadratic form evaluated in $\lambda=0$ which is diagonalizable.

In the next lemma we exhibit the trial function we need to prove our result.

######  Lemma 3.3.

Let $\gamma_{c}$ be defined by (2.13), assume
$\gamma\\!<\mspace{-1.5mu}\gamma_{c}$ and let us consider the family of trial
functions $f_{\beta}\\!\in\mspace{-1.5mu}H^{1/2}(\mathbb{R}^{3})$ such that

$\hat{f}_{\beta}(\bm{p})=\tfrac{1}{p^{2}}\exp\\!\left(-\tfrac{p^{\beta}\mspace{1.5mu}+\,p^{-\beta}}{\mspace{-12.0mu}2}\mspace{-1.5mu}\right)\\!,\qquad\beta>0.$

Then there exists $\beta_{0}>0$ such that for any $\beta\in(0,\beta_{0})$ we
have

$\big{(}\Phi^{\mspace{1.5mu}0}_{\mathrm{diag}}\\!+\Phi^{\mspace{1.5mu}0}_{\mathrm{off}}\mspace{-1.5mu}+\Phi_{\mathrm{reg}}\big{)}[f_{\beta}]\mspace{-1.5mu}<0.$

###### Proof.

We stress that our trial functions are entirely lying in the $s$-wave
subspace, therefore we have

$\displaystyle\Phi^{0}_{\mathrm{diag}}[f_{\beta}]$
$\displaystyle:=48\pi^{2}\sqrt{\tfrac{3}{4}}\\!\int_{0}^{+\infty}\mspace{-24.0mu}d\mspace{-1.5mu}p\>p^{3}\,\lvert\hat{f}_{\beta}(p)\rvert^{2},$
(3.8a) $\displaystyle\Phi^{0}_{\mathrm{off}}[f_{\beta}]$
$\displaystyle:=-96\pi\\!\int_{0}^{+\infty}\mspace{-24.0mu}d\mspace{-1.5mu}p\>p\\!\int_{0}^{+\infty}\mspace{-24.0mu}dq\>q\,\overline{\hat{f}_{\beta}(p)}\,\hat{f}_{\beta}(q)\ln\mspace{-1.5mu}\left(\frac{p^{2}+q^{2}+pq}{p^{2}+q^{2}-pq}\right)\\!,$
(3.8b) $\displaystyle\Phi_{\mathrm{reg}}[f_{\beta}]$
$\displaystyle:=24\pi\gamma\\!\\!\int_{0}^{+\infty}\mspace{-24.0mu}d\mspace{-1.5mu}p\>p\\!\int_{0}^{+\infty}\mspace{-24.0mu}dq\>q\,\overline{\hat{f}_{\beta}(p)}\,\hat{f}_{\beta}(q)\ln\mspace{-1.5mu}\left(\frac{p^{2}+q^{2}+2pq}{p^{2}+q^{2}-2pq}\right)\\!,$
(3.8c)

where we have used the identity222Equation (3.9) is an application of the
addition formula for the spherical harmonics in the $s$-wave.

$\int_{\mathbb{R}^{6}}\mspace{-7.5mu}d\bm{p}\mspace{1.5mu}d\bm{q}\>g\\!\left(p,q,\tfrac{\bm{p}\,\cdot\mspace{1.5mu}\bm{q}}{p\mspace{1.5mu}q}\right)=8\pi^{2}\\!\\!\int_{0}^{+\infty}\mspace{-24.0mu}d\mspace{-1.5mu}p\>p^{2}\\!\int_{0}^{+\infty}\mspace{-24.0mu}dq\>q^{2}\\!\int_{-1}^{1}\mspace{-7.5mu}du\>g(p,q,u)$
(3.9)

holding for any integrable function
$g:\mathbb{R}^{2}_{+}\times[-1,1]\longrightarrow\mathbb{C}$. According to,
e.g. [1, lemma 3.4], the quantities in equations (3.8) can be diagonalized
through the unitary transformation

$\begin{split}\mathcal{M}:L^{2}(\mathbb{R}_{+},\,p^{2}\textstyle{\sqrt{p^{2}+1}}\,dp)\longrightarrow
L^{2}(\mathbb{R}),\\\
\psi\longmapsto\psi^{\sharp}(x)=\frac{1}{\mspace{-4.5mu}\sqrt{\mspace{1.5mu}2\pi}\,}\\!\int_{\mathbb{R}}\mspace{-3.0mu}dt\>e^{-i\mspace{1.5mu}tx}e^{2t}\psi(e^{t})\end{split}$
(3.10)

yielding (see [1, lemmata 3.4, 3.5])

$\displaystyle\Phi^{\mspace{1.5mu}0}_{\mathrm{diag}}[f_{\beta}]=48\mspace{1.5mu}\sqrt{\tfrac{3}{4}}\,\pi^{2}\\!\int_{\mathbb{R}}\\!dx\>\lvert\hat{f}^{\,\sharp}_{\beta}(x)\rvert^{2},$
(3.11a)
$\displaystyle\Phi^{\mspace{1.5mu}0}_{\mathrm{off}}[f_{\beta}]=-48\pi^{2}\\!\mspace{-1.5mu}\int_{\mathbb{R}}\\!dx\>|\hat{f}^{\,\sharp}_{\beta}(x)\rvert^{2}\,\frac{4\sinh\\!\left(\frac{\pi}{6}x\right)}{x\cosh\\!\left(\frac{\pi}{2}x\right)},$
(3.11b)
$\displaystyle\Phi_{\mathrm{reg}}[f_{\beta}]=48\pi^{2}\\!\int_{\mathbb{R}}\\!dx\>\lvert\hat{f}^{\,\sharp}_{\beta}(x)\rvert^{2}\,\frac{\gamma\tanh\\!\left(\frac{\pi}{2}x\right)}{x}.$
(3.11c)

Let us introduce the bounded and continuous function

$S(x):=\frac{\sqrt{3}}{2}+\frac{\gamma\sinh\\!\left(\frac{\pi}{2}x\right)-4\sinh\\!\left(\tfrac{\pi}{6}x\right)}{x\mspace{1.5mu}\cosh\\!\left(\frac{\pi}{2}x\right)}$
(3.12)

so that we have

$\displaystyle\big{(}\Phi^{\mspace{1.5mu}0}_{\mathrm{diag}}\\!+\Phi^{\mspace{1.5mu}0}_{\mathrm{off}}\mspace{-1.5mu}+\Phi_{\mathrm{reg}}\big{)}[f_{\beta}]=48\pi^{2}\\!\\!\int_{\mathbb{R}}\\!dx\>\lvert\hat{f}^{\,\sharp}_{\beta}(x)\rvert^{2}\,S(x),$
(3.13)

with

$\lim_{x\rightarrow\,0}S(x)=\tfrac{\sqrt{3}}{2}-\frac{2\pi}{3}+\frac{\pi}{2}\gamma=\frac{\pi}{2}(\gamma-\gamma_{c})<0.$
(3.14)

Roughly speaking, the integral in (3.13) is negative if we choose the trial
function such that the support of $\hat{f}^{\,\sharp}_{\beta}$ is sufficiently
concentrated in a neighborhood of zero. More precisely, considering the
explicit expression of $\hat{f}_{\beta}\mspace{1.5mu}$, we have333 We stress
that
$\hat{f}^{\,\sharp}_{\beta}(x)=\sqrt{\tfrac{2}{\pi}}\,\tfrac{1}{\beta}\,\mathrm{K}_{\mspace{1.5mu}i\mspace{1.5mu}x/\beta}(1)$
because of the integral representation for the Macdonald function
$\mathrm{K}_{\nu}$ given by [9, p. 384, 3.547 4].

$\hat{f}^{\,\sharp}_{\beta}(x)=\frac{1}{\beta}\,\hat{h}\\!\left(\tfrac{x}{\beta}\right)\\!$
(3.15)

where $h(p)=e^{-\cosh{p}}$ with $h\in\mathcal{S}(\mathbb{R})$. Then

$\displaystyle\int_{\mathbb{R}}\\!dx\>\lvert\hat{f}^{\,\sharp}_{\beta}(x)\rvert^{2}S(x)=\frac{1}{\beta^{2}}\\!\int_{\mathbb{R}}\\!dx\>\lvert\hat{h}(x/\beta)\rvert^{2}\,S(x)$
$\displaystyle=\frac{1}{\beta}\\!\int_{\mathbb{R}}\\!dx\>\lvert\hat{h}(x)\rvert^{2}\,S(\beta\mspace{1.5mu}x).$

By dominated convergence theorem we obtain

$\lim_{\beta\rightarrow\,0^{+}}\int_{\mathbb{R}}\\!dx\>\lvert\hat{h}(x)\rvert^{2}\,S(\beta\mspace{1.5mu}x)=\lVert
h\rVert^{2}_{L^{2}(\mathbb{R})}\,S(0)<0.$

Hence, the lemma is proven by noticing that the previous integral is
continuous in $\beta\\!>\\!0$ and therefore, the quadratic form
$\big{(}\Phi^{\mspace{1.5mu}0}_{\mathrm{diag}}\\!+\Phi^{\mspace{1.5mu}0}_{\mathrm{off}}\mspace{-1.5mu}+\Phi_{\mathrm{reg}}\big{)}[f_{\beta}]$
is negative for any $\beta$ small enough.

∎

###### Proof of theorem 3.1.

Let $\hat{f}_{\beta}$ be the trial function given in lemma 3.3 with
$\beta<\beta_{0}$ and consider the following sequence of charges

$\hat{\eta}^{\,\beta}_{n}(\bm{p})=\frac{1}{n}\mspace{1.5mu}\hat{f}_{\beta}\\!\left(\frac{\bm{p}}{n}\right).$

By lemma 3.2, we know that

$\Phi^{\lambda}[\eta^{\,\beta}_{n}]=n^{2}(\Phi^{0}_{\mathrm{diag}}\\!+\Phi^{0}_{\mathrm{off}}\mspace{-1.5mu}+\Phi_{\mathrm{reg}})[f_{\beta}]+\mathcal{O}(n),\qquad
n\rightarrow+\infty$

and then $\Phi^{\lambda}[\eta^{\,\beta}_{n}]\longrightarrow-\infty$ as $n$
grows to infinity.

∎

##  4\. Analysis in Position Space

In this section, we give a different proof of the coercivity of
$\Phi^{\lambda}$ based on the representation of $\Phi^{\lambda}$ in position
space. In particular, this approach allows to identify the negative
contribution of the quadratic form $\Phi^{\lambda}_{\mathrm{off}}$ and
therefore to justify the choice of the regularization term
$\Phi_{\mathrm{reg}}\mspace{1.5mu}$.

In the next proposition, we write the quadratic form $\Phi^{\lambda}$ defined
in (2.11) in the position-space representation.

######  Proposition 4.1.

For any $\xi\in H^{1/2}(\mathbb{R}^{3})$ and $\lambda>0$ one has

$\displaystyle\Phi^{\lambda}_{\mathrm{diag}}[\xi]=12\pi\mspace{1.5mu}\sqrt{\lambda}\,\lVert\xi\rVert^{2}+\frac{2\sqrt{3}\mspace{1.5mu}\lambda}{\pi}\\!\int_{\mathbb{R}^{6}}\mspace{-7.5mu}d\bm{x}d\bm{y}\;\frac{\lvert\xi(\bm{x})-\xi(\bm{y})\rvert^{2}\mspace{-4.5mu}}{\lvert\bm{x}-\bm{y}\rvert^{2}}\,\mathrm{K}_{2}\\!\left(\sqrt{\tfrac{4\lambda}{3}}\mspace{1.5mu}\lvert\bm{x}-\bm{y}\rvert\right)\\!,$
(4.1a)
$\displaystyle\Phi^{\lambda}_{\mathrm{off}}[\xi]=-\frac{8\sqrt{3}\mspace{1.5mu}\lambda}{\pi}\\!\int_{\mathbb{R}^{6}}\mspace{-7.5mu}d\bm{x}d\bm{y}\;\frac{\overline{\xi(\bm{x})}\,\xi(\bm{y})}{y^{2}+x^{2}+\bm{x}\mspace{-1.5mu}\cdot\mspace{-1.5mu}\bm{y}}\,\mathrm{K}_{2}\\!\left(\sqrt{\tfrac{4\lambda}{3}}\,\textstyle{\sqrt{y^{2}+x^{2}+\bm{x}\mspace{-1.5mu}\cdot\mspace{-1.5mu}\bm{y}}}\right)\\!,$
(4.1b)
$\displaystyle\Phi_{\mathrm{reg}}[\xi]=12\pi\mspace{1.5mu}\gamma\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{x}\;\frac{\lvert\xi(\bm{x})\rvert^{2}\mspace{-4.5mu}}{\lvert\bm{x}\rvert}$
(4.1c)

where $\mathrm{K}_{\nu}:\mathbb{R}_{+}\longrightarrow\mathbb{C}$ is the
modified Bessel function of the second kind (also known as Macdonald function)
and order $\nu\in\mathbb{C}$.

###### Proof.

Identity (4.1a) is a consequence of (2.12a) and [10, section 7.12, (5)], while
(4.1c) is obtained by comparing (2.12c) with the identity

$\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{r}\;\frac{\lvert
f(\bm{r})\rvert^{2}\\!}{r}=\frac{1}{2\pi^{2}\\!}\int_{\mathbb{R}^{6}}\mspace{-7.5mu}d\bm{p}\mspace{1.5mu}d\bm{q}\;\frac{\overline{\hat{f}(\bm{p})}\,\hat{f}(\bm{q})}{\lvert\bm{p}-\bm{q}\rvert^{2}},\qquad\forall\,f\in
H^{1/2}(\mathbb{R}^{3}).$ (4.2)

Concerning the proof of (4.1b), we consider (2.12b) for
$\varphi\\-\in\\-\mathcal{S}(\mathbb{R}^{3})$ and we observe that we have
uniformly in $\lambda>0$

$\bm{\sigma}\longmapsto\frac{1}{\sigma^{2}+\tau^{2}+\bm{\tau}\\!\cdot\mspace{-1.5mu}\bm{\sigma}+\lambda}\in
L^{2}(\mathbb{R}^{3},d\bm{\sigma}),\qquad\text{for }\tau\neq 0.$

Therefore, by Plancherel’s theorem we find

$\Phi^{\lambda}_{\mathrm{off}}[\varphi]=-\frac{12}{\pi}\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{\tau}\>\overline{\hat{\varphi}(\bm{\tau})}\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{x}\;\frac{\varphi(\bm{x})}{(2\pi)^{3/2}\\!}\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{\sigma}\;\frac{e^{-i\,\bm{x}\cdot\mspace{1.5mu}\bm{\sigma}}}{\tau^{2}+\sigma^{2}+\bm{\tau}\\!\cdot\mspace{-1.5mu}\bm{\sigma}+\lambda}.$

Using the change of coordinates $\bm{\sigma}=\bm{q}-\frac{\bm{\tau}}{2}$, we
obtain

$\displaystyle\Phi^{\lambda}_{\mathrm{off}}[\varphi]$
$\displaystyle=-\frac{12}{\pi}\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{\tau}\>\overline{\hat{\varphi}(\bm{\tau})}\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{x}\;\frac{\>e^{\frac{\bm{\tau}\mspace{1.5mu}\cdot\,\bm{x}}{2}i}\\!}{(2\pi)^{3/2}\\!}\>\varphi(\bm{x})\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{q}\;\frac{e^{-i\,\bm{q}\cdot\bm{x}}}{\frac{3}{4}\tau^{2}+q^{2}+\lambda}$
$\displaystyle=-\frac{24\pi}{\>(2\pi)^{3/2}\\!}\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{\tau}\>\overline{\hat{\varphi}(\bm{\tau})}\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{x}\;\frac{\varphi(\bm{x})}{\lvert\bm{x}\rvert}\,e^{\frac{\bm{\tau}\mspace{1.5mu}\cdot\,\bm{x}}{2}i\,-\sqrt{\frac{3}{4}\tau^{2}\,+\,\lambda}\>\lvert\bm{x}\rvert}$
$\displaystyle=-\frac{24\pi}{\>(2\pi)^{3/2}\\!}\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{x}\;\frac{\varphi(\bm{x})}{\lvert\bm{x}\rvert}\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{\tau}\>\overline{\hat{\varphi}(\bm{\tau})}\,e^{\frac{\bm{\tau}\mspace{1.5mu}\cdot\,\bm{x}}{2}i\,-\sqrt{\frac{3}{4}\tau^{2}\,+\,\lambda}\>\lvert\bm{x}\rvert}.$

Since uniformly in $\lambda>0$

$\bm{\tau}\longmapsto
e^{\frac{\bm{\tau}\mspace{1.5mu}\cdot\,\bm{x}}{2}\,i\,-\sqrt{\frac{3}{4}\tau^{2}\,+\,\lambda}\>\lvert\bm{x}\rvert}\in
L^{2}(\mathbb{R}^{3},d\bm{\tau}),\qquad\text{for }x\neq 0$

we use again Plancherel’s theorem to obtain

$\displaystyle\Phi^{\lambda}_{\mathrm{off}}[\varphi]$
$\displaystyle=-\frac{24\pi}{(2\pi)^{3}}\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{x}\;\frac{\varphi(\bm{x})}{\lvert\bm{x}\rvert}\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{y}\>\overline{\varphi(\bm{y})}\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{\tau}\>e^{i\mspace{1.5mu}\bm{\tau}\cdot\left(\bm{y}\mspace{1.5mu}+\mspace{1.5mu}\frac{\bm{x}}{2}\\!\right)\,-\sqrt{\frac{3}{4}\tau^{2}\,+\,\lambda}\>\lvert\bm{x}\rvert}$
$\displaystyle=-\frac{12}{\pi}\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{x}\;\frac{\varphi(\bm{x})}{\lvert\bm{x}\rvert}\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{y}\;\frac{\overline{\varphi(\bm{y})}}{\lvert\bm{y}+\frac{\bm{x}}{2}\rvert}\int_{0}^{+\infty}\mspace{-24.0mu}d\tau\>\tau\sin\\!\left(\tau\lvert\bm{y}+\tfrac{\bm{x}}{2}\rvert\right)e^{-\sqrt{\frac{3}{4}\tau^{2}\,+\,\lambda}\>\lvert\bm{x}\rvert}.$

The last integral can be explicitly computed using the formula (see, e.g., [9,
p. 491, 3.914.6])

$\int_{0}^{+\infty}\mspace{-24.0mu}dx\>x\mspace{1.5mu}\sin(b\mspace{1.5mu}x)\,e^{-\beta\mspace{1.5mu}\sqrt{x^{2}\,+\,\gamma^{2}}}\\!=\frac{b\mspace{1.5mu}\beta\mspace{1.5mu}\gamma^{2}}{\beta^{2}\mspace{-1.5mu}+b^{2}}\,\mathrm{K}_{2}\left(\gamma\sqrt{\beta^{2}\mspace{-1.5mu}+b^{2}}\right)\\!,\quad\forall\,b\in\mathbb{R}\mspace{1.5mu}\text{
and }\beta,\gamma>0$ (4.3)

and therefore identity (4.1b) is proven for
$\varphi\\!\in\mspace{-1.5mu}\mathcal{S}(\mathbb{R}^{3})$. By a density
argument444Because of propositions 4.2 and 4.3, one can obtain a control in
the $H^{1/2}$ norm. the result is extended to any
$\varphi\mspace{-1.5mu}\in\\!H^{1/2}(\mathbb{R}^{3})$.

∎

Before proceeding, let us briefly recall some elementary properties of
$\mathrm{K}_{2}(\cdot)$:

$\displaystyle x^{2}\mspace{1.5mu}\mathrm{K}_{2}(x)\text{ is decreasing in
}x>0,$ (4.4a)
$\displaystyle\mathrm{K}_{2}(x)=\sqrt{\frac{\pi}{2}}\,e^{-x}\left[\frac{1}{\\!\sqrt{x}\>}+\mathcal{O}\\!\left(\frac{1}{x^{3/2}\\!}\right)\mspace{-1.5mu}\right]\\!,\qquad\text{as
}x\rightarrow+\infty,$ (4.4b)
$\displaystyle\mathrm{K}_{2}(x)=\frac{2}{\,x^{2}\\!}-\frac{1}{2}+\mathcal{O}\\!\left(x^{2}\ln{x}\right)\\!,\qquad\text{as
}x\rightarrow 0^{+}\\!.$ (4.4c) In particular, notice that (4.4a) and (4.4c)
imply $\mathrm{K}_{2}(x)\leq\frac{2}{\,x^{2}\\!}\,,\qquad\forall\,x>0.$ (4.4d)

In the next proposition we show the relevant fact that the negative
contribution of $\Phi^{\lambda}_{\mathrm{off}}$ can be explicitly
characterized.

######  Proposition 4.2.

For any $\varphi\in H^{1/2}(\mathbb{R}^{3})$ and $\lambda>0$ one has

$\begin{split}\Phi^{\lambda}_{\mathrm{off}}[\varphi]=&\,-24\pi\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{x}\;\frac{\;e^{-\sqrt{\lambda}\,\lvert\bm{x}\rvert}\mspace{-7.5mu}}{\lvert\bm{x}\rvert}\>\lvert\varphi(\bm{x})\rvert^{2}+\\\
&+\frac{4\sqrt{3}\mspace{1.5mu}\lambda}{\pi}\\!\int_{\mathbb{R}^{6}}\mspace{-7.5mu}d\bm{x}d\bm{y}\;\frac{\lvert\varphi(\bm{x})-\varphi(\bm{y})\rvert^{2}}{y^{2}+x^{2}+\bm{x}\mspace{-1.5mu}\cdot\mspace{-1.5mu}\bm{y}}\,\mathrm{K}_{2}\\!\left(\sqrt{\tfrac{4\lambda}{3}}\,\textstyle{\sqrt{y^{2}+x^{2}+\bm{x}\mspace{-1.5mu}\cdot\mspace{-1.5mu}\bm{y}}}\right)\\!.\end{split}$
(4.5)

###### Proof.

Let us decompose the expression given by (4.1b) as follows

$\begin{split}\Phi^{\lambda}_{\mathrm{off}}[\varphi]=-\frac{8\sqrt{3}\mspace{1.5mu}\lambda}{\pi}\,\Bigg{[}\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{y}\>\lvert\varphi(\bm{y})\rvert^{2}\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{x}\;\frac{\mathrm{K}_{2}\Big{(}\sqrt{\tfrac{4\lambda}{3}}\,\textstyle{\sqrt{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}}\Big{)}}{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}\,+\\\
+\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{y}\>\overline{\varphi(\bm{y})}\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{x}\;\frac{\varphi(\bm{x})-\varphi(\bm{y})}{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}\,\mathrm{K}_{2}\\!\left(\sqrt{\tfrac{4\lambda}{3}}\,\textstyle{\sqrt{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}}\right)\\!\Bigg{]},\end{split}$

Then, we evaluate the first term in the right hand side. In proposition 4.1 we
have seen that the function

$\begin{split}&\hat{f}^{\lambda}_{x}:\mathbb{R}^{3}\longrightarrow\mathbb{R}_{+},\qquad
x,\lambda\in\mathbb{R}_{+},\\\
&\mspace{27.0mu}\bm{\tau}\longmapsto\frac{e^{-x\mspace{1.5mu}\sqrt{\frac{3}{4}\tau^{2}\,+\,\lambda}}}{x}\end{split}$
(4.6)

is such that

$f^{\lambda}_{\lvert\bm{x}\rvert}\\!\left(\bm{y}+\frac{\bm{x}}{2}\right)=\sqrt{\frac{8}{3\pi}}\,\lambda\,\frac{\mathrm{K}_{2}\Big{(}\sqrt{\frac{4\lambda}{3}}\,\textstyle{\sqrt{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}}\,\Big{)}}{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}.$
(4.7)

Notice the symmetry in the exchange $\bm{x}\longleftrightarrow\bm{y}$. Then,

$\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{x}\>f^{\lambda}_{\lvert\bm{x}\rvert}\\!\left(\bm{y}+\tfrac{\bm{x}}{2}\right)=\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{x}\>f^{\lambda}_{\lvert\bm{y}\rvert}\\!\left(\bm{x}+\tfrac{\bm{y}}{2}\right)=\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{z}\>f^{\lambda}_{\lvert\bm{y}\rvert}(\bm{z})=(2\pi)^{3/2}\hat{f}^{\lambda}_{\lvert\bm{y}\rvert}(0).$

Therefore we find

$\frac{\lambda}{\\!\sqrt{3}\pi^{2}}\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{x}\;\frac{\mathrm{K}_{2}\Big{(}\sqrt{\tfrac{4\lambda}{3}}\,\textstyle{\sqrt{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}}\,\Big{)}}{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}=\frac{\>e^{-\sqrt{\lambda}\mspace{1.5mu}\lvert\bm{y}\rvert}\\!\\!\\!}{\lvert\bm{y}\rvert},\qquad\forall\,\lambda>0,\,\bm{y}\in\mathbb{R}^{3}\smallsetminus\\{\bm{0}\\}.$
(4.8)

According to (4.8), we obtain

$\begin{split}\Phi^{\lambda}_{\mathrm{off}}[\varphi]=&-24\pi\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{y}\>\lvert\varphi(\bm{y})\rvert^{2}\,\frac{\>e^{-\sqrt{\lambda}\mspace{1.5mu}\lvert\bm{y}\rvert}\\!\\!\\!}{\lvert\bm{y}\rvert}\;+\\\
&+\frac{8\sqrt{3}\mspace{1.5mu}\lambda}{\pi}\\!\int_{\mathbb{R}^{6}}\mspace{-7.5mu}d\bm{x}d\bm{y}\;\frac{\overline{\varphi(\bm{y})}\,[\varphi(\bm{y})-\varphi(\bm{x})]}{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}\,\mathrm{K}_{2}\Big{(}\sqrt{\tfrac{4\lambda}{3}}\,\textstyle{\sqrt{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}}\,\Big{)}.\end{split}$
(4.9)

It is now sufficient to notice that the symmetry in exchanging
$\bm{x}\longleftrightarrow\bm{y}$ allows us to write

$\displaystyle\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{y}\>\overline{\varphi(\bm{y})}\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{x}\;$
$\displaystyle\frac{\varphi(\bm{y})-\varphi(\bm{x})}{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}\,\mathrm{K}_{2}\Big{(}\sqrt{\tfrac{4\lambda}{3}}\,\textstyle{\sqrt{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}}\,\Big{)}=$
$\displaystyle=\frac{1}{2}\int_{\mathbb{R}^{6}}\mspace{-7.5mu}d\bm{x}d\bm{y}\;\frac{\lvert\varphi(\bm{y})-\varphi(\bm{x})\rvert^{2}\\!}{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}\,\mathrm{K}_{2}\Big{(}\sqrt{\tfrac{4\lambda}{3}}\,\textstyle{\sqrt{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}}\,\Big{)}.$

and the proposition is proved.

∎

Thanks to proposition 4.2, it is not hard to find lower and upper bounds for
$\Phi^{\lambda}$. To this end, it is convenient to introduce the Gagliardo
semi-norm of the Sobolev space $H^{1/2}(\mathbb{R}^{3})$, defined as

$[\mspace{1.5mu}u\mspace{1.5mu}]^{2}_{\frac{1}{2}}\\!:=\\!\int_{\mathbb{R}^{6}}\mspace{-7.5mu}d\bm{x}d\bm{y}\;\frac{\lvert
u(\bm{x})\\--u(\bm{y})\rvert^{2}\\!\\!}{\;\lvert\bm{x}-\bm{y}\rvert^{4}},\qquad
u\in H^{1/2}(\mathbb{R}^{3}),$ (4.10)

so that $\lVert u\rVert_{H^{1/2}(\mathbb{R}^{3})}^{2}\\!=\lVert
u\rVert^{2}\\!+[\mspace{1.5mu}u\mspace{1.5mu}]^{2}_{\frac{1}{2}}$. In terms of
the Fourier transform of $u$ we also have (see e.g., [4, proposition 3.4])

$[\mspace{1.5mu}u\mspace{1.5mu}]^{2}_{\frac{1}{2}}\\!=2\pi^{2}\mspace{-4.5mu}\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{k}\>\lvert\bm{k}\rvert\mspace{1.5mu}\lvert\hat{u}(\bm{k})\rvert^{2}.$
(4.11)

######  Proposition 4.3.

For any given $\varphi\in H^{1/2}(\mathbb{R}^{3})$, one has

$\displaystyle\Phi^{\lambda}[\varphi]\geq\Phi^{\lambda}_{\mathrm{diag}}[\varphi]+3\min\\{0,\gamma-2\\}[\mspace{1.5mu}\varphi\mspace{1.5mu}]^{2}_{\frac{1}{2}},$
(4.12a)
$\displaystyle\Phi^{\lambda}[\varphi]\leq\Phi^{\lambda}_{\mathrm{diag}}[\varphi]+\left(3\gamma+\tfrac{96\sqrt{3}}{\pi}\right)\\![\mspace{1.5mu}\varphi\mspace{1.5mu}]^{2}_{\frac{1}{2}}.$
(4.12b)

###### Proof.

The lower bound is obtained by neglecting the positive part of
$\Phi^{\lambda}_{\mathrm{off}}$

$\displaystyle\Phi^{\lambda}_{\mathrm{off}}[\varphi]+\Phi_{\mathrm{reg}}[\varphi]\geq
12\pi\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{y}\;\frac{\lvert\varphi(\bm{y})\rvert^{2}\\!}{\lvert\bm{y}\rvert}\left(\gamma-2\mspace{1.5mu}e^{-\sqrt{\lambda}\,\lvert\bm{y}\rvert}\right)$

and by considering the following inequalities

$\displaystyle\inf_{\bm{y}\in\,\mathbb{R}^{3}}\\!\left\\{\gamma-2\mspace{1.5mu}e^{-\sqrt{\lambda}\,\lvert\bm{y}\rvert}\right\\}\\!\geq\min\\{0,\gamma-2\\},$
(4.13)
$\displaystyle\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{x}\;\frac{\lvert\varphi(\bm{x})\rvert^{2}\\!}{\lvert\bm{x}\rvert}\leq\frac{1}{4\pi}[\mspace{1.5mu}\varphi\mspace{1.5mu}]^{2}_{\frac{1}{2}}.$
(4.14)

Notice that (4.14) is a consequence of the Hardy-Rellich inequality
(see555There is a typo in [15, equation (1.4)]: a power $2$ on the Euler Gamma
function in the numerator is missing. [15])

$\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{x}\;\frac{\lvert
u(\bm{x})\rvert^{2}\\!}{\lvert\bm{x}\rvert}\leq\frac{\pi}{2}\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{k}\>\lvert\bm{k}\rvert\mspace{1.5mu}\lvert\hat{u}(\bm{k})\rvert^{2},\qquad\forall\,u\in
H^{1/2}(\mathbb{R}^{3})$ (4.15)

compared with (4.11). In order to obtain the upper bound, we recall (4.10) to
get

$\displaystyle\int_{\mathbb{R}^{6}}\mspace{-7.5mu}d\bm{x}d\bm{y}\;\frac{\lvert\varphi(\bm{y})-\varphi(\bm{x})\rvert^{2}\\!}{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}$
$\displaystyle\,\mathrm{K}_{2}\Big{(}\sqrt{\tfrac{4\lambda}{3}}\,\textstyle{\sqrt{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}}\,\Big{)}\\!\leq$
$\displaystyle\leq[\mspace{1.5mu}\varphi\mspace{1.5mu}]_{\frac{1}{2}}^{2}\sup_{(\bm{x},\,\bm{y})\in\,\mathbb{R}^{6}}\,\frac{\lvert\bm{x}-\bm{y}\rvert^{4}}{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}\,\mathrm{K}_{2}\Big{(}\sqrt{\tfrac{4\lambda}{3}}\,\textstyle{\sqrt{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}}\,\Big{)}.$

We make use of (4.4d) and get rid of the dependence on the angles in
evaluating the $\sup$, since

$\begin{cases}x^{2}+y^{2}-2\,\bm{x}\\!\cdot\\!\bm{y}\leq
2\mspace{1.5mu}(x^{2}+y^{2}),\\\
x^{2}+y^{2}+\bm{x}\\!\cdot\\!\bm{y}\geq\frac{x^{2}+y^{2}}{2},\end{cases}\mspace{-18.0mu}\implies\frac{\lvert\bm{x}-\bm{y}\rvert^{2}}{x^{2}+y^{2}+\bm{x}\\!\cdot\\!\bm{y}}\leq
4.$

Hence,

$\displaystyle\int_{\mathbb{R}^{6}}\mspace{-7.5mu}d\bm{x}d\bm{y}\;\frac{\lvert\varphi(\bm{y})-\varphi(\bm{x})\rvert^{2}\\!\\-}{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}\,$
$\displaystyle\mathrm{K}_{2}\Big{(}\sqrt{\tfrac{4\lambda}{3}}\,\textstyle{\sqrt{y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}}}\Big{)}\leq$
$\displaystyle\leq\frac{3}{2\mspace{1.5mu}\lambda}\,[\mspace{1.5mu}\varphi\mspace{1.5mu}]_{\frac{1}{2}}^{2}\sup_{(\bm{x},\,\bm{y})\in\,\mathbb{R}^{6}}\frac{\lvert\bm{x}-\bm{y}\rvert^{4}}{\\!\left(y^{2}+x^{2}+\bm{x}\\!\cdot\\!\bm{y}\right)^{\\!2}\\!}=\frac{24}{\lambda}\,[\mspace{1.5mu}\varphi\mspace{1.5mu}]_{\frac{1}{2}}^{2}.$

We stress that this estimate is optimal, since the argument of the supremum in
$\mathbb{R}^{6}$ attains the previous upper bound along the hyperplane
$\bm{x}+\bm{y}=\bm{0}$.

So far, we have obtained for any $\varphi\in H^{1/2}(\mathbb{R}^{3})$

$\Phi^{\lambda}_{\mathrm{off}}[\varphi]\leq-24\pi\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{y}\;\frac{\lvert\varphi(\bm{y})\rvert^{2}\\!}{\lvert\bm{y}\rvert}\,e^{-\sqrt{\lambda}\,\lvert\bm{y}\rvert}+\frac{96\sqrt{3}}{\pi}\,[\mspace{1.5mu}\varphi\mspace{1.5mu}]_{\frac{1}{2}}^{2}.$
(4.16)

We complete the proof simply by neglecting the negative contribution.

∎

The major difficulties in the proof of the coercivity of $\Phi^{\lambda}$ in
momentum space obtained in [1] lie in the search of a lower bound. On the
other hand, in position space such estimate, provided in proposition 4.3,
turns out to be much easier. However, some accuracy is lost in this framework.
Indeed, adopting (4.12a) to obtain an estimate from below for
$\Phi^{\lambda}$, one gets

$\Phi^{\lambda}[\xi]\geq\Phi^{\lambda}_{\mathrm{diag}}[\xi]-3\max\\{0,2-\gamma\\}[\mspace{1.5mu}\xi\mspace{1.5mu}]^{2}_{\frac{1}{2}}+12\pi\,\mathop{\mathrm{essinf}}_{\bm{y}\,\in\,\mathbb{R}^{3}}\beta(\bm{y})\,\lVert\xi\rVert^{2}_{L^{2}(\mathbb{R}^{3})}$
(4.17)

or, equivalently

$\Phi^{\lambda}[\xi]\mspace{-1.5mu}\geq\\!\int_{\mathbb{R}^{3}}\mspace{-7.5mu}d\bm{p}\left[12\pi\sqrt{\tfrac{3}{4}p^{2}+\lambda}-6\pi^{2}\max\\{0,2\mspace{-1.5mu}-\mspace{-1.5mu}\gamma\\}p+12\pi\,\mathop{\mathrm{essinf}}_{\bm{y}\,\in\,\mathbb{R}^{3}}\beta(\bm{y})\right]\mspace{-1.5mu}\lvert\xi(\bm{p})\rvert^{2}.$
(4.18)

The function in square brackets attains its minimum at

$p_{\mathrm{min}}=\frac{\pi\max\\{0,2\mspace{-1.5mu}-\mspace{-1.5mu}\gamma\\}\sqrt{\lambda}}{\\!\sqrt{9-3\pi^{2}\max\\{0,2\mspace{-1.5mu}-\mspace{-1.5mu}\gamma\\}^{2}}\,},$

provided

$3-\pi^{2}\max\\{0,2-\gamma\\}^{2}>0\iff\gamma>2-\frac{\sqrt{3}}{\pi}=:\gamma_{c}^{\prime}.$
(4.19)

Indeed, plugging the value of $p_{\mathrm{min}}$ in the right hand side of
(4.18), one gets

$\Phi^{\lambda}[\xi]\mspace{-1.5mu}\geq\left(4\sqrt{3}\pi\mspace{1.5mu}\sqrt{\lambda}\sqrt{3-\pi^{2}\max\\{0,2\mspace{-1.5mu}-\mspace{-1.5mu}\gamma\\}^{2}}+12\pi\mathop{\mathrm{essinf}}_{\bm{y}\,\in\,\mathbb{R}^{3}}\beta(\bm{y})\right)\lVert\xi\rVert^{2}_{L^{2}(\mathbb{R}^{3})}$
(4.20)

that is positive for

$\lambda>\frac{3\min\\{0,\mathop{\mathrm{essinf}}\beta\\}^{2}}{3-\pi^{2}\max\\{0,2\mspace{-1.5mu}-\mspace{-1.5mu}\gamma\\}^{2}}.$
(4.21)

As mentioned above, $\gamma^{\prime}_{c}\approx 1.44867$ is not optimal, since
$\gamma^{\prime}_{c}\gneq\gamma_{c}\mspace{1.5mu}$.

## References

  * [1] G. Basti, C. Cacciapuoti, D. Finco and A. Teta, Three-body Hamiltonian with regularized zero-range interactions in dimension three. arXiv:2107.07188 [math-ph] Cited 15 Jul 2021. To appear in Annales H. Poincaré.
  * [2] M. Correggi, G. Dell’Antonio, D. Finco, A. Michelangeli and A. Teta, Stability for a system of N fermions plus a different particle with zero-range interactions. Rev. Math. Phys., 24, 1250017 (2012).
  * [3] G. F. Dell’Antonio, R. Figari and A. Teta, Hamiltonians for systems of N particles interacting through point interactions. Ann. Inst. H. Poincare, Phys. Theor., 60, 253–290 (1994).
  * [4] E. Di Nezza, G. Palatucci and E. Valdinoci, Hitchhiker’s guide to the fractional Sobolev spaces. Bull. Sci. math. 136, 521–573 (2012).
  * [5] J. Dimock and S. G. Rajeev, Multi-particle Schrödinger operators with point interactions in the plane. J. Phys. A, 37, 9157–9173 (2004).
  * [6] D. Ferretti and A. Teta, Regularized zero-range Hamiltonian for a Bose gas with an impurity. arXiv:2202.12765 [math-ph] Cited 25 Feb 2022.
  * [7] R. Figari and A. Teta, On the Hamiltonian for three bosons with point interactions. arXiv:2001.10462v1 [math-ph] Cited 28 Jan 2020. To appear in “Interplays between Mathematics and Physics through Stochastics and Infinite Dimensional Analysis: Sergio Albeverio’s contribution”, Springer (2022).
  * [8] D. Finco and A. Teta, Quadratic forms for the fermionic unitary gas model. Rep. Math. Phys., 69, 131–159 (2012).
  * [9] I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series, and Products \- 7th edition. Academic Press (2007).
  * [10] E. H. Lieb and M. Loss, Analysis \- 2nd edition. American Mathematical Society (2001).
  * [11] A. Michelangeli, Models of zero-range interaction for the bosonic trimer at unitarity. Rev. Math. Phys., 33, 2150010 (2021).
  * [12] R. A. Minlos and L. Faddeev, On the point interaction for a three-particle system in Quantum Mechanics. Soviet Phys. Dokl., 6, 1072–1074 (1962).
  * [13] R. A. Minlos and L. Faddeev, Comment on the problem of three particles with point interactions. Soviet Phys. Jetp., 14, 1315–1316 (1962).
  * [14] A. Posilicano, A Krein-like formula for singular perturbations of self-adjoint operators and applications. J. Funct. Anal. 183, 109–147 (2001).
  * [15] D. Yafaev, Sharp constants in the Hardy-Rellich inequalities. J. Funct. Anal., 168, 121–144 (1999).


# AuTO: A Framework for Automatic differentiation in Topology Optimization

Aaditya Chandrasekhar
Department of Mechanical Engineering
University of Wisconsin-Madison
Madison, WI
<EMAIL_ADDRESS>
&Saketh Sridhara
Department of Mechanical Engineering
University of Wisconsin-Madison
Madison, WI
<EMAIL_ADDRESS>
&Krishnan Suresh
Department of Mechanical Engineering
University of Wisconsin-Madison
Madison, WI
<EMAIL_ADDRESS>

###### Abstract

A critical step in topology optimization (TO) is finding sensitivities. Manual
derivation and implementation of the sensitivities can be quite laborious and
error-prone, especially for non-trivial objectives, constraints and material
models. An alternate approach is to utilize automatic differentiation (AD).
While AD has been around for decades, and has also been applied in TO, wider
adoption has largely been absent.

In this educational paper, we aim to reintroduce AD for TO, and make it easily
accessible through illustrative codes. In particular, we employ JAX, a high-
performance Python library for automatically computing sensitivities from a
user defined TO problem. The resulting framework, referred to here as AuTO, is
illustrated through several examples in compliance minimization, compliant
mechanism design and microstructural design.

##  1 Introduction

####

Fueled by improvements in manufacturing capabilities and computational
modeling, the field of topology optimization (TO) has witnessed tremendous
growth in recent years. To further accelerate the development of TO, we
consider here automating a critical step in TO, namely computing the
sensitivities, i.e., computing the derivatives of objectives, constraints,
material models, projections and filters, with respect to the design
variables, typically the elemental pseudo-densities, in the popular density-
based TO.

####

Conceptually, there are four different methods for computing sensitivities
[1]: (1) numerical, (2) symbolic, (3) manual, and (4) automatic
differentiation. Numerical, i.e., finite difference, based sensitivity
computation suffers from truncation and floating-point errors, and is
therefore not recommended. Symbolic differentiation using software packages
such as SymPy [2] or Mathematica [3] is a reasonable choice for simple
expressions. However, it is impractical when the quantity of interest involves
loops (such as when assembling stiffness matrices), and/or flow-control (if-
then-else). The default method today for computing sensitivities is manual.
While theoretically straightforward, the manual process is unfortunately
cumbersome and error prone; it is often the bottle-neck in the development of
new TO modules and exploratory studies. In this educational paper, we
therefore _illustrate and promote the use of automatic differentiation for
computing sensitivities in TO_.

####

Automatic differentiation (AD), is a collection of methods for efficiently and
accurately computing derivatives of numeric functions expressed as computer
programs [1]. AD has been around for decades [4] and has been exploited in a
wide range of problems ranging from molecular dynamics simulations [5] to the
design of photonic crystals [6]; see [7] for a critical review. AD in the
context of finite element analysis is reviewed in [8] and [9]. More recently,
AD was demonstrated for for shape optimization in [10] using the Firedrake
framework [11]. AD has also been exploited in TO of turbulent fluid flow
systems [12], [13].

####

Despite the pioneering research, AD is not widely used in TO. The objective of
this educational paper _is to accelerate the adoption of AD in TO by providing
standalone codes for popular TO problems._ In particular, we employ JAX [14],
a high-performance Python library for end-to-end AD. Its NumPy [15] like
syntax, low memory footprint and support of just-in-time (JIT) compilation for
accelerated code performance makes it an ideal candidate for the task. We
demonstrate the use of AD within the popular density-based TO framework [16],
by replicating existing educational TO codes for compliance minimization [17],
compliant mechanism design [18] and microstructural design [19]. Critical code
snippets are highlighted in this article; the complete codes are available at
https://github.com/UW-ERSL/AuTO

##  2  Compliance minimization

###  2.1 Problem Formulation

First we consider compliance minimization as modeled [17] in subject to a
volume constraint; this TO problem is very popular due to its self-adjoint
nature. In a mesh discretized form, the problem can be posed as:

$\displaystyle\underset{\bm{\rho}}{\text{minimize}}$ $\displaystyle
J=\bm{u}^{\mathsf{T}}\bm{K}(\bm{\rho})\bm{u}$ (1a) subject to
$\displaystyle\bm{K}(\bm{\rho})\bm{u}=\bm{f}$ (1b)
$\displaystyle\sum_{e}\rho_{e}v_{e}\leq V^{*}$ (1c)

where $\bm{u}$ is the displacement (in structural problems) or temperature (in
thermal problems), $\bm{K}$ is the stiffness matrix, $\bm{\rho}$ is the
pseudo-density design variables, $\bm{f}$ is the structural/thermal load and
$V^{*}$ is the volume constraint. To solve this problem, one must define the
material model (see below), rely on finite element analysis to solve Equation
1b, and use design update schemes such as MMA [20] or Optimality Criteria
[21].

####

A critical ingredient for the design update schemes is the sensitivity, i.e.,
derivative, of the objective and constraint with respect to the pseudo-density
variables. As mentioned earlier, this is typically carried out manually. For
the above self-adjoint problem, the sensitivity of the compliance, for the
solid isotropic material with penalization (SIMP) [21] material model, can be
easily derived:

$\frac{\partial
J}{\partial\rho_{e}}=-\bm{u}^{T}\frac{\partial\bm{K}}{\partial\rho_{e}}\bm{u}=-p\rho^{p-1}{u_{e}}^{T}K{u_{e}}$
(2)

_However, in this paper, we will rely on automatic differentiation (AD)
framework for sensitivity analysis._

###  2.2 AuTO Framework

####

The algorithm for solving the above compliance minimization is illustrated in
1. Code snippets that illustrate the use of AD are discussed below.

Algorithm 1 Compliance Minimization

1:procedure complianceMin(mesh, material, filter, BC, $V^{*}$)

2: $i=0$ $\triangleright$ Iteration index

3: $\rho=V^{*}$ $\triangleright$ Design variable initialization

4: $\Delta=1.0$ $\triangleright$ Design change

5: while $\Delta>\epsilon\;\text{and}\;i\leq\text{MaxIter}$ do

6: $i\leftarrow i+1$

7: $E\leftarrow\rho$ $\triangleright$ Material model

8: $\bm{K}\leftarrow E$ $\triangleright$ Compute stiffness matrix and assemble

9: $u\;\text{via}\;Ku=f$ $\triangleright$ Solve with imposed BC

10: $J\leftarrow(K,u)$ $\triangleright$ Objective

11: $\frac{\partial\mathcal{J}}{\partial\rho}\leftarrow AD(\rho\rightarrow J)$
$\triangleright$ Automatic differentiation of objective

12: $g\leftarrow(\bar{\rho},V^{*})$ $\triangleright$ Vol. Constraint

13: $\frac{\partial g}{\partial\rho}\leftarrow AD(\rho\rightarrow g)$
$\triangleright$ Automatic differentiation of constraint

14: $\phi^{i}\leftarrow(J,g,\frac{\partial J}{\partial\rho},\frac{\partial
g}{\partial\rho})$ $\triangleright$ MMA Solver [20]

15: $\Delta=(||\rho^{i}-\rho^{i-1}||)$

16: end while

17:end procedure

Steps 7-10 are captured through the following Python code, where the @jit
directive refers to the "just-in-compilation", i.e., the compiler translates
the Python functions to optimized machine code at run-time, approaching the
speeds of C or FORTRAN [22]. @jit def computeCompliance(rho): E =
MaterialModel(rho) K = assembleK(E) u = solveKuf(K) J = jnp.dot(u.T,
jnp.dot(K,u)) return J

SIMP [21] is a typical material model, and implemented as follows (the @jit
directive has been removed here to avoid repetition). def MaterialModel(rho):
E = Emin + (Emax-Emin)*rho**penal # SIMP return E The stiffness matrix is
assembled in a compact manner as follows. def assembleK(E): K =
jnp.zeros((ndof,ndof)) sK = D0.flatten()[np.newaxis]*E.T.flatten() K =
jax.ops.index_add(K, idx, sK) return K; where
$D0=\int\limits_{\Omega_{e}}[B]^{T}[C_{0}][B]d\Omega_{e}$ is the element base
stiffness matrix [23] with $E=1$ and prescribed $\nu$; idx reflects the global
numbering of the element nodes. The underlying linear system is solved using a
direct solver. def solveKuf(K): u_free =
jax.scipy.linalg.solve(K[free,:][:,free],force[free]) u = jnp.zeros((ndof)) u
= jax.ops.index_add(u, free,u_free.reshape(-1)) return u;

Finally, to compute the compliance and its sensitivity in step 11, we simply
request for the function and its derivative as follows. The JAX environment
automatically traces the chain of function calls, and ensures an end-to-end
automatic differentiation. J, gradJ = value_and_grad(computeCompliance)(rho)

The global volume constraint (in step 11) is defined as follows, @jit def
globalVolumeConstraint(rho): vc = jnp.mean(rho)/vf - 1. return vc;

As before, the value and its gradient (via AD) can be computed via

g, gradg = value_and_grad(globalVolumeConstraint)(rho)

As summarized in step 14 of the algorithm, the computed objective, objective
gradient, constraint, constraint gradient are then passed to standard
optimizers (MMA in our case) [20]. The reader is referred to the complete code
provided.

###  2.3 Illustrative Examples

####

We illustrate the above AD framework using two popular examples of compliance
minimization [16]: (a) minimizing structural compliance of a tip-loaded
cantilever (see Figure 1a) and (b) minimizing thermal compliance of a square
plate under a uniform heat load (see Figure 1b). The mesh was chosen to be 60
$\times$ 30 grid for the structural problem, and 60 $\times$ 60 grid for the
thermal problem. The target volume fraction in both problems is $V^{*}=0.5$.
The material properties are $E=1$, $\nu=0.3$, $k=1$, and MMA was used as the
design update scheme, with default parameters. The computed designs
illustrated in 1a and 1b matches those in the literature [16].

Figure 1: Compliance minimization examples: (a) Tip loaded cantilever and
optimized topology at $V^{*}=0.5$ (b) Heat conduction on a square plate and
optimized topology at $V^{*}=0.5$

The totla time taken for optimization using analytical derivatives and AD are
compared in Figure 2. We observe that AD is marginally more expensive, but we
will observe later that this is not always the case.

Figure 2: Computational cost for compliance minimization using AD and
analytical implementation.

###  2.4 Advantages of AD

####

While SIMP is a popular material model, other models have been proposed [24].
The advantage of AD is that one can easily replace SIMP with, for example,
RAMP [25], by simply changing the material model. def MaterialModel(rho): E =
Emax*rho/(1.+S*(1.-rho)) # RAMP return Y All downstream sensitivity
computations are handled automatically.

####

Often additional filters and projections are used in TO. For instance, they
can be used to remove checkerboard patterns [26], impose minimum length scale
[27], limit gray elements [28], etc. The filters apart from being complex in
their own right, they are often used in tandem. For instance, in [29], eight
such schemes were compounded to obtain shell-infill type structures. This
results in highly complicated sensitivity expressions that can be laborious to
derive. However, using an AD framework, the user simply needs to include the
desired projections in the pipeline and the sensitivity is taken care of. For
instance, we can introduce the following filter to reduce grayness in design,
just before computing the material model. def projectionFilter(rho):
if(projection[’isOn’]): nmr = np.tanh(c0*beta) + jnp.tanh(beta*(rho-c0)) dnmr
= np.tanh(c0*beta) + np.tanh(beta*(1-c0)) rho_tilde = nmr/dnmr return
rho_tilde else: return rho

####

Finally, manufacturing constraints [30], [31] are often imposed in TO; these
include limiting overhang of structures [32], connectivity [33], material
utilization [34], and length scale control [35]. Such constraints are easy to
impose within the AD framework. For example, the local volume constraint
proposed in [35] may be implemented as follows. def
maxLengthScaleConstraint(rho): v = jnp.matmul(L, (1.01-rho)**n); # L averaged
prior cons = 1 - jnp.power(jnp.sum(v**p),1./p)/vstar; return cons; As before,
one can calculate the value and gradient of the constraint via vc, gradvc =
value_and_grad(maxLengthScaleConstraint)(rho) The computed constraint and
gradient can then be passed on to MMA. To illustrate, for the tip cantilever
problem in Figure 1(a), with the additional max length scale radius of $r=30$,
and maximum void volume at $0.75\pi r^{2}$, the resulting topology is
illustrated in Figure 3.

Figure 3: Tip cantilever beam with length scale control.

##  3 Compliant Mechanism Design

We next illustrate the AuTO framework using compliant mechanisms (CMs) [36];
see [37] for a comprehensive review on TO for CMs.

###  3.1 Problem Formulation

Consider the displacement inverter considered in the 104-line educational
MATLAB code [16]. The objective is to maximize the output displacement
$u_{out}$ at the point of interest when a force $f_{in}$ is applied, as
illustrated in Figure 4. The spring constants are specified by the user to
control the behavior of the CM.

Figure 4: The displacement inverter compliant mechanism.

This TO problem can be written as:

$\displaystyle\underset{\bm{\rho}}{\text{maximize}}$
$\displaystyle\bm{u}_{out}$ (3a) subject to
$\displaystyle\bm{K}(\bm{\rho})\bm{u}=\bm{f_{in}}$ (3b)
$\displaystyle\sum_{e}\rho_{e}v_{e}\leq V^{*}$ (3c)

The standard implementation entails computing the elemental sensitivity
$\frac{\partial u_{out}}{\partial\rho_{e}}$ given by:

$\frac{\partial
u_{out}}{\partial\rho_{e}}=\bm{\lambda}^{T}\frac{\partial\bm{K}}{\partial\rho_{e}}\bm{u}=p\rho^{p-1}{\lambda_{e}}^{T}K{u_{e}}$
(4)

where $\bm{\lambda}$ is the solution to the adjoint load problem
$\bm{K\lambda}=-\bm{l}$. $\bm{l}$ is a vector with the value 1 at the degree
of freedom corresponding to the output point, and with zeros elsewhere.
Observe that, in the manual method, two sets of analysis (one to compute
$\bm{u}$, and the other to compute $\bm{\lambda}$) are required per iteration
for evaluating sensitivities.

###  3.2 AuTO Framework

The implementation in AuTO for CM design is similar to the compliance
minimization problem, with two minor changes: (a) the stiffness matrix
assembly includes the spring constants, and (b) the objective is the
displacement at the output node. The relevant code snippets are provided
below. def assembleKWithSprings(E): K = jnp.zeros((ndof,ndof)) sK =
D0.flatten()[np.newaxis]*E.T.flatten() K = jax.ops.index_add(K, idx, sK) #
springs at input and output nodes K = jax.ops.index_add(K,
jax.ops.index[nodeIn, nodeIn], kspringIn) K = jax.ops.index_add(K,
jax.ops.index[nodeOut, nodeOut], kspringOut) return K; def
CompliantMechanism(rho): E = MaterialModel(rho) K = assembleKWithSprings(E) u
= solveKuf(K) return u[bc[’nodeOut’]]

To compute the objective and its gradient, we rely on JAX as follows. J, gradJ
= value_and_grad(CompliantMechanism)(rho)

The design update using MMA is as per Section 2. Using the problem
specification in [18], the resulting topology for the inverter is illustrated
in Figure 5a; this is in agreement with the result in [18].

###  3.3 Advantages of AD

For the design of CMs using TO, a key advantage of AD stems from the following
observation [37] "no universally accepted objective formulation exists". For
example, consider two additional objectives:

  1. 1.

$\min:-\omega MSE+(1-\omega)SE$ [38]

  2. 2.

$\min:-MSE/SE$ [39]

where $MSE=\bm{v}^{T}\bm{K}\bm{u}$ is the mutual strain energy, which
describes the flexibility of the designed mechanism and
$SE=\bm{u}^{T}\bm{K}\bm{u}$ is the strain energy, $\bm{v}$ is the output
displacement when a unit dummy load applied at the degree of freedom
corresponding to the output point.

In the AuTO framework, one can easily explore various objectives as follows.
def CompliantMechanism(rho): E = MaterialModel(rho) K =
assembleKWithSprings(E) u = solveKuf(K) v = solve_dummy(K) MSE = jnp.dot(v.T,
jnp.dot(K,u)) SE = jnp.dot(u.T, jnp.dot(K,u)); J = -MSE/SE # or # w = 0.9 # J
= -w*MSE + (1 - w)*SE return J The topologies obtained with the two additional
objectives are illustrated in Figure 5b and Figure 5c.

Figure 5: Displacement inverter design using three formulations at
$V^{*}=0.35$

For the objective of maximizing output displacement, the computational costs
using analytical and AD methods are illustrated in Figure 6. Observe that the
the analytical method is more expensive since one must solve an adjoint
problem explicitly. On the other hand, JAX internally optimizes the code for
computing sensitivities via AD.

Figure 6: Computational cost of optimization using AD vs analytical
implementation.

##  4 Design of Materials

In this section, we replicate the educational article [19] for the design of
microstructures using AuTO. In particular, we consider (a) maximizing bulk
modulus (b) maximizing shear modulus, and (c) designing microstructures with
negative Poisson’s ratio.

###  4.1 Problem setup

The mathematical formulation is as follows [19]:

$\displaystyle\underset{\bm{\rho}}{\text{minimize}}$ $\displaystyle
c(E_{ijkl}^{H}(\bm{\rho}))$ (5a) subject to
$\displaystyle\bm{K}(\bm{\rho})\bm{U}^{A(kl)}=\bm{F}^{(kl)},k,l=1,2,\ldots,d$
(5b) $\displaystyle\sum_{e}\rho_{e}v_{e}\leq V^{*}$ (5c)

where the objective $c(E_{ijkl}^{H})$ represents the material property we
intend to minimize, $\bm{K}$ is the stiffness matrix, $\bm{U}^{A(kl)}$ and
$\bm{F}^{(kl)}$ are the displacement vector and the external force vector for
test case $(kl)$ respectively. The different test cases correspond to the unit
strain tests along different directions, where $d$ is the spatial dimension.

####

In 2D, maximization of bulk modulus corresponds to:

$c=-(E_{1111}+E_{1122}+E_{2211}+E_{2222})$ (6)

and maximization of shear modulus corresponds to:

$c=-E_{1212}$ (7)

Finally, for the design of materials with negative Poisson’s ratio, the
following was proposed [19]:

$c=-E_{1122}-\beta^{l}(E_{1111}+E_{2222})$ (8)

where $\beta\in(0,1)$ is a user-defined fixed parameter and $l$ is the design
iteration number. Observe that, in the manual method, computing the
sensitivity requires solving for the adjoint [19],

###  4.2 Implementation on AuTO

In AuTO, the bulk modulus objective, for example, can be captured as follows:

def MicrosructuralDesign(rho): E = MaterialModel(rho) K = assembleK(E) Kr, F =
computeSubMatrices(K) U = performFE(Kr, F) EMatrix = homogenizedMatrix(U, rho)
bulkModulus = -EMatrix[’0_0’]-EMatrix[’0_1’]-EMatrix[’1_1’]-EMatrix[’1_0’]
return bulkModulus Other objectives can be similarly captured. Figure 7
illustrates three different microstuctures for the three different objectives,
for a volume fraction of 0.25.

Figure 7: Maximization of bulk modulus, shear modulus and design of material
with negative Poisson’s ratio, with $v_{f}^{*}=0.25$.

##  5 Conclusion

####

In this paper, we demonstrated the simplicity and benefits of AD in TO.
Possible extensions include multi-physics [40, 41, 42] and non-linear problems
[43, 44]. In the current implementation, direct solvers were employed. For
large scale problems, sparse pre-conditioned iterative solvers [45], [46] will
be critical (but not fully supported by JAX). One of the advantages of the
manual approach (that is lost in AD) is that the expressions can provide key
insights to the problem.

## Replication of results

The Python code pertinent to this paper is available at https://github.com/UW-
ERSL/AuTO.

## Acknowledgments

The authors would like to thank the support of National Science Foundation
through grant CMMI 1561899.

## Compliance with ethical standards

The authors declare that they have no conflict of interest.

## References

  * [1] Atılım Günes Baydin, Barak A Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind.  Automatic differentiation in machine learning: a survey.  The Journal of Machine Learning Research, 18(1):5595–5637, 2017\.
  * [2] Aaron Meurer, Christopher P. Smith, Mateusz Paprocki, Ondřej Čertík, Sergey B. Kirpichev, Matthew Rocklin, AMiT Kumar, Sergiu Ivanov, Jason K. Moore, Sartaj Singh, Thilina Rathnayake, Sean Vig, Brian E. Granger, Richard P. Muller, Francesco Bonazzi, Harsh Gupta, Shivam Vats, Fredrik Johansson, Fabian Pedregosa, Matthew J. Curry, Andy R. Terrel, Štěpán Roučka, Ashutosh Saboo, Isuru Fernando, Sumith Kulal, Robert Cimrman, and Anthony Scopatz.  Sympy: symbolic computing in python.  PeerJ Computer Science, 3:e103, January 2017.
  * [3] Wolfram Research, Inc.  Mathematica, Version 12.2.  Champaign, IL, 2020.
  * [4] David E. Rumelhart, Geoffrey E. Hinton, and Ronald J. Williams.  Learning representations by back-propagating errors.  Nature, 323(6088):533–536, 1986.
  * [5] Samuel S Schoenholz and Ekin D Cubuk.  Jax md: End-to-end differentiable, hardware accelerated, molecular dynamics in pure python.  2019\.
  * [6] Momchil Minkov, Ian AD Williamson, Lucio C Andreani, Dario Gerace, Beicheng Lou, Alex Y Song, Tyler W Hughes, and Shanhui Fan.  Inverse design of photonic crystals through automatic differentiation.  ACS Photonics, 7(7):1729–1741, 2020.
  * [7] Louis B Rall.  Perspectives on automatic differentiation: past, present, and future?  In Automatic Differentiation: Applications, Theory, and Implementations, pages 1–14. Springer, 2006.
  * [8] I Ozaki, F Kimura, and M Berz.  Higher-order sensitivity analysis of finite element method by automatic differentiation.  Computational mechanics, 16(4):223–234, 1995.
  * [9] F Van Keulen, RT Haftka, and NH Kim.  Review of options for structural design sensitivity analysis. part 1: Linear systems.  Computer methods in applied mechanics and engineering, 194(30-33):3213–3243, 2005.
  * [10] Alberto Paganini and Florian Wechsung.  Fireshape: a shape optimization toolbox for firedrake.  Structural and Multidisciplinary Optimization, pages 1–17, 2021\.
  * [11] Peter Gangl, Kevin Sturm, Michael Neunteufel, and Joachim Schöberl.  Fully and semi-automated shape differentiation in ngsolve.  Structural and multidisciplinary optimization, pages 1–29, 2020\.
  * [12] Cetin B Dilgen, Sumer B Dilgen, David R Fuhrman, Ole Sigmund, and Boyan S Lazarov.  Topology optimization of turbulent flows.  Computer Methods in Applied Mechanics and Engineering, 331:363–393, 2018.
  * [13] Sumer B Dilgen, Cetin B Dilgen, David R Fuhrman, Ole Sigmund, and Boyan S Lazarov.  Density based topology optimization of turbulent flow heat transfer systems.  Structural and Multidisciplinary Optimization, 57(5):1905–1918, 2018.
  * [14] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang.  JAX: composable transformations of Python+NumPy programs, 2018\.
  * [15] Charles R. Harris, K. Jarrod Millman, St’efan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fern’andez del R’ıo, Mark Wiebe, Pearu Peterson, Pierre G’erard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant.  Array programming with NumPy.  Nature, 585(7825):357–362, September 2020.
  * [16] Martin Philip Bendsoe and Ole Sigmund.  Topology optimization: theory, methods, and applications.  Springer Science & Business Media, 2013.
  * [17] Ole Sigmund.  A 99 line topology optimization code written in matlab.  Structural and multidisciplinary optimization, 21(2):120–127, 2001\.
  * [18] M P Bendsoe and Ole. Sigmund.  Topology optimization: theory, methods, and applications.  Springer Berlin Heidelberg, 2 edition, 2003.
  * [19] Liang Xia and Piotr Breitkopf.  Design of materials using topology optimization and energy-based homogenization approach in matlab.  Structural and multidisciplinary optimization, 52(6):1229–1241, 2015.
  * [20] Krister Svanberg.  The method of moving asymptotes—a new method for structural optimization.  International journal for numerical methods in engineering, 24(2):359–373, 1987.
  * [21] Martin P Bendsøe and Ole Sigmund.  Optimization of structural topology, shape, and material, volume 414.  Springer, 1995.
  * [22] Siu Kwan Lam, Antoine Pitrou, and Stanley Seibert.  Numba: A llvm-based python jit compiler.  In Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC, pages 1–6, 2015.
  * [23] Klaus-Jürgen Bathe.  Finite element procedures.  Klaus-Jurgen Bathe, 2006.
  * [24] Grzegorz Dzierżanowski.  On the comparison of material interpolation schemes and optimal composite properties in plane shape optimization.  Structural and Multidisciplinary Optimization, 46(5):693–710, 2012\.
  * [25] Mathias Stolpe and Krister Svanberg.  An alternative interpolation scheme for minimum compliance topology optimization.  Structural and Multidisciplinary Optimization, 22(2):116–124, 2001\.
  * [26] O. Sigmund and J. Petersson.  Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, mesh-dependencies and local minima.  Structural Optimization, 16(1):68–75, 1998.
  * [27] James K Guest, Jean H Prévost, and Ted Belytschko.  Achieving minimum length scale in topology optimization using nodal design variables and projection functions.  International journal for numerical methods in engineering, 61(2):238–254, 2004.
  * [28] Jun Wu, Niels Aage, Ruediger Westermann, and Ole Sigmund.  Infill Optimization for Additive Manufacturing – Approaching Bone-like Porous Structures.  IEEE Transactions on Visualization and Computer Graphics, 24(2):1127–1140, 2018.
  * [29] Jun Wu, Anders Clausen, and Ole Sigmund.  Minimum compliance topology optimization of shell–infill composites for additive manufacturing.  Computer Methods in Applied Mechanics and Engineering, 326:358–375, Nov 2017.
  * [30] Sandro L. Vatanabe, Tiago N. Lippi, Cícero R.de Lima, Glaucio H. Paulino, and Emílio C.N. Silva.  Topology optimization with manufacturing constraints: A unified projection-based approach.  Advances in Engineering Software, 100:97–112, oct 2016.
  * [31] Jikai Liu, Andrew T Gaynor, Shikui Chen, Zhan Kang, Krishnan Suresh, Akihiro Takezawa, Lei Li, Junji Kato, Jinyuan Tang, Charlie CL Wang, et al.  Current and future trends in topology optimization for additive manufacturing.  Structural and Multidisciplinary Optimization, 57(6):2457–2483, 2018.
  * [32] Xiaoping Qian.  Undercut and overhang angle control in topology optimization: A density gradient based integral approach.  International Journal for Numerical Methods in Engineering, 111(3):247–272, Jul 2017.
  * [33] Quhao Li, Wenjiong Chen, Shutian Liu, and Liyong Tong.  Structural topology optimization considering connectivity constraint.  Structural and Multidisciplinary Optimization, 54(4):971–984, 2016\.
  * [34] Emily D. Sanders, Miguel A. Aguiló, and Glaucio H. Paulino.  Multi-material continuum topology optimization with arbitrary volume and mass constraints.  Computer Methods in Applied Mechanics and Engineering, 340:798–823, oct 2018.
  * [35] James K. Guest.  Imposing maximum length scale in topology optimization.  Structural and Multidisciplinary Optimization, 37(5):463–473, 2009\.
  * [36] Larry L Howell.  Compliant mechanisms.  In 21st century kinematics, pages 189–216. Springer, 2013.
  * [37] Benliang Zhu, Xianmin Zhang, Hongchuan Zhang, Junwen Liang, Haoyan Zang, Hai Li, and Rixin Wang.  Design of compliant mechanisms using continuum topology optimization: a review.  Mechanism and Machine Theory, 143:103622, 2020.
  * [38] Shinji Nishiwaki, Mary I Frecker, Seungjae Min, and Noboru Kikuchi.  Topology optimization of compliant mechanisms using the homogenization method.  International journal for numerical methods in engineering, 42(3):535–559, 1998.
  * [39] A Saxena and GK Ananthasuresh.  On an optimal property of compliant topologies.  Structural and multidisciplinary optimization, 19(1):36–49, 2000\.
  * [40] Joe Alexandersen and Casper Schousboe Andreasen.  A review of topology optimisation for fluid-based problems.  Fluids, 5(1):29, 2020.
  * [41] Johannes Semmler, Lukas Pflug, and Michael Stingl.  Material optimization in transverse electromagnetic scattering applications.  SIAM Journal on Scientific Computing, 40(1):B85–B109, 2018.
  * [42] Shiguang Deng and Krishnan Suresh.  Topology optimization under thermo-elastic buckling.  Structural and Multidisciplinary Optimization, 55(5):1759–1772, 2017.
  * [43] Fengwen Wang, Ole Sigmund, and Jakob Søndergaard Jensen.  Design of materials with prescribed nonlinear properties.  Journal of the Mechanics and Physics of Solids, 69:156–174, 2014\.
  * [44] Anders Clausen, Fengwen Wang, Jakob S Jensen, Ole Sigmund, and Jennifer A Lewis.  Topology optimized architectures with programmable poisson’s ratio over large deformations.  Adv. Mater, 27(37):5523–5527, 2015.
  * [45] Martin S Andersen, Joachim Dahl, and Lieven Vandenberghe.  Cvxopt: A python package for convex optimization.  abel. ee. ucla. edu/cvxopt, 2013.
  * [46] Praveen Yadav and Krishnan Suresh.  Assembly-Free Large-Scale Modal Analysis on the Graphics-Programmable Unit.  Journal of Computing and Information Science in Engineering, 13(1):011003, Jan 2013.