update

RapportStage/Rapport.bib

@@ -756,4 +756,40 @@
 	journaltitle = {The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science},
 	author = {Coulson, C. A. and Fischer, I.},
 	date = {1949-04-01},
+}
+
+@article{Burton_2019,
+	title = {Complex adiabatic connection: {A} hidden non-{Hermitian} path from ground to excited states},
+	volume = {150},
+	doi = {10.1063/1.5085121},
+	number = {4},
+	journal = {The Journal of Chemical Physics},
+	author = {Burton, Hugh G. A. and Thom, Alex J. W. and Loos, Pierre-François},
+	month = jan,
+	year = {2019},
+	pages = {041103},
+}
+
+@article{Burton_2019a,
+	title = {Parity-{Time} {Symmetry} in {Hartree}–{Fock} {Theory}},
+	volume = {15},
+	doi = {10.1021/acs.jctc.9b00289},
+	number = {8},
+	journal = {Journal of Chemical Theory and Computation},
+	author = {Burton, Hugh G. A. and Thom, Alex J. W. and Loos, Pierre-François},
+	month = aug,
+	year = {2019},
+	pages = {4374--4385},
+}
+
+@article{Hiscock_2014,
+	title = {Holomorphic {Hartree}–{Fock} {Theory} and {Configuration} {Interaction}},
+	volume = {10},
+	doi = {10.1021/ct5007696},
+	number = {11},
+	journal = {Journal of Chemical Theory and Computation},
+	author = {Hiscock, Hamish G. and Thom, Alex J. W.},
+	month = nov,
+	year = {2014},
+	pages = {4795--4800},
 }

RapportStage/Rapport.tex

@@ -287,7 +287,7 @@ In the Born-Oppenheimer approximation, the equation \eqref{eq:ExactHamiltonian}
     \bH=\sum\limits_{i=1}^{n}\left[ -\frac{1}{2}\grad_i^2 - \sum\limits_{k=1}^{N} \frac{Z_k}{|\vb{r}_i-\vb{R}_k|}+\sum\limits_{i<j}^{n}\frac{1}{|\vb{r}_i-\vb{r}_j|}\right]
 \end{equation}
 
-In the Hartree-Fock (HF) approximation the wave function is approximated as a single Slater determinant (which is an anti-symmetric combination of $n$ one electron spin-orbitals). Rather than solving the equation \eqref{eq:SchrEq}, the Hartree-Fock theory use the variational principle to find an approximation to $\Psi$. Hence the Slater determinants are not eigenfuctions of the exact Hamiltonian $\bH$. However they are eigenfunctions of an approximated Hamiltonian $\bH^{(0)}$, called the Hartree-Fock Hamiltonian, which is the sum of the one-electron Fock operators.
+In the Hartree-Fock (HF) approximation the wave function is approximated as a single Slater determinant (which is an anti-symmetric combination of $n$ one electron spin-orbitals). Rather than solving the equation \eqref{eq:SchrEq}, the Hartree-Fock theory uses the variational principle to find an approximation to $\Psi$. Hence the Slater determinants are not eigenfuctions of the exact Hamiltonian $\bH$. However, they are eigenfunctions of an approximated Hamiltonian $\bH^{(0)}$, called the Hartree-Fock Hamiltonian, which is the sum of the one-electron Fock operators.
 
 \begin{equation}\label{eq:HFHamiltonian}
 \bH^{(0)}= \sum\limits_{i=1}^{n} f(i)
