# Nodal and constant sign solutions for singular elliptic problems

Dumitru Motreanu  Université de Perpignan, Département de Mathématiques, 66860
Perpignan, France<EMAIL_ADDRESS>and  Abdelkrim Moussaoui  Applied
Mathematics Laboratory, Faculty of Exact Sciences
and Biology Department, Faculty of Natural and Life Sciences
A. Mira Bejaia University, Targa Ouzemour, 06000 Bejaia, Algeria
<EMAIL_ADDRESS>
###### Abstract.

We establish the existence of multiple solutions for singular quasilinear
elliptic problems with a precise sign information: two opposite constant sign
solutions and a nodal solution. The approach combines sub-supersolutions
method and Leray-Schauder topological degree involving perturbation argument.

###### Key words and phrases:

Singularity; perturbation; nodal solutions; sub-supersolutions; topological
degree.

###### 1991 Mathematics Subject Classification:

35J75, 35J62, 35J92.

##  1\. Introduction

Let $\Omega\subset\mathbb{R}^{N}$ $\left(N\geq 2\right)$ be a bounded domain
with $C^{1}$ boundary $\partial\Omega$. We consider the singular elliptic
problem

$(\mathrm{P})\qquad\left\\{\begin{array}[]{ll}-\Delta_{p}u=f(x,u)&\text{in}\;\Omega,\\\
u=0&\text{on}\;\partial\Omega,\end{array}\right.$

where $\Delta_{p}$ stands for the $p$-Laplacian operator on
$W_{0}^{1,p}(\Omega)$ with $1<p<\infty$ and the reaction term $f(x,u$) is
described by a Carathéodory function
$f:\Omega\times(\mathbb{R}\backslash\\{0\\})\rightarrow\mathbb{R}$ that can
exhibit singularity when $u$ approaches zero. This occurs under the following
hypothesis on the nonlinearity $f(x,s)$ reflecting the singular character of
problem $(\mathrm{P})$:

$\lim_{s\rightarrow 0}f(x,s)=+\infty\text{ \ uniformly for a.e. }x\in\Omega.$

For more information on the study of singular problems we refer to [5, 6, 7,
8, 10, 12, 11, 13, 14, 15, 16] and the references therein.

The presence of a singularity in $(\mathrm{P})$ represents a major obstacle to
overcome, especially since nonpositive solutions are addressed. Namely, beside
positive and negative solutions, we also discuss the existence of nodal
solutions for $(\mathrm{P})$ that are sign changing functions and inevitably
go through the singularity. As far as we know, such solutions for problems
involving singularities have been rarely investigated in the literature.
Actually, [16] is the only paper that has considered this issue. The existence
of nodal solutions is established for a class of semilinear singular system by
means of the trapping region formed by pairs of appropriately defined sign
changing sub-supersolutions. Exploiting spectral properties of Laplacian
operator as well as adequate truncation, the author establishes that the nodal
solution vanishes on negligible sets. This is an essential point enabling
nodal solutions investigation for singular problems.

In line with [16], by a solution of problem $(\mathrm{P})$ we mean any $u\in
W_{0}^{1,p}\left(\Omega\right)$ such that $f(x,u)\varphi\in L^{1}(\Omega)$ and

$\int_{\Omega}\left|\nabla u\right|^{p-2}\nabla u\nabla\varphi\
dx=\int_{\Omega}f(x,u)\varphi\ dx,$

for all $\varphi\in W_{0}^{1,p}(\Omega)$.

Our main goal is to establish multiplicity result for singular quasilinear
elliptic problem $(\mathrm{P})$, with a precise sign information. We provide
three solutions for $(\mathrm{P})$, two of which are of opposite constant-sign
solutions while the third one is nodal. Our approach combines sub-
supersolution method and topological degree theory. We first establish the
existence of opposite constant-sign solutions
$u_{+},u_{-}\in\mathcal{C}^{1,\tau}(\overline{\Omega})$. They are located in
positive and negative intervals formed by two opposite constant sign sub-
supersolutions pairs, constructed through a choice of suitable functions with
an adjustment of adequate constants. Then, focusing on the aforementioned
intervals, we show the existence of a minimal positive solution $u_{+}^{\ast}$
and a maximal negative solution $u_{-}^{\ast}$ for problem $(\mathrm{P})$. The
argument is essentially based on the Zorn’s Lemma as well as the
$S_{+}$-property of the negative $p$-Laplacian operator on
$W_{0}^{1,p}(\Omega)$.

A significant feature in the present work concerns the finding of a nodal
solution for singular problem $(\mathrm{P})$. It is a third solution for
$(\mathrm{P})$ achieved via the Leray-Schauder degree theory. The latter is
applied to a perturbed system $(\mathrm{P}_{\varepsilon}),$ depending on a
parameter $\varepsilon>0,$ whose study is relevant for problem $(\mathrm{P})$.
Precisely, we prove that the degree on a ball
$\mathcal{B}_{R_{\varepsilon}}(0)$, centered at $0$ of radius
$R_{\varepsilon}$, encompassing all potential solutions of
$(\mathrm{P}_{\varepsilon})$ is $0,$ while the degree on
$\mathcal{B}_{R_{\varepsilon}}(0)\backslash\mathcal{A}_{R_{\varepsilon}},$
where $\mathcal{A}_{R_{\varepsilon}}$ is a set containing only solutions lying
inside the interval $]u_{-}^{\ast},u_{+}^{\ast}[,$ is equal to $1$. By the
excision property of Leray-Schauder degree, this leads to the existence of a
nontrivial solution $u_{\varepsilon}^{\ast}$ for $(\mathrm{P}_{\varepsilon})$
in $\mathcal{A}_{R_{\varepsilon}}.$ Then, a solution
$u^{\ast}\in]u_{-}^{\ast},u_{+}^{\ast}[$ of $(\mathrm{P})$ is derived by
passing to the limit as $\varepsilon\rightarrow 0.$ The argument is based on a
priori estimates, dominated convergence theorem as well as $S_{+}$-property of
the negative $p$-Laplacian. Due to the extremality of constant-sign solutions
$u_{-}^{\ast}$ and $u_{+}^{\ast},$ then $u^{\ast}$ is a third solution of
$(\mathrm{P})$ that changes sign.

The rest of this article is organized as follows. Section 2 deals with
opposite constant sign solutions while Section 3 provides a nodal solution for
problem $(\mathrm{P})$.

##  2\. Constant sign solutions

Given $1<p<\infty$, the spaces $L^{p}(\Omega)$ and $W_{0}^{1,p}(\Omega)$ are
endowed with the usual norms $\|u\|_{p}=(\int_{\Omega}|u|^{p}\ dx)^{1/p}$ and
$\|u\|_{1,p}=(\int_{\Omega}|\nabla u|^{p}\ dx)^{1/p}$, respectively. For a
later use, set $p^{\prime}=\frac{p}{p-1}$. We will also use the spaces
$C(\overline{\Omega})$ and $C_{0}^{1,\beta}(\overline{\Omega})=\\{u\in
C^{1,\beta}(\overline{\Omega}):u=0\ \mbox{on $\partial\Omega$}\\}$ with
$\beta\in(0,1)$.

Let us denote by $\lambda_{1,p}$ the first eigenvalue of $-\Delta_{p}$ on
$W_{0}^{1,p}(\Omega)$ which is given by

(2.1) $\lambda_{1,p}=\inf_{u\in
W_{0}^{1,p}(\Omega)\setminus\\{0\\}}\frac{\int_{\Omega}|\nabla
u|^{p}\,\mathrm{d}x}{\int_{\Omega}|u|^{p}\,\mathrm{d}x}.$

Let $\phi_{1,p}$ be the positive normalized eigenfunction of $-\Delta_{p}$
corresponding to $\lambda_{1,p}$, that is,

$-\Delta_{p}\phi_{1,p}=\lambda_{1,p}\phi_{1,p}^{p-1}\ \text{in }\Omega,\ \
\phi_{1,p}=0\ \text{on }\partial\Omega,$

$\phi_{1,p}>0$ in $\Omega$, and $\|\phi_{1,p}\|_{p}=1$. Recall that there
exists a constant $c_{0}>0$ such that

(2.2) $\phi_{1,p}(x)\geq c_{0}d(x)\text{ for all }x\in\Omega,$

where $d(x)$ denotes the distance from $x\in\overline{\Omega}$ to the boundary
$\partial\Omega$.

For $u,v\in C^{1}(\overline{\Omega})$, the notation $u\ll v$ means

$\begin{array}[]{c}u(x)<v(x),\,\,\forall
x\in\Omega,\,\,\mbox{and}\,\,\,\frac{\partial v}{\partial\eta}<\frac{\partial
u}{\partial\eta}\,\,\,\mbox{on}\,\,\,\partial\Omega,\end{array}$

where $\eta$ is the outward unit normal to $\partial\Omega$.

######  Remark 1.

Under assumption $(\mathrm{H.1})$, any solution $u\in W_{0}^{1,p}(\Omega)$ of
$(\mathrm{P})$ satisfies $u(x)\neq 0$ for a.e. $x\in\Omega$. Otherwise, if $u$
vanishes on a set of positive measure, then $(\mathrm{P})$ and hypothesis
$(\mathrm{H}.1)$ are contradictory.

###  2.1. Opposite constant sign solutions

For a fixed $\delta>0$, set
$\Omega_{\delta}:=\left\\{x\in\Omega:d(x)<\delta\right\\}$. Let
$y,y_{\delta}\in\mathcal{C}^{1}\left(\overline{\Omega}\right)$ be the unique
solutions of the Dirichlet problems

(2.3) $-\Delta_{p}y=d(x)^{\alpha}\text{ in }\Omega,\text{ \ }y=0\text{ on
}\partial\Omega,$ (2.4)
$-\Delta_{p}y_{\delta}=\left\\{\begin{array}[]{ll}d(x)^{\alpha}&\text{in \
}\Omega\backslash\overline{\Omega}_{\delta},\\\ -1&\text{in \
}\Omega_{\delta},\end{array}\right.,\text{ \ }y_{\delta}=0\text{ on
}\partial\Omega,$

where $\alpha\in(-1,0)$. They satisfy

(2.5) $c^{-1}d(x)\leq y_{\delta}(x)\leq y(x)\leq cd(x)\text{ in }\Omega,$

with a constant $c>1$ (see [5, Lemma 3.1]).

For a constant $C>1$, we also pose

(2.6) $\overline{u}=Cy\text{ \ and \ }\underline{u}=C^{-1}y_{\delta}.$

Due to (2.5), this gives

(2.7) $\overline{u}\geq\underline{u}\ \text{ in }\ \overline{\Omega},$

so we can consider the ordered interval

$\left[\underline{u},\overline{u}\right]=\\{v\in
W_{0}^{1,p}(\Omega):\underline{u}(x)\leq v(x)\leq\overline{u}(x)\ \text{a.e.}\
x\in\Omega\\}.$

We further assume.

$(\mathrm{H}.2)$:

There exist constants $M>0$ and $-1<\alpha<0$ such that

$f(x,s)\leq M(1+s^{\alpha})\text{ \ and \ }f(x,-s)\geq-M(1+|s|^{\alpha}),$

for a.e. $x\in\Omega,$ all $s>0.$

$(\mathrm{H}.3)$:

There exist constants $\beta<0<m$ such that

$\alpha\geq\beta>-\min\\{1,p-1\\},$ $f(x,s)\geq ms^{\beta}\text{ \ and \
}f(x,-s)\leq-m|s|^{\beta},$

for a.e. $x\in\Omega,$ all $s>0.$

######  Theorem 1.

Assume that $(\mathrm{H.2})$ and $(\mathrm{H.3})$ hold. Then problem
$(\mathrm{P})$ admits a positive solution
$u_{+}\in\mathcal{C}^{1,\tau}(\overline{\Omega})$ and a negative solution
$u_{-}\in\mathcal{C}^{1,\tau}(\overline{\Omega})$ within the ordered intervals
$\left[\underline{u},\overline{u}\right]$ and
$\left[-\overline{u},-\underline{u}\right],$ respectively. Moreover, every
positive solution $u_{+}$ and negative solution $u_{-}$ of $(\mathrm{P})$
within $\left[0,\overline{u}\right]\cap\mathcal{C}^{1}(\overline{\Omega})$ and
$\left[-\overline{u},0\right]\cap\mathcal{C}^{1}(\overline{\Omega})$,
respectively, satisfy

(2.8) $\underline{u}(x)\leq u_{+}(x)\text{ \ and \
}u_{-}(x)\leq-\underline{u}(x),\text{ }\forall x\in\Omega.$

###### Proof.

By (2.5) one has

$\displaystyle C^{p-1}d(x)^{\alpha}\geq
M\left(\left\|d\right\|_{\infty}^{-\alpha}+(Cc^{-1})^{\alpha}\right)d(x)^{\alpha}$
$\displaystyle\geq
M\left(d(x)^{-\alpha}+(Cc^{-1})^{\alpha}\right)d(x)^{\alpha}$
$\displaystyle\geq M(1+(Cc^{-1}d(x))^{\alpha})\geq M(1+(Cy(x))^{\alpha})\text{
\ for all }x\in\overline{\Omega},$

provided $C>0$ is sufficiently large. Thus, from $(\mathrm{H.2})$, (2.3) and
(2.6), it follows that

(2.9) $\displaystyle-\Delta_{p}\overline{u}\geq M(1+\overline{u}^{\alpha})\geq
f(x,\overline{u})\text{ in }\Omega$

and

(2.10) $-\Delta_{p}(-\overline{u})\leq-M(1+\overline{u}^{\alpha})\leq
f(x,-\overline{u})\text{ in }\Omega.$

From (2.4), (2.5), (2.6), and $(\mathrm{H.3})$, we get

(2.11) $\displaystyle-\Delta_{p}\underline{u}\leq C^{-(p-1)}d(x)^{\alpha}\leq
C^{-\beta}m(cd(x))^{\beta}$ $\displaystyle\leq
C^{-\beta}my(x)^{\beta}=m\underline{u}(x)^{\beta}\leq
f(x,\underline{u}(x))\text{ in }\Omega$

and

(2.12) $-\Delta_{p}(-\underline{u})\geq-m\underline{u}(x)^{\beta}\geq
f(x,-\underline{u}(x))\text{ in }\Omega$

provided $C>1$ is sufficiently large.

In view of (2.7), (2.9)-(2.12), we deduce that $(\underline{u},\overline{u})$
and $(-\overline{u},-\underline{u})$ form sub-supersolutions pairs for
$(\mathrm{P})$. Then, the sub-supersolution method (see [8, Theorem 2]) leads
to the existence of a positive solution
$u_{+}\in\left[\underline{u},\overline{u}\right]\cap\mathcal{C}^{1,\tau}(\overline{\Omega})$
and a negative solution
$u_{-}\in\left[-\overline{u},-\underline{u}\right]\cap\mathcal{C}^{1,\tau}(\overline{\Omega})$,
with some $\tau\in(0,1)$.

Let
$u_{+}\in\left[0,\overline{u}\right]\cap\mathcal{C}^{1}(\overline{\Omega})$ be
a positive solution of $(\mathrm{P})$. By $(\mathrm{H.3})$, (2.5) and (2.6),
we get

$\displaystyle-\Delta_{p}u_{+}$ $\displaystyle=$ $\displaystyle f(x,u_{+})\geq
mu_{+}^{\beta}\geq m\overline{u}^{\beta}\geq m(Ccd(x))^{\beta}$
$\displaystyle\geq$ $\displaystyle
C^{-(p-1)}d(x)^{\alpha}\geq-\Delta_{p}\underline{u}(x)\text{ \ in }\Omega,$

provided $C>1$ is sufficiently large and $\delta>0$ sufficiently small in
(2.4). By the weak comparison principle, we infer that property (2.8) holds
true. Using a similar argument, we show the corresponding property for a
negative solution. ∎

###  2.2. Extremal solutions

######  Theorem 2.

Under assumptions $(\mathrm{H.2})$ and $(\mathrm{H.3})$, problem
$(\mathrm{P})$ admits a smallest positive solution
$u_{+}^{\ast}\in\mathcal{C}^{1}(\overline{\Omega})$ within
$[\underline{u},\overline{u}]$ and a biggest negative solution
$u_{-}^{\ast}\in\mathcal{C}^{1}(\overline{\Omega})$ within
$[-\overline{u},-\underline{u}]$. In addition, it holds

(2.13) $\begin{array}[]{c}\underline{u}\ll u_{+}^{\ast}\text{ \ and
}-\underline{u}\gg u_{-}^{\ast}\quad\text{in}\ \Omega.\end{array}$

###### Proof.

We only prove the existence of a smallest positive solution
$u_{+}^{\ast}\in\mathcal{C}^{1}(\overline{\Omega})$ within
$[\underline{u},\overline{u}]$. The proof for the biggest negative solution in
$[-\overline{u},-\underline{u}]$ can be carried out in a similar way.

Denote by $\mathcal{S}$ the set of all functions
$w\in[\underline{u},\overline{u}]$ that are solutions of $(\mathrm{P})$. From
Theorem 1 it is known that $\mathcal{S}$ is not empty. We claim that
$\mathcal{S}$ is downward directed. To this end, let
$u_{1},u_{2}\in\mathcal{S}$. Since $(\underline{u},u_{1})$ and
$(\underline{u},u_{2})$ form pairs of sub-supersolutions for $(\mathrm{P})$,
if we set $\tilde{u}=\min\\{u_{1},u_{2}\\}$, by virtue of [3, Theorem 3.20],
$(\underline{u},\tilde{u})$ is a pair of sub-supersolution, too. Then, owing
to [8, Theorem 2], there exists a solution of $(\mathrm{P})$ in
$[\underline{u},\tilde{u}]\cap\mathcal{C}^{1}(\overline{\Omega})$, which
proves the claim.

Now we show that Zorn’s lemma can be applied in $\mathcal{S}$. In this
respect, let us consider a chain $\mathcal{C}$ in $\mathcal{S}$ with respect
to the order $\geq$. There is a sequence $\\{u_{k}\\}_{k\geq
1}\subset\mathcal{C}$ such that $\inf C=\inf_{k\geq 1}u_{k}$ (see [4, page
336]) and we can assume that the sequence $\\{u_{k}\\}_{k\geq 1}$ is
decreasing. If $\hat{u}:=\inf\mathcal{C}$, one has $u_{k}\rightarrow\hat{u}$
a.e. in $\Omega$ and $\hat{u}\in[\underline{u},\overline{u}]$. Having that
each $u_{k}$ is a solution of $(\mathrm{P}),$ from $(\mathrm{H.2})$, (2.5) and
(2.6), we obtain

$\left\|\nabla u_{k}\right\|_{p}^{p}=\int_{\Omega}f(x,u_{k})u_{k}\ dx\leq
M\int_{\Omega}(u_{k}+u_{k}^{\alpha+1})\text{ }dx\leq
C_{0}(\left\|\overline{u}\right\|_{\infty}+\left\|\overline{u}\right\|_{\infty}^{1+\alpha}),$

with a constant $C_{0}>0$. Therefore $\\{u_{k}\\}$ is bounded in
$W_{0}^{1,p}(\Omega)$ and so up to a subsequence we have
$u_{k}\rightharpoonup\hat{u}\text{ in }W_{0}^{1,p}(\Omega)$. Using that

$\langle-\Delta_{p}u_{k},u_{k}-\hat{u}\rangle=\int_{\Omega}f(x,u_{k})(u_{k}-\hat{u})\
dx,$

$(\mathrm{H.2})$, (2.5), (2.6), $u_{k}\in\mathcal{S}$, we find

$\displaystyle|f(x,u_{k})(u_{k}-\hat{u})|$ $\displaystyle\leq$ $\displaystyle
M(1+u_{k}^{\alpha})|u_{k}-\hat{u}|\leq
M(1+(Cc)^{-\alpha}d(x)^{\alpha})|u_{k}-\hat{u}|$ $\displaystyle\leq$
$\displaystyle
2M(1+(Cc)^{-\alpha}d(x)^{\alpha})\left\|\overline{u}\right\|_{\infty}\text{ \
a.e. in }\Omega.$

In view of [10, Lemma, page 726], we have $d(x)^{\alpha}\in L^{1}(\Omega)$.
Then Lebesgue’s dominated convergence theorem implies that

$\underset{k\rightarrow\infty}{\lim}\langle-\Delta_{p}u_{k},u_{k}-\hat{u}\rangle=0,$

which allows us to invoke the $S_{+}$-property of $-\Delta_{p}$ on
$W_{0}^{1,p}(\Omega)$ that gives $u_{k}\longrightarrow\hat{u}$ in
$W_{0}^{1,p}(\Omega)$. We conclude that
$\hat{u}=\inf\mathcal{C}\in\mathcal{S}$. Zorn’s Lemma can be applied which
provides a minimal element $u_{+}^{\ast}$ of $\mathcal{S}$. Since
$\mathcal{S}\subset[\underline{u},\overline{u}]$, by $(\mathrm{H.2})$, (2.5)
and (2.6), we infer from the regularity result in [7, Lemma 3.1] that
$u_{+}^{\ast}\in\mathcal{C}^{1,\tau}(\overline{\Omega})$ for some
$\tau\in(0,1)$.

We claim that $u_{+}^{\ast}$ is the smallest solution of $(\mathrm{P})$ in
$\mathcal{S}$. Indeed, if $u\in\mathcal{S}$, there is
$\tilde{u}\in\mathcal{S}$ with $\tilde{u}\leq u_{+}^{\ast}$ and $\tilde{u}\leq
u$ because $\mathcal{S}$ is downward directed. Recalling that $u_{+}^{\ast}$
is a minimal element of $\mathcal{S}$, this results in
$u_{+}^{\ast}=\tilde{u}\leq u$. This proves the claim.

What remains to show is the first inequality in (2.13). The second one follows
similarly. By $(\mathrm{H.3})$, (2.5) and (2.6), for each compact set
$\mathrm{K}\subset\Omega,$ there is a constant $\sigma=\sigma(\mathrm{K})>0$
such that

$\displaystyle-\Delta_{p}u_{+}^{\ast}$ $\displaystyle=$ $\displaystyle
f(x,u_{+}^{\ast})\geq m(u_{+}^{\ast})^{\beta}\geq m\overline{u}^{\beta}\geq
m(Ccd(x))^{\beta}$ $\displaystyle>$ $\displaystyle
C^{-(p-1)}d(x)^{\alpha}+\sigma\geq-\Delta_{p}\underline{u}(x)+\sigma\text{ \
in }\mathrm{K}$

provided $C>1$ is sufficiently large and $\delta>0$ sufficiently small in
(2.4). By the strong comparison principle [1, Proposition $2.6$] , we infer
that property (2.13) holds true. A quite similar argument shows the
corresponding property for a biggest negative solution $u_{-}^{\ast}$. This
completes the proof. ∎

##  3\. Nodal solutions

Our main result is stated as follows.

######  Theorem 3.

Under assumptions $(\mathrm{H.1})$, $(\mathrm{H.2})$ and $(\mathrm{H.3})$,
problem $(\mathrm{P})$ possesses a solution $u_{\ast}\in
W_{0}^{1,p}\left(\Omega\right)$ such that $u_{\ast}\neq u_{+}$ and
$u_{\ast}\neq u_{-}$. Moreover, $u_{\ast}$ is a nodal (sign-changing)
solution.

We will make use of the Leray-Schauder degree theory. For any
$\varepsilon\in(0,1)$, we introduce the function

(3.1) $\gamma_{\varepsilon}(s):=\varepsilon\text{ }sgn(s),\quad\forall
s\in\mathbb{R},$

and corresponding to it the regularized problem

$(\mathrm{P}_{\varepsilon})\qquad\left\\{\begin{array}[]{ll}-\Delta_{p}u=f(x,u+\gamma_{\varepsilon}(u))&\text{in}\;\Omega,\\\
u=0&\text{on}\;\partial\Omega.\end{array}\right.$

Our goal is to prove that $(\mathrm{P}_{\varepsilon})$ admits a solution
within $[-\underline{u},\underline{u}]$ and then, passing to the limit as
$\varepsilon\rightarrow 0$, to get the existence of the desired solution
$u_{\ast}$ for problem $(\mathrm{P})$.

Theorem 1 provides two opposite constant-sign solutions $u_{\varepsilon,+}$
and $u_{\varepsilon,-}$ for problem $(\mathrm{P}_{\varepsilon})$ in
$\mathcal{C}^{1,\tau}(\overline{\Omega})$ with some $\tau\in(0,1)$ satisfying
(2.8). In particular, the set of solutions
$u_{\varepsilon}\in\mathcal{C}^{1,\tau}(\overline{\Omega})$ to
$(\mathrm{P}_{\varepsilon})$ is nonempty.

For any $R>0$, set

$\mathcal{A}_{R}=\left\\{u\in\mathcal{B}_{R}(0)\,:-\underline{u}\leq
u\leq\underline{u}\right\\}$

and

$\mathcal{B}_{R}(0)=\left\\{u\in\mathcal{C}^{1}(\overline{\Omega}):\,\left\|u\right\|_{\mathcal{C}^{1}(\overline{\Omega})}<R\right\\}.$

We are going to establish the solvability of problem
$(\mathrm{P}_{\varepsilon})$.

######  Theorem 4.

Under assumptions $(\mathrm{H.1})$, $(\mathrm{H.2})$ and $(\mathrm{H.3})$, the
regularized problem $(\mathrm{P}_{\varepsilon})$ possesses a solution
$u_{\varepsilon}\in]u_{-}^{\ast},u_{+}^{\ast}[$ for all $\varepsilon\in(0,1)$,
where $u_{-}^{\ast}$ and $u_{+}^{\ast}$ are the extremal constant-sign
solutions found in Theorem 2.

####  3.0.1. The topological degree on a ball
$\mathcal{B}_{R_{\varepsilon}}(0)$

We introduce the following homotopy class of problems

$(\mathrm{P}_{\varepsilon,t})\qquad-\Delta_{p}u=\mathrm{F}{{}_{\varepsilon,t}(x},u)\text{
in }\Omega,\text{ \ }u=0\text{ \ on }\partial\Omega,$

with

$\mathrm{F}{{}_{\varepsilon,t}(x},u)=tf(x,u+\gamma_{\varepsilon}(u))+(1-t)(\lambda_{1,p}(u^{+})^{p-1}+1),$

for $t\in[0,1]$ and $\varepsilon\in(0,1)$, where $s^{+}:=\max\\{0,s\\}.$ It is
known that the problem

$(\mathrm{P}_{\varepsilon,0})\qquad\left\\{\begin{array}[]{l}-\Delta_{p}u=\lambda_{1,p}(u^{+})^{p-1}+1\text{
in }\Omega\\\ u=0\text{ \ on }\partial\Omega\end{array}\right.$

does not admit solutions $u\in W_{0}^{1,p}(\Omega)$ (see [12, Proposition
9.64]).

######  Proposition 1.

Assume that $(\mathrm{H.2})$ and $(\mathrm{H.3})$ hold. Then each solution
$u_{\varepsilon}$ of $(\mathrm{P}_{\varepsilon,t})$ belongs to
$\mathcal{C}^{1}(\overline{\Omega})$ and satisfies

(3.2)
$\left\|u_{\varepsilon}\right\|_{\mathcal{C}^{1}(\overline{\Omega})}<R_{\varepsilon},$

for all $t\in(0,1]$ and $\varepsilon\in(0,1)$, provided the constant
$R_{\varepsilon}>0$ is sufficiently large.

###### Proof.

The nonlinear regularity theory implies that
$u_{\varepsilon}\in\mathcal{C}^{1}(\overline{\Omega})$. Due to the fact that
$(\mathrm{P}_{\varepsilon,0})$ has no solution, it suffices to prove that the
assertion holds for every $t\in[\delta,1]$ with any $\delta\in(0,1)$. Suppose
by contradiction that for every $n$ there exist $t_{n}\in[\delta,1]$ and a
solution $u_{\varepsilon,n}$ of $(\mathrm{P}_{\varepsilon,t_{n}})$ such that

$t_{n}\rightarrow t\in[\delta,1]\text{ \ and \
}\theta_{n}:=\|u_{\varepsilon,n}\|_{\mathcal{C}^{1}(\overline{\Omega})}\rightarrow\infty\text{
\ as }n\rightarrow\infty.$

Setting

$v_{\varepsilon,n}:=\frac{u_{\varepsilon,n}}{\theta_{n}}\in\mathcal{C}^{1}(\overline{\Omega}),$

problem $(\mathrm{P}_{\varepsilon,t_{n}})$ results in

(3.3)
$-\Delta_{p}v_{\varepsilon,n}=\frac{t_{n}}{\theta_{n}^{p-1}}f(x,u_{\varepsilon,n}+\gamma_{\varepsilon}(u_{\varepsilon,n}))+\frac{1-t_{n}}{\theta_{n}^{p-1}}(\lambda_{1,p}(u_{\varepsilon,n}^{+})^{p-1}+1).$

By $(\mathrm{H.2})$ and (3.1), we have the estimate

$\displaystyle\left|\frac{t_{n}}{\theta_{n}^{p-1}}f(x,u_{\varepsilon,n}+\gamma_{\varepsilon}(u_{\varepsilon,n}))+\frac{1-t_{n}}{\theta_{n}^{p-1}}(\lambda_{1,p}(u_{\varepsilon,n}^{+})^{p-1}+1)\right|$
$\displaystyle\leq
M(1+|u_{\varepsilon,n}+\gamma_{\varepsilon}(u_{\varepsilon,n})|^{\alpha})+\lambda_{1,p}({\normalsize
v}_{\varepsilon,n}^{+})^{p-1}+1$ $\displaystyle\leq
M(1+\varepsilon^{\alpha})+\lambda_{1,p}({\normalsize
v}_{\varepsilon,n}^{+})^{p-1}+1$ $\displaystyle\leq C(1+\|{\normalsize
v}_{\varepsilon,n}\|_{C^{1}(\overline{\Omega})}^{p-1})$

in $\Omega$, with some constant $C>0$ independent of $n$. Thanks to the
regularity up to the boundary (refer to [9]), we derive from (3.3) that
${\normalsize v}_{\varepsilon,n}$ is bounded in
$\mathcal{C}^{1,\tau}(\overline{\Omega})$ with certain $\tau\in(0,1)$. The
compactness of the embedding
$\mathcal{C}^{1,\tau}(\overline{\Omega})\subset\mathcal{C}^{1}(\overline{\Omega})$
implies along a subsequence that

(3.4) ${\normalsize v}_{\varepsilon,n}\rightarrow{\normalsize
v}_{\varepsilon}\text{ \ in }\mathcal{C}^{1}(\overline{\Omega})\text{.}$

Letting $n\rightarrow\infty$ in (3.3), on the basis of (3.4) we arrive at

$\left\\{\begin{array}[]{ll}-\Delta_{p}{\normalsize
v}_{\varepsilon}=(1-t)\lambda_{1,p}({\normalsize
v}_{\varepsilon}^{+})^{p-1}&\text{ in }\Omega,\\\ {\normalsize
v}_{\varepsilon}=0&\text{on }\partial\Omega.\end{array}\right.$

This ensures that ${\normalsize v}_{\varepsilon}={\normalsize
v}_{\varepsilon}^{+}$, which is nonzero because $\|{\normalsize
v}_{\varepsilon}\|_{\mathcal{C}^{1}(\overline{\Omega})}=1$. Acting with
${\normalsize v}_{\varepsilon}$ produces

$\int_{\Omega}|\nabla{\normalsize
v}_{\varepsilon}|^{p}\,\mathrm{d}x=(1-t)\lambda_{1,p}\int_{\Omega}{\normalsize
v}_{\varepsilon}^{p}\,\mathrm{d}x.$

We cannot have $t=1$ because ${\normalsize v}_{\varepsilon}\neq 0$. If $t<1$,
by (2.1) we reach the contradiction

$\lambda_{1,p}\int_{\Omega}{\normalsize
v}_{\varepsilon}^{p}\,\mathrm{d}x\leq(1-t)\lambda_{1,p}\int_{\Omega}{\normalsize
v}_{\varepsilon}^{p}\,\mathrm{d}x,$

which completes the proof. ∎

For every $\varepsilon\in(0,1)$, let us define the homotopy
$\mathrm{H}_{\varepsilon}:[0,1]\times\overline{\mathcal{B}}_{R_{\varepsilon}}(0)\rightarrow\mathcal{C}^{1}(\overline{\Omega})$
by

$\mathrm{H}_{\varepsilon}(t,u)=u-(-\Delta_{p})^{-1}\mathrm{F}{\
{}_{\varepsilon,t}(x},u)$

that is admissible for the Leray-Schauder topological degree by Proposition 1
and because the operator $(-\Delta_{p})^{-1}$ is compact.

The nonexistence of solutions to problem $(\mathrm{P}_{\varepsilon,0})$ yields

$\deg\left(\mathrm{H}_{\varepsilon}(0,\cdot),\mathcal{B}_{R_{\varepsilon}}(0),0\right)=0.$

Consequently, the homotopy invariance property implies

(3.5)
$\begin{array}[]{c}\deg\left(\mathrm{H}_{\varepsilon}(1,\cdot),\mathcal{B}_{R_{\varepsilon}}(0),0\right)=0\text{
for all }\varepsilon\in(0,1).\end{array}$

####  3.0.2. The degree on
$\mathcal{B}_{R_{\varepsilon}}(0)\backslash\overline{\mathcal{A}}_{R_{\varepsilon}}$.

Consider the problem

$(\mathrm{\tilde{P}}_{\varepsilon,t})\qquad-\Delta_{p}u=\mathrm{\tilde{F}}{{}_{\varepsilon,t}}({x,}u)\text{
in }\Omega,\text{ \ }u=0\text{ \ on }\partial\Omega,$

for all $t\in[0,1]$ and $\varepsilon\in(0,1)$, where

$\mathrm{\tilde{F}}{{}_{\varepsilon,t}}({x,}u)=tf(x,u+\gamma_{\varepsilon}(u))+(1-t)\text{
}sgn(u)\text{ }\frac{\phi_{1,p}}{\left\|\phi_{1,p}\right\|_{\infty}}.$

######  Proposition 2.

Assume that $(\mathrm{H.2})$ and $(\mathrm{H.3})$ are fulfilled. Then every
solution $u_{\varepsilon}$ of $(\mathrm{\tilde{P}}_{\varepsilon,t})$ belongs
to $\mathcal{C}^{1}(\overline{\Omega})$ and satisfies (3.2) for all
$t\in[0,1]$ and $\varepsilon\in(0,1)$, provided the constant
$R_{\varepsilon}>0$ is sufficiently large. In addition, every positive (resp.
negative) solution $u_{\varepsilon,+}$ (resp. $u_{\varepsilon,-}$) of
$(\mathrm{\tilde{P}}_{\varepsilon,t})$ verifies

(3.6) $\begin{array}[]{c}\underline{u}\ll u_{\varepsilon,+}\text{ \ (resp.
}-\underline{u}\gg u_{\varepsilon,-}\text{)}\quad\text{in}\
\Omega.\end{array}$

###### Proof.

By assumption $(\mathrm{H.2})$, one has

$\begin{array}[]{l}|\mathrm{\tilde{F}}{{}_{\varepsilon,t}}({x,}u)|\leq
M(1+\varepsilon^{\alpha})+1.\end{array}$

The Moser iteration can be implemented ensuring that $u_{\varepsilon}\in
L^{\infty}(\Omega)$. Next, the regularity theory up to the boundary (see [9])
entails the bound in (3.2).

We are going to show the first inequality in (3.6). The second one follows
similarly. Let $u_{\varepsilon,+}$ be a positive solution of
$(\mathrm{\tilde{P}}_{\varepsilon,t})$. Given $\delta>0$, on account of
$(\mathrm{H.3})$, (2.2), (3.1) and (3.2), we have

$\displaystyle C^{-(p-1)}d(x)^{\alpha}\leq C^{-(p-1)}\delta^{\alpha}$
$\displaystyle<$ $\displaystyle
tm(R_{\varepsilon}+1)^{\beta}+(1-t)\frac{c_{0}\delta}{\left\|\phi_{1,p}\right\|_{\infty}}$
$\displaystyle\leq$ $\displaystyle
tm(u_{\varepsilon,+}+\gamma_{\varepsilon}(u_{\varepsilon,+}))^{\beta}+(1-t)\frac{\phi_{1,p}}{\left\|\phi_{1,p}\right\|_{\infty}}$
$\displaystyle\leq$
$\displaystyle\mathrm{\tilde{F}}{{}_{\varepsilon,t}}({x,}u_{\varepsilon,+})\text{
\ in\ }\Omega\backslash\overline{\Omega}_{\delta}$

and

$-C^{-(p-1)}<0\leq\mathrm{\tilde{F}}{{}_{\varepsilon,t}}({x,}u_{\varepsilon,+})\text{
in }\Omega_{\delta},$

for all $t\in[0,1]$ and $\varepsilon\in(0,1)$, provided the constant $C>0$ is
large enough. Thus, for each compact set $\mathrm{K}\subset\subset\Omega$,
there is a constant $\sigma=\tau(\mathrm{K})>0$ such that

$\sigma+\left\\{\begin{array}[]{ll}d(x)^{\alpha}&\text{in \
}\Omega\setminus\overline{\Omega}_{\delta},\\\ -1&\text{in \
}\Omega_{\delta},\end{array}\right.<\mathrm{\tilde{F}}{{}_{\varepsilon,t}}({x,}u_{\varepsilon,+})\text{
in }\ \mathrm{K}.$

Hence, via (2.4), (2.6) and the strong comparison principle (see [1,
Proposition $2.6$]), we infer that $\underline{u}\ll u_{\varepsilon,+}$ in
$\overline{\Omega}$. Then Theorem 2 ends the proof. ∎

Let us define the homotopy
$\mathrm{\tilde{H}}_{\varepsilon}:[0,1]\times\overline{\mathcal{B}_{R_{\varepsilon}}(0)\backslash\mathcal{A}_{R_{\varepsilon}}}\rightarrow\mathcal{C}^{1}(\overline{\Omega})$
by

$\mathrm{\tilde{H}}_{\varepsilon}(t,u)=u-(-\Delta_{p})^{-1}\mathrm{\tilde{F}}{{}_{\varepsilon,t}(x},u)$

for all $t\in[0,1]$ and $\varepsilon\in(0,1)$. The homotopy
$\mathrm{\tilde{H}}_{\varepsilon}$ is admissible for the Leray-Schauder
topological degree by Proposition 2 and because the operator
$(-\Delta_{p})^{-1}$ is compact.

The Minty-Browder Theorem (see, e.g., [2, Theorem 5.16]) and the regularity
result in [9] guarantee that the problem

$-\Delta_{p}u=\frac{\phi_{1,p}}{\left\|\phi_{1,p}\right\|_{\infty}}\text{ in
}\Omega,\text{ }u=0\text{ on }\partial\Omega,$

admits a unique solution $u\in\mathcal{C}^{1}(\overline{\Omega})$. Then the
homotopy invariance property of the Leray-Schauder degree shows

$\begin{array}[]{ll}\deg(\mathrm{\tilde{H}}_{\varepsilon}(1,\cdot),\mathcal{B}_{R_{\varepsilon}}(0)\backslash\overline{\mathcal{A}}_{R_{\varepsilon}},0)&=\deg(\mathrm{\tilde{H}}_{\varepsilon}(0,\cdot),\mathcal{B}_{R_{\varepsilon}}(0)\backslash\overline{\mathcal{A}}_{R_{\varepsilon}},0)\\\
&=1.\end{array}$

Since

$\mathrm{H}_{\varepsilon}(1,\cdot)=\mathrm{\tilde{H}}_{\varepsilon}(1,\cdot)\quad\text{in}\,\,\,\mathcal{B}_{R_{\varepsilon}}(0)\backslash\mathcal{A}_{R_{\varepsilon}},$

we deduce that

(3.7)
$\begin{array}[]{c}\deg(\mathrm{H}_{\varepsilon}(1,\cdot),\mathcal{B}_{R_{\varepsilon}}(0)\backslash\overline{\mathcal{A}}_{R_{\varepsilon}},0)=1.\end{array}$

####  3.0.3. Proof of Theorem 4.

We may assume that

$\mathrm{H}_{\varepsilon}(1,u)\not=0\quad\text{for all
}u\in\partial\mathcal{A}_{R_{\varepsilon}}\text{ and }\varepsilon\in(0,1).$

Otherwise, $u\in\partial\mathcal{A}_{R_{\varepsilon}}$ would be a solution of
$(\mathrm{P}_{\varepsilon})$ within
$\left[-\underline{u},\underline{u}\right]$ and thus Theorem 4 is proven.

By virtue of the domain additivity property of Leray-Schauder degree it
follows that

$\displaystyle\deg(\mathrm{H}_{\varepsilon}(1,\cdot),\mathcal{B}_{R_{\varepsilon}}(0)\backslash\overline{\mathcal{A}}_{R_{\varepsilon}},0)+\deg(\mathrm{H}_{\varepsilon}(1,\cdot),\mathcal{A}_{R_{\varepsilon}},0)$
$\displaystyle=$
$\displaystyle\deg(\mathrm{H}_{\varepsilon}(1,\cdot),\mathcal{B}_{R_{\varepsilon}}(0),0).$

Hence, by (3.5) and (3.7), we deduce

$\deg(\mathrm{H}_{\varepsilon}(1,\cdot),\mathcal{A}_{R_{\varepsilon}},0)=-1,$

obtaining that problem $(\mathrm{P}_{\varepsilon})$ has a solution
$u_{\varepsilon}\in\mathcal{A}_{R_{\varepsilon}},$ for all
$\varepsilon\in(0,1)$. Furthermore, the nonlinear regularity theory up to the
boundary (see [9]) guarantees that
$u_{\varepsilon}\in\mathcal{C}^{1,\tau}(\overline{\Omega})$ with some
$\tau\in(0,1)$.

###  3.1. Proof of Theorem 3.

Set $\varepsilon=\frac{1}{n}$ in $(\mathrm{P}_{\varepsilon})$ with any integer
$n\geq 1$. Theorem 4 provides a solution $u_{n}:=u_{\frac{1}{n}}$ of problem
$(\mathrm{P}_{\frac{1}{n}})$ such that

(3.8) $u_{n}\in[-\underline{u},\underline{u}],\quad\forall n.$

There it holds

(3.9) $\int_{\Omega}|\nabla u_{n}|^{p-2}\nabla
u_{n}\text{\thinspace}\nabla\varphi\text{
}\mathrm{d}x=\int_{\Omega}f(x,u_{n}+\gamma_{\frac{1}{n}})\varphi\ \mathrm{d}x$

for all $\varphi\in W_{0}^{1,p}(\Omega)$. Acting with $\varphi=u_{n}$, we find
from $(\mathrm{H.2})$ and (3.8) the estimate

$\displaystyle\int_{\Omega}|\nabla u_{n}|^{p}\text{
}\mathrm{d}x=\int_{\Omega}f(x,u_{n}+\gamma_{\frac{1}{n}}(u_{n}))u_{n}\
\mathrm{d}x$ $\displaystyle\leq$ $\displaystyle
M\int_{\Omega}(1+|u_{n}+\gamma_{\frac{1}{n}}(u_{n})|^{\alpha})|u_{n}|\
\mathrm{d}x$ $\displaystyle\leq$ $\displaystyle
M\int_{\Omega}(|u_{n}|+|u_{n}|^{\alpha+1})\ \mathrm{d}x$ $\displaystyle\leq$
$\displaystyle
M|\Omega|(\left\|\overline{u}\right\|_{\infty}+\left\|\overline{u}\right\|_{\infty}^{\alpha+1})<\infty.$

Hence, along a subsequence it holds $u_{n}\rightharpoonup u_{\ast}$ in
$W_{0}^{1,p}(\Omega)$ with some $u_{\ast}\in W_{0}^{1,p}(\Omega)$. By (3.8),
we note that

(3.10) $-\underline{u}\leq u_{\ast}\leq\underline{u}\ \text{ in }\Omega.$

Inserting $\varphi=u_{n}-u_{\ast}$ yields

$\int_{\Omega}|\nabla u_{n}|^{p-2}\nabla
u_{n}\text{\thinspace}\nabla(u_{n}-u_{\ast})\text{
}\mathrm{d}x=\int_{\Omega}f(x,u_{n}+\gamma_{\frac{1}{n}}(u_{n}))(u_{n}-u_{\ast})\
\mathrm{d}x.$

Then, from $(\mathrm{H.2})$ and (3.8), we see that

$\underset{n\rightarrow\infty}{\lim}\int_{\Omega}|\nabla u_{n}|^{p-2}\nabla
u_{n}\text{\thinspace}\nabla(u_{n}-u_{\ast})\text{ }\mathrm{d}x=0.$

At this point, the $S_{+}$-property of $-\Delta_{p}$ on $W_{0}^{1,p}(\Omega)$
(see, e.g., [12, Proposition 3.5]) implies

$u_{n}\rightarrow u_{\ast}\text{ in }W_{0}^{1,p}(\Omega).$

By virtue of $(\mathrm{H.2})$, Lebesgue’s dominated convergence theorem
entails

$\lim_{n\rightarrow\infty}\int_{\Omega}f(x,u_{n}+\gamma_{\frac{1}{n}}(u_{n}))\varphi\
\mathrm{d}x=\int_{\Omega}f(x,u_{\ast})\varphi\ \mathrm{d}x,$

for all $\varphi\in W_{0}^{1,p}(\Omega)$. Consequently, we may pass to the
limit in (3.9) finding that $u_{\ast}$ is a solution of problem
$(\mathrm{P})$. Remark 1 renders that the solution $u_{\ast}$ is nontrivial.
Moreover, due to (3.6), (3.10) and Theorems 1 and 2, we infer that the
solution $u_{\ast}$ of problem $(\mathrm{P})$ cannot be of constant sign
(otherwise the opposite constant-sign solutions $u_{-}^{\ast}$ and
$u_{+}^{\ast}$ were not extremal as known from Theorem 2). Therefore the
solution $u_{\ast}$ is nodal. This completes the proof.

## References

  * [1] D. Arcoya & D. Ruiz, _The Ambrosetti-Prodi problem for the p-Laplace operator_ , Comm. Partial Diff. Eqts. 31 (2006), 849-865.
  * [2] H. Brézis, _Functional Analysis, Sobolev Spaces and Partial Differential Equations_ , Springer, New York, 2010.
  * [3] S. Carl, V. K. Le & D. Motreanu, Nonsmooth Variational Problems and Their Inequalities. Comparison Principles and Applications, Springer, New York, 2007.
  * [4] N. Dunford & J. T. Schwartz, _Linear Operators. I. General Theory_ , Interscience Publishers, Inc., New York, (1958).
  * [5] H. Dellouche & A. Moussaoui, _Singular quasilinear elliptic systems with gradient dependence_ , Positivity 26 (2022), 10. doi:10.1007/s11117-022-00868-3.
  * [6] H. Didi & A. Moussaoui, _Multiple positive solutions for a class of quasilinear singular elliptic systems_ , Rend. Circ. Mat. Palermo, II. Ser. 69 (2020), 977-994.
  * [7] D. D. Hai, _On a class of singular p-Laplacian boundary value problems_ , J. Math. Anal. Appl. 383 (2011), 619-626.
  * [8] B. Khodja & A. Moussaoui, _Positive solutions for infinite semipositone/positone quasilinear elliptic systems with singular and superlinear terms_ , Diff. Eqts. App. 8(4) (2016), 535-546.
  * [9] G. M. Lieberman, _Boundary regularity for solutions of degenerate elliptic equations_ , Nonlinear Anal. 12 (1988), 1203-1219.
  * [10] A. C. Lazer & P. J. Mckenna, _On a singular nonlinear elliptic boundary-value problem_ , Proc. American Math. Soc. 3 (111), 1991.
  * [11] S.A. Marano, G. Marino, A. Moussaoui, _Singular quasilinear elliptic systems in_ $\mathbb{R}^{N}$, Ann. Mat. Pura Appl. 198(4) (2019), 1581-1594.
  * [12] D. Motreanu, V.V. Motreanu & N. Papageorgiou, _Topological and Variational methods with applications to Nonlinear Boundary Value Problems_ , Springer, New York, 2014.
  * [13] D. Motreanu, A. Moussaoui, _An existence result for a class of quasilinear singular competitive elliptic systems_ , Appl. Math. Lett.__ 38 (2014), 33-37.
  * [14] D. Motreanu, A. Moussaoui, A quasilinear singular elliptic system without cooperative structure, Acta Math. Sci.__ 34 (B) (2014), 905-916.
  * [15] D. Motreanu, A. Moussaoui, _Existence and boundedness of solutions for a singular cooperative quasilinear elliptic system_ , Complex Var. Elliptic Equ. 59 (2014), 285-296.
  * [16] A. Moussaoui, _Nodal solutions for singular semilinear elliptic systems_ , FILOMAT 15 (37) (2023), 5107-5117.


# Correlation energy and quantum correlations in a solvable model

Javier Faba  Center for Computational Simulation, Universidad Politécnica de
Madrid, Campus Montegancedo, 28660 Boadilla del Monte, Madrid, Spain
Departamento de Física Teórica and CIAFF, Universidad Autónoma de Madrid,
E-28049 Madrid, Spain  Vicente Martín  Center for Computational Simulation,
Universidad Politécnica de Madrid, Campus Montegancedo, 28660 Boadilla del
Monte, Madrid, Spain  Luis Robledo  Departamento de Física Teórica and CIAFF,
Universidad Autónoma de Madrid, E-28049 Madrid, Spain  Center for
Computational Simulation, Universidad Politécnica de Madrid, Campus
Montegancedo, 28660 Boadilla del Monte, Madrid, Spain

###### Abstract

Typically in many-body systems the correlation energy, which is defined as the
difference between the exact ground state energy and the mean-field solution,
has been a measure of the system’s total correlations. However, under the
quantum information context, it is possible to define some quantities in terms
of the system’s constituents that measure the classical and quantum
correlations, such as the entanglement entropy, mutual information, quantum
discord, one-body entropy, etc. In this work, we apply concepts of quantum
information in fermionic systems in order to study traditional correlation
measures (the relative correlation energy) from a novel approach. Concretely,
we analyze the two and three level Lipkin models, which are exactly solvable
(but non trivial) models very used in the context of the many-body problem.

###### pacs:

Valid PACS appear here

##  I Introduction

The atomic nucleus is a mesoscopic system made of protons and neutrons with
strong interactions among its constituents. Due to the complexity of the
nuclear interaction and the large number of particles involved, the dynamic
governing the nucleus is very rich, giving rise to a huge amount of different
situations involving single particle and/or collective excitations Ring and
Schuck (2004); Bender  _et al._ (2003). As a consequence of the underlying
mean field which implies the existence of well defined orbits, low energy
nuclear properties can dramatically change by changing a few units of the
nucleus’ proton and neutron numbers as a consequence of the filling of
different orbits. At low excitation energies, the so-called collective
excitations show more regular patterns than the single particle ones. The
reason it that they are associated to more macroscopic-like degrees of freedom
as the shape of the nucleus and are intimately connected with the mechanism of
spontaneous symmetry breaking and symmetry restoration Robledo  _et al._
(2019); Sheikh  _et al._ (2019). In finite systems this mechanism can be
viewed as an artifact of the underlying mean field description to capture
correlations in a simple way. Nevertheless the breaking for symmetries at the
mean field level is intimately connected to properties of the exact wave
functions of the system. The subsequent symmetry restoration of the symmetry-
broken mean-field wave functions gives rise to collective bands (being
rotational bands the most prominent example) that represent a prominent part
of the nuclear spectrum with very specific and universal properties like the
$I(I+1)$ energy rule of rotational bands Bender  _et al._ (2003); Robledo  _et
al._ (2019); Sheikh  _et al._ (2019). To improve upon the mean field plus
symmetry restoration paradigm, one usually add an additional layer where
fluctuations on the collective degrees of freedom are explicitly treated. This
is usually done in the framework of the Generator Coordinate Method Bender
_et al._ (2003); Robledo  _et al._ (2019). A question that arises very often
is how to quantify the balance between the correlations associated to symmetry
restoration and quantum fluctuations. The answer to this question might help
to devise new approaches to solve the nuclear many body problem. On the other
hand, the connection between the exact solution of the problem and the
approximate mean field plus symmetry restoration plus fluctuations approach is
not straightforward and there has been quite a lot of work to extract from the
exact shell model solution Caurier  _et al._ (2005) the underlying symmetry
breaking mean field. Therefore, it is also interesting to find a quantity to
be computed with the exact solution of the problem that is able to pin-point
the quantum phase transitions observed in the mean field description of the
nucleus. This is an approach also pursued in other fields like quantum
chemistry Legeza and Sólyom (2006); Szalay  _et al._ (2015), superconductors
in condensed matter Zeng  _et al._ (2014), atomic physics Tichy  _et al._
(2011) and even nuclear physics Kanada-En’yo (2015); Kruppa  _et al._ (2021).
With these two goals in mind we analyze in this paper some quantum information
related quantities as the overall entropy expressed in the basis of the
natural states, the quantum discord and the correlation energy. We will carry
out our study in the realm of a simple, albeit rich, exactly solvable nuclear
physics problem: the Lipkin model with two Lipkin  _et al._ (1965) and three
Li  _et al._ (1970); Holzwarth and Yukawa (1974a) active orbits. Both models
show quantum phase transitions as a function of the interaction parameter
strength that mimic the spontaneous symmetry breaking mechanism discussed
above.

##  II Theoretical background

In this section we will introduce briefly some concepts that we will use in
the next sections. When dealing with correlations in a many-body system, one
has to clarify two fundamental issues: what are we defining as subsystem, and
how to quantify the correlations among them? If our Hilbert space is defined
as a tensor product of Hilbert spaces, then the notion of subsystem arises
naturally. For example, the Hilbert space of a system formed by $N$ qubits is
simply the tensor product of each qubit’s Hilbert space. However, if we are
dealing with indistinguishable particles (fermions in our case) in the context
of second quantization, we cannot define the Hilbert space as a tensor product
of each particle’s Hilbert space because of the (anti)symmetry of the
wavefunction. A lot of effort has been made to disentangle the correlations
associated to the symmetrization principle or super-selection rules from the
ones coming from the dynamic of the system Legeza and Sólyom (2003); Bañuls
_et al._ (2007); Ding  _et al._ (2021) and quantities like the fermionic
partial trace between modes Friis  _et al._ (2013) or the von Neumann entropy
of the one-body density matrix Kanada-En’yo (2015) have been defined. Those
quantities have been thoroughly used in the literature Legeza and Sólyom
(2003, 2006); Szalay  _et al._ (2015); Gigena and Rossignoli (2015); Kruppa
_et al._ (2021); Robin  _et al._ (2021). In this work, we will discuss
quantities that make use of both concepts.

There are many possibilities in order to characterize and quantify
correlations in a quantum system. Typically, if our Hilbert space can be
written as $\mathcal{H}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}$111As discussed
above, when dealing with indistinguishable particles in the second
quantization formalism we don’t have a tensor product structure. However, if
we define the subsystems as the single particle states (also called orbitals
or modes in this work), we can treat the system as a tensor product if we take
into account some subtleties which arise from the fermionic anticommutation
rules Friis  _et al._ (2013). we can measure the entanglement between the $A$
and $B$ subsystems for a given pure state $\ket{\psi}\in\mathcal{H}$ through
the von Neumann entropy of the reduced states, namely:

$S(\rho^{(A)})=-\Tr(\rho^{(A)}\ln\rho^{(A)})$ (1)

where $\rho^{(A)}=\Tr_{B}(\ket{\psi}\bra{\psi})$ Nielsen and Chuang (2011).
However, if we are dealing with mixed states, this method is no longer valid
as an entanglement measure. Furthermore, entanglement is not only the only
type of correlation present in a quantum system: it can also have classical
correlations, and quantum correlations beyond entanglement.

The quantum discord Ollivier and Zurek (2001) is a measurement-based quantity
of the total quantum correlations (including entanglement and beyond) between
two subsystems. It is defined as

$\delta(A,B)=I(A,B)-J(A,B).$

where $I(A,B)=S(\rho^{(A)})+S(\rho^{(B)})-S(\rho^{(A,B)})$ is the mutual
information, and $J(A,B)$ is defined as

$J(A,B)=\max_{\\{\Pi_{k}^{(B)}\\}}S(\rho^{(A)})-S(\rho^{(A,B)}|\\{\Pi_{k}^{(B)}\\}).$
(2)

While $I(A,B)$ is a measure of all kind of correlations, $J(A,B)$ quantifies
only the classical part. The measurement-based conditional entropy in Eq. (2)
is defined as

$S(\rho^{(A,B)}|\\{\Pi_{k}^{(B)}\\})=\sum_{k}p_{k}S(\rho_{k}^{(A,B)})$

where $\rho_{k}^{(A,B)}=\frac{1}{p_{k}}\Pi_{k}^{(B)}\rho^{(A,B)}\Pi_{k}^{(B)}$
is the measured-projected total state and
$p_{k}=\tr(\Pi_{k}^{(B)}\rho^{(A,B)}\Pi_{k}^{(B)})$ is the associated
probability. The measurement and the associated projector $\Pi_{k}^{(B)}$ are
defined only in the sector $B$ of the bi-partition. For pure states, the
quantum discord reduces to the entanglement between subsystems with
$J(A,B)=\delta(A,B)$ Luo (2008). However, for mixed states this is not true in
general. This quantity is interesting since, as we will see in the next
sections, it can be a useful measure in order to study quantum phase
transitions in many-body systems Sarandy (2009); Allegra  _et al._ (2011);
Dillenschneider (2008).

However Eq. (2) requires a variational procedure involving all possible
$B$-subsystem projectors, so that computing quantum discord is in general
analytically and computationally intractable Huang (2014). Fortunately, if we
are dealing with fermionic systems, no optimization process is needed in order
to compute quantum discord between two arbitrary orbitals Faba  _et al._
(2021).

Other useful measure of system’s correlations is the overall entropy, defined
as

$S_{ov}=\sum_{i}S(\rho^{(i)})$

where $\rho^{(i)}$ is the reduced density matrix for the $i$-th orbital. Its
value is a measure of the total system’s correlations, if the total state is
pure Szalay  _et al._ (2015). It is closely related to the one-body entropy,
defined as the von Neumann entropy of the one-body density matrix Kanada-En’yo
(2015), whose elements are $\gamma_{ij}=\langle c^{\dagger}_{j}c_{i}\rangle$.
Because of the parity super-selection rule Wick  _et al._ (1952) we have

$\rho^{(i)}=\begin{pmatrix}1-\langle c^{\dagger}_{i}c_{i}\rangle&0\\\
0&\langle c^{\dagger}_{i}c_{i}\rangle\\\ \end{pmatrix}$

where the operator $c^{\dagger}_{i}$ ($c_{i}$) creates (annihilates) a
particle in the $i$-th orbital and the usual fermionic anticommutation rules
$\\{c^{\dagger}_{i},c_{j}\\}=\delta_{ij}$, $\\{c_{i},c_{j}\\}=0$ are
fulfilled. If the overall entropy is evaluated in the natural orbital basis
$\\{a^{\dagger}_{i}\\}$ (which is defined as the one that diagonalizes the
one-body density matrix, and it has been shown that is the basis that
minimizes the overall entropy Gigena and Rossignoli (2015)), then

$\begin{split}S_{ov}^{nat}&=\sum_{i}f(\langle a^{\dagger}_{i}a_{i}\rangle)\\\
S(\gamma)&=\sum_{i}g(\langle a^{\dagger}_{i}a_{i}\rangle)\end{split}$ (3)

where the functions $f$ and $g$ are defined as $f(x)=-(1-x)\log(1-x)-x\log(x)$
and $g(x)=-x\log(x)$. Since both $f$ and $g$ satisfy $f(0)=g(0)=f(1)=g(1)=0$
and are real valued smooth and strictly concave functions, the information and
behaviour of $S_{ov}^{nat}$ and $S(\gamma)$ are essentially the same.

As we will discuss in the following sections it will be useful to compare this
quantity, which quantifies the total system correlation (under a quantum
information perspective), with the relative correlation energy Löwdin (1955);
Wigner (1934), defined as

$\epsilon_{corr}=\frac{E_{exact}-E_{HF}}{E_{exact}}$

where $E_{exact}$ is the exact ground state energy and $E_{HF}$ is the ground
state energy obtained at the mean field (HF) level. Traditionally,
$\epsilon_{corr}$ has been used to quantify the amount of correlations in a
system, since it compares the exact ground state energy which contains all the
correlations in the system with the mean field one which is taken here as an
uncorrelated reference. Moreover, the correlation energy is closely related
with the overlap between the exact ground state and the Hartree-Fock one
Benavides-Riveros  _et al._ (2017).

##  III Two-level Lipkin model

In this section we will discuss the quantities previously defined, under the
context of the two-level Lipkin model.

The so called Lipkin model Lipkin  _et al._ (1965) (proposed by Lipkin,
Meshkov and Glick in 1964) consists of a $N$-fermion two level system
separated by an energy gap $\epsilon$, each level having a $N$-fold degeneracy
(we assume that all fermions are of the same type and have no spin, for
simplicity). We label the upper/lower level with the quantum number $\sigma=+$
or $\sigma=-$ respectively, and the degeneracy with the quantum number
$p=1,2,\ldots,N$. The quantum number $\sigma$ can also be interpreted as a
parity quantum number (see below). The Hamiltonian is given in terms of
fermionic creation and annihilation operators by

$H=\epsilon K_{0}-\frac{1}{2}V(K_{+}K_{+}+K_{-}K_{-})$ (4)

with

$\begin{split}K_{0}&=\frac{1}{2}\sum_{p=1}^{N}(c^{\dagger}_{+,p}c_{+,p}-c^{\dagger}_{-,p}c_{-,p})\\\
K_{+}&=\sum_{p=1}^{N}c^{\dagger}_{+,p}c_{-,p}\qquad
K_{-}=(K_{+})^{\dagger}\end{split}$

As the interaction is of the monopole-monopole type, the quantum number $p$ is
conserved in the model.

The advantage of this model is that it is exactly solvable, since the
operators introduced in Eq. (4) are the generators of the algebra of
$SU(2)$222See Refs. Robledo (1992); Di Tullio  _et al._ (2019) and references
therein for a detailed discussion of the exact solution..

The mean field (HF) solution can be easily obtained Lipkin  _et al._ (1965)
because the HF energy depends on a single variational parameter. Defining the
dimensionless interaction strength $\chi=\frac{(N-1)V}{\epsilon}$ it is
observed that, for certain values of $\chi$ the HF solution breaks the parity
symmetry of the Hamiltonian in Eq. (4) (to be associated with the $\sigma$
quantum number). With the above definitions, the parity operator is defined as

$P_{z}=\exp(i\pi\sum_{p=1}^{N}c^{\dagger}_{-,p}c_{-,p}).$

The HF states will have a well defined parity if they are eigenstates of the
parity operator $P_{z}$. When $\chi\leq 1$ the HF solution preserves the
parity symmetry (spherical phase) whereas the symmetry is broken (deformed
phase) when $\chi>1$ Robledo (1992). The correlation energy of the ground
state can be easily computed by comparing the exact solution with the Hartree-
Fock (HF) solution. This quantity as well as the the overall entropy in the
natural orbital basis depend on the strength parameter $\chi$ and they are
strongly correlated as can be seen in Fig. 1.

Figure 1: Relative correlation energy as a function of the overall entropy,
for the exact ground state of the two level Lipkin model and different
particle number $N$.

We can distinguish three regions in this plot. For low enough values of the
overall entropy, the relative correlation energy grows quasi-linearly. Then,
after a sudden discontinuity of the second derivative (see Fig. 2 below) the
relative correlation energy reaches a maximum and then bends down to gently
decrease until the overall entropy saturates. This change of tendency is due
to the phase transition observed in the HF solution at $\chi=1$. Fig. 1 can
also be interpreted in terms of the values of $\chi$. In the spherical phase
($\chi\leq 1$) the mean field solution catches as many correlations as
possible while preserving the non-interacting picture and preserving the
system’s symmetry. As the correlation/interaction grows (quantified by the
overall entropy/the parameter $\chi$) the relative correlation energy grows
too, showing that the mean field approach is less accurate since the
difference between $E_{HF}$ and $E_{exact}$ is getting bigger. When $\chi\geq
1$, the system’s correlations are too strong and the mean field solution
breaks the parity symmetry in order to catch as much as possible of them (see
Fig. 3). In this way, the relative correlation energy shows a decreasing
behaviour until the saturation of the overall entropy. This change in the
behaviour of the system can also be seen in Fig. 2, where the second
derivative of the relative correlation energy is plotted as a function of the
overall entropy. A sudden jump is observed in this quantity when $\chi=1$
signaling the quantum phase transition.

Figure 2: Second derivative of the relative correlation energy (Fig. 1) for
different values of $N$. We observe a discontinuity in $\chi=1$, related to
the phase transition of the model.

However, Fig. 2 hides some subtleties. Although the phase transition at
$\chi=1$ is clear by the presence of the discontinuity, Fig. 2 must not be
interpreted as a ‘genuine phase transition indicator’. In a genuine phase
transition we observe a change in the system’s behaviour which becomes more
evident as the size of the system grows. In Fig. 2 we see the opposite
behaviour: the discontinuity is less abrupt when the system’s size (the number
of particles) is higher. This is due to the nature of the HF approximation: it
is more accurate for higher values of $N$ Ring and Schuck (2004). In this way,
we observe in Fig. 1 lower values for the relative correlation energy as $N$
increases, and therefore the discontinuity in the second derivative is less
abrupt.

We conclude from the previous discussion that the phase transition present in
the HF solution is not only a ‘feature’ of the mean field method but it also
reflects a structural change in the exact wave function of the system. This
statement is consistent with the fact that the mean field approximation
becomes more and more accurate as the number of particles in the system
increase and also with the fact that the phase transition is better defined as
the number of particles increases. To exemplify the later results, we show in
Fig. 3 the averaged333The averaged overall entropy is simply the ‘overall
entropy per particle’, this is, $\frac{S_{ov}}{N}$. overall entropy (for the
exact solution) as a function of the interaction parameter $\chi$ for some
values of $N$.

Figure 3: Averaged one body entropy of the exact ground state as a function of
the parameter $\kappa$. As the number of particles increases, the transition
between the spherical phase (low correlation regime) and the deformed one
(high correlated regime) is sharper.

We observe a sudden change in the averaged one-body entropy when $\chi=1$,
which is sharper as the number of particles increases. As the value of the
averaged one-body entropy is a measure of the correlations in the system we
conclude that the spherical phase ($\chi<1$) corresponds to a low-correlated
regime in the exact solution, while the deformed phase ($\chi\geq 1$)
corresponds to a high-correlated regime. Therefore, the behaviour of the
overall entropy, which quantifies the total correlation, can help us to
distinguish between different phases in the exact solution.

Other interesting quantity related to the overall entropy, which is computed
from the mean field state, is the two-orbital quantum discord Faba  _et al._
(2021) between a couple of modes with same $p$ and opposite $\sigma$ for the
HF ground state. Because of the symmetries of this model, the reduced density
matrix of those modes are still pure (see Appendix A). For this reason, all
the quantum correlations are entanglement and the quantum discord reduces to
the entanglement entropy between modes. However, as we will see in Sec. IV,
this will not be the case for the three level Lipkin model. If we plot the
quantum discord Faba  _et al._ (2021)

$\delta(\sigma,p;-\sigma,p)=\begin{cases}0,&\text{if }\chi\leq 1\\\
h(\chi),&\text{if }1<\chi\end{cases}$ (5)

with
$h(x)=-\frac{1}{2}(1-\frac{1}{x})\ln\frac{1}{2}(1-\frac{1}{x})-\frac{1}{2}(1+\frac{1}{x})\ln\frac{1}{2}(1+\frac{1}{x})$
as a function of the interaction parameter $\chi$ we obtain Fig. 4.

Figure 4: Quantum discord between up and down levels of same quantum number
$p$ (same degeneration level) for the HF ground state.

As in Fig. 3, we see clearly the quantum phase transition at $\chi=1$. In
fact, Fig. 4 is very similar to Fig. 3 when the particle number is large. This
is to be expected as the two orbital reduced state is pure and therefore the
single orbital entropy represents the entanglement between the two orbitals.
Thus, the overall entropy is twice the sum of the entanglement between the
orbital pairs. On the other hand, if we use Eq. (3) and we take into account
that in the natural basis $\langle a^{\dagger}_{-p}a_{-p}\rangle+\langle
a^{\dagger}_{+p}a_{+p}\rangle=1$, then $S_{ov}^{nat}=2S(\gamma)$. For this
reason, Fig. 3 and 4 are almost the same in the limit $N\rightarrow{\infty}$.
It is relevant to note that the quantum discord depicted in Fig. 4 does not
depend on the particle number since it is a ‘microscopic’ quantity (i.e, it is
defined between a couple of orbitals) of a mean-field state. However, we can
see clearly the quantum phase transition in the behaviour of this quantity.
Moreover, as discussed in Di Tullio  _et al._ (2019), the nonzero quantum
discord (entanglement in this model) showed in Fig. 4 is a direct consequence
of the symmetry breaking at the mean field level: for the exact ground state,
the reduced density matrix for two levels with the same $p$ and opposite
$\sigma$ does not have coherent elements and therefore entanglement. However,
for the HF ground state, the reduced density matrix is pure and entangled.

##  IV Three level Lipkin model

This model is a generalization Li  _et al._ (1970) of the $N$-particle two
level Lipkin model discussed in the previous section. There are three energy
levels in the model, each one with a $N$-fold degeneracy and, analogously to
the two level Lipkin model, the interaction term can’t change the degeneracy
quantum number $p=1,2,...,N$. If we assume that the interaction is the same
for the three levels, which are equally spaced, we can write the Hamiltonian
as

$H=\epsilon(K_{22}-K_{00})-\frac{V}{2}(K^{2}_{10}+K^{2}_{20}+K^{2}_{21}+h.c)$
(6)

with

$K_{\sigma\sigma^{\prime}}=\sum_{p=1}^{N}c^{\dagger}_{\sigma
p}c_{\sigma^{\prime}p}$

As explained in Li  _et al._ (1970); Holzwarth and Yukawa (1974b); Hagino and
Bertsch (2000) the exact ground state of Eq. (6) can be easily computed
numerically in the basis $\ket{n_{1},n_{2}}$, where $n_{i}$ is the number of
particles in the $i$-th level. The basis elements are built upon the action of
the operators $K_{10}^{n_{1}}$ and $K_{20}^{n_{2}}$ acting on the states with
all the orbits in level 0 occupied. The given set of states is a basis to
diagonalize $H$ because the operators $K_{\sigma\sigma^{\prime}}$ are the
generators of the algebra of $SU(3)$.

Figure 5: Relative correlation energy as a function of the overall entropy,
for the exact ground state of the three level Lipkin model.

If we compute the HF solution of the three level Lipkin model Holzwarth and
Yukawa (1974b); Hagino and Bertsch (2000), it can be seen that there are two
phase transitions, each one corresponding to the breaking of a level’s
symmetry. More precisely, the first phase transition is located at $\chi=1$
and corresponds to a parity-like breaking of the $\sigma=1$ level, while the
second one is located in $\chi=3$ and corresponds to a parity-like breaking of
the $\sigma=2$ level. This behaviour is reflected in Fig. 5 where the relative
correlation energy as a function of the overall entropy, for the exact ground
state is depicted in a similar way as in Fig 1. When the system’s correlation
is low enough, the relative correlation energy grows quasi-linearly until
reaching the second derivative discontinuity at $\chi=1$ (Fig. 6). This is
required in order to catch the maximum correlations as possible while
maintaining the non-interacting ansatz. From this point on, the relative
correlation energy increases more slowly until the the second quantum phase
transition takes place at $\chi=3$. From there on, the relative correlation
energy decreases while the overall entropy increase reflecting the fact the
mean field solution approximates better the exact solution. Finally, as in
Fig. 1, the overall entropy saturates.

Figure 6: Second derivative of Fig. 5. We observe a discontinuity in $\chi=1$
and $\chi=3$ related to the two phase transitions of this model.

As discussed in Sec. III, the relative correlation energy acquire lower values
as the particle number increases, in agreement with the general idea that the
mean field picture increases its accuracy in the thermodinamic limit (infinite
number of particles). This is the reason why the discontinuity in the second
derivative depicted in Fig. 6 for different values of particle number $N$ is
less and less pronounced as $N$ gets higher and higher. The behaviour is
essentially the same as in the two level Lipkin model except for the double
quantum phase transition.

Following the same analysis as in Sec. III, we can study the quantum phase
transitions through the overall entropy as a function of the interaction
parameter (Fig. 7), and compare it with the quantum discord between levels
with different $\sigma$ and same $p$ (Fig. 8).

Figure 7: Averaged one body entropy of the exact ground state as a function of
the parameter $\chi$. Figure 8: Quantum discord between levels of the same
degeneration number $p$ for the HF ground state. The sum of the quantum
discord between the three different configurations is represented with the
solid red line.

The behaviour of the overall entropy is very similar in both models. If $N$ is
small (for example $N=5$ in Fig. 7) the shape of the overall entropy is almost
the same for the two and three Lipkin models: since a quantum phase transition
is a global property there is no difference between phases. However, as the
particle number increases, the distinction between the three regions
(spherical phase in $\chi\leq 1$, first parity break in $1\leq\chi\leq 3$ and
second parity break in $3\leq\chi$) is sharper, and the differences between
the two and three level Lipkin models arises. They can be clearly observed
through the quantum discord between levels of the same degeneration number $p$
for the HF ground state (Fig. 8)444Unlike in the two level Lipkin model, here
the reduced state is in general mixed, and we can’t compute the entanglement
as in Eq. (1).. For the spherical region, there is no quantum correlation
between any level since the HF orbitals are related to the original ones
through the identity matrix. That is, the mean field state is simply the non
interacting ground state of Eq. (6). When the first symmetry breaking occurs,
the quantum correlations between states with $\sigma=0$ and $1$ increases
abruptly, while it remains zero between $\sigma=0$ and $2$ and $\sigma=1$ and
2. Indeed, its value is exactly the same as the two level Lipkin model (see
Eq. (5) and (8)). Since the $\sigma=2$ level remains unfilled (the $\sigma=0$
and $1$ levels are mixed while the $\sigma=2$ level is not), there is no
difference between quantum correlations of the two and three level Lipkin
models within the mean field description. However, when the second symmetry
breaking occurs, the three $\sigma$ levels are completely mixed. The quantum
correlations between $\sigma=0$ and $1$ levels spontaneously decreases due to
the redistribution of the occupation between all levels, while the quantum
discord between $\sigma=0$ and $2$ and $\sigma=1$ and $2$ grows in a very
similar fashion (being the quantum correlations between $\sigma=1$ and $2$
always lower). Finally, it is interesting to note that if we compare the sum
of the quantum discord between the three possible orbital combinations (solid
red line in Fig. 8) with the one-body entropy in Fig. 7, we see that if $N$ is
high enough, both line’s shapes follow the same ‘double jump’ trend.

##  V Conclusions

The relative correlation energy has been typically used in order to quantify
the amount of correlation in a state, since it is defined as the relative
difference between the exact energy and the mean-field one. On the other hand,
with the fast growth in the last decades of the quantum information field,
there are a variety of methods nowadays in order to quantify the correlation
in a system in terms of their subsystems. An example are the entanglement
entropy, mutual information, quantum discord or the one-body entropy. In this
work we have analyzed the relative correlation energy and some quantum
information measures in the context of the two and three level Lipkin model.
We found that the relative correlation energy is not a good estimator of the
total correlation of a system, but it is a good estimator of the accuracy of
the mean field approximation. Comparing the overall entropy (which is a
measure of the total system’s correlation under the quantum information
context) and the relative correlation energy we don’t find quasi-linear or
monotonously increasing behaviour. Indeed, we find regions in the parameter
space in which the overall entropy grows but the relative correlation energy
tends to decrease, and regions in which both tend to grow. Those regions are
defined by quantum phase transitions, which can be analyzed by computing the
quantum discord between orbitals at HF level, without the need of computing
the exact ground state.

Future work includes the analysis of different models, both analytically or
numerically solvable, such as $N$-level Lipkin, picket fence or single-$j$
shell models. Also, a more exhaustive analysis can be performed in more
complex systems by computing the quantum discord between bigger orbital
subsystems of interest, or extending the mean field picture to a quasiparticle
vacuum.

###### Acknowledgements.

The authors want to thank the Madrid regional government, Comunidad Autónoma
de Madrid, for the project Quantum Information Technologies: QUITEMAD-CM
P2018/TCS-4342. The work of LMR was supported by Spanish Ministry of Economy
and Competitiveness (MINECO) Grants No. PGC2018-094583-B-I00. We would like to
thank Jorge Tabanera for enlightening discussions.

##  Appendix A  Purity of the two orbital reduced density matrix for the
Hartree-Fock ground state of the two level Lipkin model

In this section we will compute the purity of the two orbital reduced density
matrix for the HF ground state of the two level Lipkin model. As explained in
Davis and Heiss (1986)555Here the authors work with the Agassi model, which is
an extension of the two level Lipkin model. we can write the one-body density
matrix of the HF ground state as

$\gamma_{\sigma
p,\sigma^{\prime}p^{\prime}}=\begin{cases}\frac{1}{2}(1-\sigma\cos\varphi)\delta_{p,p^{\prime}},&\text{if
}\sigma=\sigma^{\prime}\\\
-\frac{1}{2}\sin\varphi\delta_{p,p^{\prime}},&\text{if
}\sigma=-\sigma^{\prime}\end{cases}$

with

$\cos\varphi=\begin{cases}1,&\text{if }\chi\leq 1\\\ \frac{1}{\chi},&\text{if
}\chi>1\end{cases}$

Following the results in Faba  _et al._ (2021), the two orbital reduced
density matrix is

$\rho^{(A,B)}=\frac{1}{2}\begin{pmatrix}0&0&0&0\\\
0&1+\cos\varphi&-\sin\varphi&0\\\ 0&-\sin\varphi&1-\cos\varphi&0\\\ 0&0&0&0\\\
\end{pmatrix}$

whose eigenvalues are $0$ and $1$.

##  Appendix B  Quantum discord for the Hartree-Fock state in the three level
Lipkin model

In this section we will briefly develop the analytic expression for the two
orbital quantum discord in the HF state of the three level Lipkin model.
Following reference Faba  _et al._ (2021), we only need to compute the one-
body elements and the two-body diagonal elements for each orbital. If we
assume that the system is in the HF ground state, i.e,
$\ket{HF}=\prod_{q=1}^{N}a^{\dagger}_{0q}\ket{0}$ (with $\ket{0}$ the vacuum
state) then, using Wick’s theorem,

$\begin{split}\bra{HF}c^{\dagger}_{\alpha i}c_{\beta
j}\ket{HF}&=U^{\dagger}_{\alpha 0}U_{0\beta}\delta_{ij}\\\
\bra{HF}c^{\dagger}_{\alpha i}c^{\dagger}_{\beta j}c_{\beta j}c_{\alpha
i}\ket{HF}&=\absolutevalue{U_{0\alpha}U_{0\beta}}^{2}(1-\delta_{ij})\end{split}$

with $a^{\dagger}_{\alpha
i}=\sum_{\beta=0}^{2}U_{\alpha\beta}c^{\dagger}_{\beta i}$ and
$UU^{\dagger}=1$. Following the results in Holzwarth and Yukawa (1974b);
Hagino and Bertsch (2000), the mean field solution can be written as

$U=\begin{pmatrix}c_{\alpha}&c_{\beta}s_{\alpha}&s_{\beta}s_{\alpha}\\\
-c_{\beta}s_{\alpha}&1+c^{2}_{\beta}(c_{\alpha}-1)&s_{\beta}c_{\beta}(c_{\alpha}-1)\\\
-s_{\beta}s_{\alpha}&s_{\beta}c_{\beta}(c_{\alpha}-1)&1+s^{2}_{\beta}(c_{\alpha}-1)\end{pmatrix}$
(7)

with $c_{\alpha}\coloneqq\cos\alpha$, $s_{\alpha}\coloneqq\sin\alpha$ and

$\cos^{2}\alpha=\begin{cases}1,&\text{if }\chi\leq 1\\\
\frac{1}{2}(1+\frac{1}{\chi}),&\text{if }1<\chi\leq 3\\\
\frac{\chi+3}{3\chi},&\text{if }3<\chi\\\ \end{cases}$
$\cos^{2}\beta=\begin{cases}1,&\text{if }1<\chi\leq 3\\\
\frac{1}{2}(\frac{3}{2\chi-3}+1),&\text{if }3<\chi\\\ \end{cases}$

Using those results and Eq. (5) in Faba  _et al._ (2021), we easily obtain the
analytic expression for the quantum discord between any orbital pair:

$\begin{split}\delta(0,p;1,p)=&\begin{cases}0,&\text{if }\chi\leq 1\\\
s(\frac{1}{2}(1+\frac{1}{\chi}))+s(\frac{1}{2}(1-\frac{1}{\chi})),&\text{if
}1<\chi\leq 3\\\
-s(\frac{2}{3}+\frac{1}{\chi})+s(\frac{1}{3}+\frac{1}{\chi})+s(\frac{1}{3}),&\text{if
}3<\chi\\\ \end{cases}\\\ \delta(0,p;2,p)=&\begin{cases}0,&\text{if }\chi\leq
3\\\
s(\frac{1}{3}+\frac{1}{\chi})+s(\frac{1}{3}-\frac{1}{\chi})-s(\frac{2}{3}),&\text{if
}3<\chi\\\ \end{cases}\\\ \delta(1,p;2,p)=&\begin{cases}0,&\text{if }\chi\leq
3\\\
-s(\frac{2}{3}-\frac{1}{\chi})+s(\frac{1}{3}-\frac{1}{\chi})+s(\frac{1}{3}),&\text{if
}3<\chi\\\ \end{cases}\end{split}$ (8)

with $s(x)=-x\log x$.

## References

  * Ring and Schuck (2004) P. Ring and P. Schuck,  _The Nuclear Many-Body Problem_, Physics and astronomy online library (Springer, 2004).
  * Bender  _et al._ (2003) M. Bender, P.-H. Heenen, and P.-G. Reinhard, Rev. Mod. Phys. 75, 121 (2003).
  * Robledo  _et al._ (2019) L. M. Robledo, T. R. Rodríguez, and R. R. Rodríguez-Guzmán, Journal of Physics G: Nuclear and Particle Physics 46, 013001 (2019).
  * Sheikh  _et al._ (2019) J. A. Sheikh, J. Dobaczewski, P. Ring, L. M. Robledo, and C. Yannouleas, “Symmetry restoration in mean-field approaches,” (2019), arXiv:1901.06992 [nucl-th] .
  * Caurier  _et al._ (2005) E. Caurier, G. Martínez-Pinedo, F. Nowacki, A. Poves, and A. P. Zuker, Rev. Mod. Phys. 77, 427 (2005).
  * Legeza and Sólyom (2006) O. Legeza and J. Sólyom, Phys. Rev. Lett. 96, 116401 (2006).
  * Szalay  _et al._ (2015) S. Szalay, M. Pfeffer, V. Murg, G. Barcza, F. Verstraete, R. Schneider, and Ö. Legeza, International Journal of Quantum Chemistry 115, 1342 (2015), https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.24898 .
  * Zeng  _et al._ (2014) G.-M. Zeng, L.-A. Wu, and H.-J. Xing, Scientific Reports 4, 6377 (2014).
  * Tichy  _et al._ (2011) M. C. Tichy, F. Mintert, and A. Buchleitner, Journal of Physics B: Atomic, Molecular and Optical Physics 44, 192001 (2011).
  * Kanada-En’yo (2015) Y. Kanada-En’yo, Progress of Theoretical and Experimental Physics 2015 (2015), 10.1093/ptep/ptv050, 043D04, https://academic.oup.com/ptep/article-pdf/2015/4/043D04/19302046/ptv050.pdf .
  * Kruppa  _et al._ (2021) A. T. Kruppa, J. Kovács, P. Salamon, and O. Legeza, Journal of Physics G: Nuclear and Particle Physics 48, 025107 (2021).
  * Lipkin  _et al._ (1965) H. Lipkin, N. Meshkov, and A. Glick, Nuclear Physics 62, 188 (1965).
  * Li  _et al._ (1970) S. Y. Li, A. Klein, and R. M. Dreizler, Journal of Mathematical Physics 11, 975 (1970), https://doi.org/10.1063/1.1665234 .
  * Holzwarth and Yukawa (1974a) G. Holzwarth and T. Yukawa, Nuclear Physics A 219, 125 (1974a).
  * Legeza and Sólyom (2003) O. Legeza and J. Sólyom, Phys. Rev. B 68, 195116 (2003).
  * Bañuls  _et al._ (2007) M.-C. Bañuls, J. I. Cirac, and M. M. Wolf, Phys. Rev. A 76, 022311 (2007).
  * Ding  _et al._ (2021) L. Ding, S. Mardazad, S. Das, S. Szalay, U. Schollwöck, Z. Zimborás, and C. Schilling, Journal of Chemical Theory and Computation 17, 79 (2021).
  * Friis  _et al._ (2013) N. Friis, A. R. Lee, and D. E. Bruschi, Phys. Rev. A 87, 022338 (2013).
  * Gigena and Rossignoli (2015) N. Gigena and R. Rossignoli, Phys. Rev. A 92, 042326 (2015).
  * Robin  _et al._ (2021) C. Robin, M. J. Savage, and N. Pillet, Phys. Rev. C 103, 034325 (2021).
  * Nielsen and Chuang (2011) M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ , 10th ed. (Cambridge University Press, USA, 2011).
  * Ollivier and Zurek (2001) H. Ollivier and W. H. Zurek, Phys. Rev. Lett. 88, 017901 (2001).
  * Luo (2008) S. Luo, Phys. Rev. A 77, 042303 (2008).
  * Sarandy (2009) M. S. Sarandy, Phys. Rev. A 80, 022108 (2009).
  * Allegra  _et al._ (2011) M. Allegra, P. Giorda, and A. Montorsi, Phys. Rev. B 84, 245133 (2011).
  * Dillenschneider (2008) R. Dillenschneider, Phys. Rev. B 78, 224413 (2008).
  * Huang (2014) Y. Huang, New Journal of Physics 16, 033027 (2014).
  * Faba  _et al._ (2021) J. Faba, V. Martín, and L. Robledo, Phys. Rev. A 103, 032426 (2021).
  * Wick  _et al._ (1952) G. C. Wick, A. S. Wightman, and E. P. Wigner, Phys. Rev. 88, 101 (1952).
  * Löwdin (1955) P.-O. Löwdin, Phys. Rev. 97, 1509 (1955).
  * Wigner (1934) E. Wigner, Phys. Rev. 46, 1002 (1934).
  * Benavides-Riveros  _et al._ (2017) C. L. Benavides-Riveros, N. N. Lathiotakis, C. Schilling, and M. A. L. Marques, Phys. Rev. A 95, 032507 (2017).
  * Robledo (1992) L. M. Robledo, Phys. Rev. C 46, 238 (1992).
  * Di Tullio  _et al._ (2019) M. Di Tullio, R. Rossignoli, M. Cerezo, and N. Gigena, Phys. Rev. A 100, 062104 (2019).
  * Holzwarth and Yukawa (1974b) G. Holzwarth and T. Yukawa, Nuclear Physics A 219, 125 (1974b).
  * Hagino and Bertsch (2000) K. Hagino and G. F. Bertsch, Phys. Rev. C 61, 024307 (2000).
  * Davis and Heiss (1986) E. D. Davis and W. D. Heiss, Journal of Physics G: Nuclear Physics 12, 805 (1986).