@@ -313,7 +313,7 @@ K_b(1)\phi_a(1)=\left[\int\dd\vb{x}_2\phi_b^*(2)\frac{1}{r_{12}}\phi_a(2) \right
 The Hartree-Fock Hamiltonian \eqref{eq:HFHamiltonian} can be used as the zeroth-order Hamiltonian of the equation \eqref{eq:SchrEq-PT}. This partitioning of the Hamiltonian leads to the so-called M{\o}ller-Plesset (MP) perturbation theory \cite{Moller_1934}. The discovery of a partitioning of the Hamiltonian that allowed chemists to recover a part of the correlation energy (i.e. the difference between the exact energy and the Hartree-Fock energy) using perturbation theory has been a major step in the development of post-Hartree-Fock methods. This yields the Hamiltonian $\bH(\lambda)$ of the equation \eqref{eq:MPHamiltonian} where the two-electron part of the Fock operator $f(i)$ has been written $V_i^{HF}$ for convenience.
 
 \begin{equation}\label{eq:MPHamiltonian}
-    H(\lambda)=\sum\limits_{i=1}^{n}\left[-\frac{1}{2}\grad_i^2 - \sum\limits_{k=1}^{N} \frac{Z_k}{|\vb{r}_i-\vb{R}_k|} + (1-\lambda)V_i^{HF}+\lambda\sum\limits_{i<j}^{n}\frac{1}{|\vb{r}_i-\vb{r}_j|} \right]
+    \bH(\lambda)=\sum\limits_{i=1}^{n}\left[-\frac{1}{2}\grad_i^2 - \sum\limits_{k=1}^{N} \frac{Z_k}{|\vb{r}_i-\vb{R}_k|} + (1-\lambda)V_i^{\text{HF}}+\lambda\sum\limits_{i<j}^{n}\frac{1}{|\vb{r}_i-\vb{r}_j|} \right]
 \end{equation}
 
 In the perturbation theory the energy is a power series of $\lambda$ and the physical energy is obtained by taking $\lambda$ equal to 1. We will refer to the energy up to the $n$-th order as the MP$n$ energy. The MP0 energy overestimates the energy by double counting the electron-electron interaction, the MP1 corrects this effect and the MP1 energy is equal to the Hartree-Fock energy. The MP2 energy starts to recover a part of the correlation energy \cite{SzaboBook}.
@@ -372,7 +372,7 @@ Handy and co-workers exhibited the same behaviors of the series (oscillating and
 Cremer and He analyzed 29 atomic and molecular systems at the FCI level \cite{Cremer_1996} and regrouped all the systems in two classes. The class A systems which have a monotonic convergence to the FCI value and the class B which converge erratically after initial oscillations. The sample of systems contains stretched molecules and also some at equilibrium geometry, there are also some systems in various basis sets. They highlighted that \cite{Cremer_1996}\begin{quote}
 \textit{Class A systems are characterized by electronic structures with well-separated electron pairs while class B systems are characterized by electronic structures with electron clustering in one or more regions.}
 \end{quote}
-They proved that using different extrapolation formulas of the first terms of the M{\o}ller-Plesset of the series for class A and class B systems improve the precision of those formulas. The mean deviation from FCI correlation energies is 0.3 mhartrees whereas with the formula that do not distinguish the system it is 12 mhartrees. Even if there were still shaded areas and that this classification was incomplete, this work showed that understanding the origin of the different modes of convergence would lead to a more rationalized use of the M{\o]ller-Plesset perturbation theory and to more accurate correlation energies.
+They proved that using different extrapolation formulas of the first terms of the M{\o}ller-Plesset of the series for class A and class B systems improve the precision of those formulas. The mean deviation from FCI correlation energies is 0.3 mhartrees whereas with the formula that do not distinguish the system it is 12 mhartrees. Even if there were still shaded areas and that this classification was incomplete, this work showed that understanding the origin of the different modes of convergence would lead to a more rationalized use of the M{\o}ller-Plesset perturbation theory and to more accurate correlation energies.
 
 \subsection{Cases of divergence}
 
@@ -407,19 +407,23 @@ Moreover they proved that the extrapolation formulas of Cremer and He \cite{Crem
 
 In the 2000's Sergeev and Goodson \cite{Sergeev_2005, Sergeev_2006}  analyzed this problem from a more mathematical point of view by looking at the whole singularity structure where Olsen and his co-workers were trying to find the dominant singularity causing the divergence. They regrouped singularities in two classes: the $\alpha$ singularities which have unit order imaginary parts and the $\beta$ singularities which have very small imaginary parts. The singularities $\alpha$ are related to large avoided crossing between the ground state and a low-lying excited states. Whereas the singularities $\beta$ come from a sharp avoided crossing between the ground state and a highly diffuse state. They succeeded to explain the divergence of the series caused by $\beta$ singularities using a previous work of Stillinger \cite{Stillinger_2000}. 
 
-The M{\o}ller-Plesset Hamiltonian is defined as below and by reassembling the term we get the expression \eqref{eq:HamiltonianStillinger}.
+To understand the convergence properties at $\lambda=1$ of the perturbation series one need to look at the whole complex plane. In particular for real negative value of $\lambda$. If $\lambda$ is negative, the Coulomb interaction becomes attractive but the mean field stays repulsive as it is proportional to $(1-\lambda)$.
 
-The first two terms, the kinetic energy and the electron-nucleus attraction, form the mono-electronic core Hamiltonian which is independant of $\lambda$. The third term is the mean field repulsion of the Hartree-Fock calculation done to get $H_0$ and the last term is the Coulomb repulsion. If $\lambda$ is negative, the Coulomb interaction becomes attractive but the mean field stays repulsive as it is proportional to $(1-\lambda)$. If $\lambda$ becomes more and more negative the mean field becomes more and more repulsive so the nucleus cannot bind the electrons anymore because the electron-nucleus attraction is not scaled with $\lambda$. The repulsive mean field is localized around nucleus whereas the electrons interactions persist away from nucleus. There is a real negative value $\lambda_c$ where the electrons form a bound cluster and goes to infinity. According to Baker this value is a critical point of the system and by analogy with thermodynamics the energy $E(\lambda)$ exhibits a singularity at $\lambda_c$ \cite{Baker_1971}. At this point the system undergo a phase transition and a symmetry breaking. Beyond $\lambda_c$ there is a continuum of eigenstates with electrons dissociated from the nucleus.
+\begin{equation}\label{eq:HamiltonianStillinger}
+    \bH(\lambda)=\sum\limits_{i=1}^{n}\left[ \underbrace{-\frac{1}{2}\grad_i^2 - \sum\limits_{k=1}^{N} \frac{Z_k}{|\vb{r}_i-\vb{R}_k|}}_{\text{Independent of }\lambda} + \overbrace{(1-\lambda)V_i^{\text{HF}}}^{\textcolor{red}{Repulsive}}+\underbrace{\lambda\sum\limits_{i<j}^{n}\frac{1}{|\vb{r}_i-\vb{r}_j|}}_{\textcolor{blue}{Attractive}}  \right]
+\end{equation}
 
-This reasoning is done on the exact Hamiltonian and energy, this is the exact energy which exhibits this singularity on the negative real axis. But in finite basis set, one can prove that for a Hermitian Hamiltonian the singularities of $E(\lambda)$ occurs in complex conjugate pair with non-zero imaginary parts. Sergeev and Goodson proved, as predicted by Stillinger, that in a finite basis set the critical point on the real axis is modeled by a cluster of sharp avoided crossings with diffuse functions, equivalently by a cluster of $\beta$ singularities in the negative half plane. They explain that Olsen et al., because they used a $2\times2$ model, only observed the first singularity of this cluster of singularities causing the divergence.
+The major difference between those two terms is that the repulsive mean field is localized around nucleus whereas the electrons interactions persist away from nuclei. If $\lambda$ becomes more and more negative the mean field becomes more and more repulsive so the nuclei can not bind the electrons anymore because the electron-nucleus attraction is not scaled with $\lambda$. There is a real negative value $\lambda_c$ where the electrons form a bound cluster and goes to infinity. According to Baker this value is a critical point of the system and by analogy with thermodynamics the energy $E(\lambda)$ exhibits a singularity at $\lambda_c$ \cite{Baker_1971}. At this point the system undergo a phase transition and a symmetry breaking. Beyond $\lambda_c$ there is a continuum of eigenstates with electrons dissociated from the nucleus.
 
-Finally, it was shown that $\beta$ singularities are very sensitive to the basis sets but not to the stretching of the system. On the contrary $\alpha$ singularities are relatively insensitive to the basis sets but very sensitive to bond stretching. According to Goodson the singularity structure from molecules stretched from the equilibrium geometry is difficult \cite{Goodson_2004}, this is consistent with the observation of Olsen and co-workers on the \ce{HF} molecule at equilibrium geometry and stretched geometry \cite{Olsen_2000}. To our knowledge the effect of bond stretching on singularities, its link with spin contamination and symmetry breaking of the wave function hasn't been as well understood as the ionization effect and its link with diffuse function. In this work we try to improve our understanding of the effect of symmetry breaking on the singularities of $E(\lambda)$ and we hope that it will lead to a deeper understanding of perturbation theory.
+This reasoning is done on the exact Hamiltonian and energy, this is the exact energy which exhibits this singularity on the negative real axis. But in a finite basis set, one can prove that for a Hermitian Hamiltonian the singularities of $E(\lambda)$ occurs in complex conjugate pairs with non-zero imaginary parts. Sergeev and Goodson proved \cite{Sergeev_2005}, as predicted by Stillinger, that in a finite basis set the critical point on the real axis is modeled by a cluster of sharp avoided crossings with diffuse functions, equivalently by a cluster of $\beta$ singularities in the negative half plane. They explain that Olsen et al., because they used a $2\times2$ model, only observed the first singularity of this cluster of singularities causing the divergence.
+
+Finally, it was shown that $\beta$ singularities are very sensitive to changes of the basis set but not to the stretching of the system. On the contrary $\alpha$ singularities are relatively insensitive to the basis sets but very sensitive to bond stretching. According to Goodson the singularity structure from molecules stretched from the equilibrium geometry is difficult \cite{Goodson_2004}, this is consistent with the observation of Olsen and co-workers on the \ce{HF} molecule at equilibrium geometry and stretched geometry \cite{Olsen_2000}. To our knowledge the effect of bond stretching on singularities, its link with spin contamination and symmetry breaking of the wave function hasn't been as well understood as the ionization effect and its link with diffuse function. In this work we try to improve our understanding of the effect of symmetry breaking on the singularities of $E(\lambda)$ and we hope that it will lead to a deeper understanding of perturbation theory.
 
 \subsection{The physics of quantum phase transition}
 
 In the previous section, we saw that a reasoning on the Hamiltonian allows us to predict the existence of a critical point. In a finite basis set this critical point is model by a cluster of singularity $\beta$. It is now well-known that this phenomenon is a specific case of a more general phenomenon. Indeed, theoretical physicists proved that EPs close to the real axis are connected to quantum phase transitions \cite{Heiss_1988, Heiss_2002, Cejnar_2005, Cejnar_2007, Cejnar_2009, Borisov_2015, Sindelka_2017}. In quantum mechanics, the Hamiltonian is almost always dependent of at least one parameter, in some cases the variation of a parameter can lead to abrupt changes at a critical point. Those quantum phase transitions exist both for ground and excited states \cite{Cejnar_2009, Sachdev_2011, Cejnar_2015, Cejnar_2016, Caprio_2008, Macek_2019}. A ground-state quantum phase transition is characterized by the successive derivative of the ground-state energy with respect to a non-thermal control parameter \cite{Cejnar_2009, Sachdev_2011}. The transition is called discontinuous and of first order if the first derivative is discontinuous at the critical parameter value. Otherwise, it is called continuous and of n-th order if the n-th derivative is discontinuous. A quantum phase transition can also be identify by the discontinuity of an appropriate order parameter (or one of its derivative). 
 
-The presence of an EP close to the real axis is characteristic of a sharp avoided crossing. Yet at such an avoided crossing eigenstates change abruptly. Although it is now well understood that EPs are closely related to quantum phase transitions, the link between the type of QPT (ground state or excited state, first or superior order) and EPs still need to be clarify. One of the major challenge in order to do this resides in our ability to compute the distribution of EPs. The numerical assignment of an EP to two energies on the real axis is very difficult in large dimensions so methods are required to get information on the location of EPs. Cejnar et al. developed a method based on a Coulomb analogy giving access to the density of EP close to the real axis \cite{Cejnar_2005, Cejnar_2007}. More recently Stransky and co-workers proved that the distribution of EPs is not the same around a QPT of first or second order \cite{Stransky_2018}. Moreover, that when the dimension of the system increases they tends towards the real axis in a different manner, meaning respectively exponentially and algebraically.
+The presence of an EP close to the real axis is characteristic of a sharp avoided crossing. Yet at such an avoided crossing eigenstates change abruptly. Although it is now well understood that EPs are closely related to quantum phase transitions, the link between the type of QPT (ground state or excited state, first or superior order) and EPs still need to be clarified. One of the major challenges in order to do this resides in our ability to compute the distribution of EPs. The numerical assignment of an EP to two energies on the real axis is very difficult in large dimensions so methods are required to get information on the location of EPs. Cejnar et al. developed a method based on a Coulomb analogy giving access to the density of EP close to the real axis \cite{Cejnar_2005, Cejnar_2007}. More recently Stransky and co-workers proved that the distribution of EPs is not the same around a QPT of first or second order \cite{Stransky_2018}. Moreover, that when the dimension of the system increases they tends towards the real axis in a different manner, meaning respectively exponentially and algebraically.
 
 It seems like our understanding of the physics of spatial and/or spin symmetry breaking in the Hartree-Fock theory can be enlightened by quantum phase transition theory. Indeed, the second derivative of the energy is discontinuous at the Coulson-Fischer point which means that the system undergo a second order quantum phase transition. Moreover, the $\beta$ singularities introduced by Sergeev to describe the EPs modeling the formation of a bound cluster of electrons are actually a more general class of singularities. The EPs close to the real axis (the so-called $\beta$ singularities) are connected to QPT because they result from a sharp avoided crossings at which the eigenstates change quickly. However the $\alpha$ singularities arise from large avoided crossings therefore they can not be connected to QPT. The avoided crossings generating $\alpha$ singularities generally involve the ground state and low-lying doubly-excited states. Those excited states have a non-negligible contribution to the exact FCI solution because they have the same spatial and spin symmetry as the ground state. We think that $\alpha$ singularities are connected to the multi-reference behavior of the wave function in the same way as $\beta$ singularities are linked to symmetry breaking and phase transition of the wave function. 
 
@@ -440,10 +444,6 @@ The laplacian operators are the kinetic operators for each electrons and $r_{12}
 	\item Radius of the spherium that ultimately dictates the correlation regime.
 \end{itemize}
 
-\subsection{Weak correlation regime}
-
-\subsubsection{Restricted and unrestricted equation for the  spherium model}
-
 In the restricted Hartree-Fock formalism, the wave function cannot model properly the physics of the system at large R because the spatial orbitals are restricted to be the same. Then a fortiori it cannot represent two electrons on opposite side of the sphere. In the unrestricted formalism there is a critical value of R, called the Coulson-Fischer point \cite{Coulson_1949}, at which a second unrestricted Hartree-Fock solution appear. This solution is symmetry-broken as the two electrons tends to localize on opposite side of the sphere. By analogy with the case of \ce{H_2} \cite{SzaboBook}, the unrestricted Hartree-Fock wave function is defined as:
 
 \begin{equation}
@@ -481,8 +481,6 @@ E_{\text{p}} = \sum\limits_{i,j,k,l=0}^{\infty}C_{\alpha,i}C_{\alpha,j}C_{\beta,
 S_{i,j,k,l}=\sqrt{(2i+1)(2j+1)(2k+1)(2l+1)}
 \end{equation*}
 
-\subsubsection{The minimal basis example}
-
 We obtained the equation \eqref{eq:EUHF} for the general form of the wave function, but to be associated with a physical wave function the energy need to be stationary with respect to the coefficient. The general method is to use the Hartree-Fock self-consistent field method to get the coefficients of the wave functions corresponding to physical solutions. We will work in a minimal basis to illustrate the difference between the RHF and UHF solutions. In this basis there is a shortcut to find the stationary solutions. One can define the mono-electronic wave functions $\phi(\theta)$ using a mixing angle between the two basis functions. Hence we just need to minimize the energy with respect to the two mixing angles $\chi_\alpha$ and $\chi_\beta$.
 
 \begin{equation}
@@ -498,9 +496,13 @@ The minimization gives the 3 anticipated solutions (valid for all value of R):
 
 Vérifier minimum, maximum, point selle
 
-\subsubsection{Symmetry-broken solutions}
+In addition, there is also the well-known symmetry-broken UHF solution. For $R>3/2$ an other stationary UHF solution appear, this solution is a minimum of the Hartree-Fock equations. This solution corresponds to the configuration with the electron $\alpha$ on one side of the sphere and the electron $\beta$ on the opposite side and the other way round. The electrons can be on opposite sides of the sphere because the choice of $Y_{10}$ as a basis function induced a privileged axis on the sphere for the electrons. This solution have the energy \eqref{eq:EUHF} for $R>3/2$.
 
-In addition, there is also the well-known symmetry-broken UHF solution. For $R>3/2$ an other stationary UHF solution appear, this solution is a minimum. This solution corresponds to the configuration with the electron $\alpha$ on one side of the sphere and the electron $\beta$ on the opposite side and the other way round. The electrons can be on opposite sides of the sphere because the choice of $Y_{10}$ induced a privileged axis on the sphere for the electrons. The exact solution for the ground state is a singlet so this wave function does not have the true symmetry. However this solution gives more accurate results for the energy at large R as shown in \autoref{tab:ERHFvsEUHF}. In fact at the Coulson-Fischer point, it becomes more efficient to minimize the Coulomb repulsion than the kinetic energy leading to this symmetry breaking. There is a competition between those two effects: keeping the symmetry of the exact wave function and minimize the energy. This type of symmetry breaking is called a spin density wave because the system oscillate between the two symmetry-broken configurations.
+\begin{equation}
+E_{\text{sb-UHF}}=-\frac{75}{112R^3}+\frac{25}{28R^2}+\frac{59}{84R}
+\end{equation}
+
+The exact solution for the ground state is a singlet so this wave function does not have the true symmetry. However this solution gives more accurate results for the energy at large R as shown in \autoref{tab:ERHFvsEUHF}. In fact at the Coulson-Fischer point, it becomes more efficient to minimize the Coulomb repulsion than the kinetic energy leading to this symmetry breaking. This type of symmetry breaking is called a spin density wave because the system oscillate between the two symmetry-broken configurations.
 
 \begin{table}[h!]
 \centering
@@ -517,13 +519,22 @@ Exact & 9.783874 & 0.852781 & 0.391959 & 0.247898 & 0.139471 & 0.064525 & 0.0054
 \label{tab:ERHFvsEUHF}
 \end{table}
 
-There is also another symmetry-broken solution for $R>75/38$ but this one corresponds to a maximum. This solution is associated with another type of symmetry breaking somewhat well-known. Indeed it corresponds to a configuration where both electrons are on the same side of the sphere, in the same orbital. This solution is called symmetry-broken RHF. At the critical value of R, the repulsion of the two electrons on the same side of the sphere maximizes more the energy than the kinetic energy of the $Y_{10}$ orbitals. This symmetry breaking is associated with a charge density wave: the system oscillate between the situations with the electrons on each side.
+There is also another symmetry-broken solution for $R>75/38$ but this one corresponds to a maximum of the . This solution is associated with another type of symmetry breaking somewhat less known. It corresponds to a configuration where both electrons are on the same side of the sphere, in the same orbital. This solution is called symmetry-broken RHF. At the critical value of R, the repulsion of the two electrons on the same side of the sphere maximizes more the energy than the kinetic energy of the $Y_{10}$ orbitals. This symmetry breaking is associated with a charge density wave: the system oscillate between the situations with the electrons on each side.
 
-We can also consider negative value of R. This corresponds to the situation where one of the electrons is replaced by a positron. There are also a sb-RHF and a sb-UHF solution for some values of R but in this case the sb-RHF solution is a minimum and the sb-UHF is a maximum. Indeed, the sb-RHF state minimize the energy by placing the electron and the positron on the same side of the sphere. And the sb-UHF state maximize the energy because the two particles are on opposite side of the sphere.
+\begin{equation}
+E_{\text{sb-RHF}}=\frac{75}{88R^3}+\frac{25}{22R^2}+\frac{91}{66R}
+\end{equation}
 
-NRJ graphics
+We can also consider negative value of R. This corresponds to the situation where one of the electrons is replaced by a positron. There are also a sb-RHF ($R<-3/2$) and a sb-UHF ($R<-75/38$) solution for negative values of R (see Fig.\ref{fig:SpheriumNrj} but in this case the sb-RHF solution is a minimum and the sb-UHF is a maximum. Indeed, the sb-RHF state minimizes the energy by placing the electron and the positron on the same side of the sphere. And the sb-UHF state maximizes the energy because the two particles are on opposite sides of the sphere.
 
-\subsection{Strongly correlated regime}
+In addition, we can also consider the symmetry-broken solutions beyond their respective Coulson-Fischer by analytically continuing the respective energies leading to the so-called holomorphic solutions \cite{Hiscock_2014, Burton_2019, Burton_2019a}. All those energies are plot in \autoref{fig:SpheriumNrj}. The dotted curves corresponds to the holomorphic domain of the energies.
+
+\begin{figure}[h!]
+    \centering
+    \includegraphics[width=0.9\textwidth]{EsbHF.pdf}
+    \caption{\centering Energies of the 5 solutions of the HF equations (multiplied by $R^2$). The dotted curves correspond to the analytic continuation of the symmetry-broken solutions.}
+    \label{fig:SpheriumNrj}
+\end{figure}
 
 \section{Radius of convergence and exceptional points}
 
@@ -537,7 +548,7 @@ Size of the basis set
 
 Strong coupling ???
 
-\subsection{Exceptional points}
+\subsection{Exceptional points in the UHF formalism}
 
 RHF vs UHF