# Many-body perturbation theory for strongly correlated effective Hamiltonians using effective field theory methods

Raphaël Photopoulos¹,²*† and Antoine Boulet¹‡

1 Institut Supérieur des Matériaux du Mans, ISMANS CESI École d'ingénieurs, office 51, 44 Avenue Frédéric Auguste Bartholdi, 72000 Le Mans, France

2 IUT Grand Ouest Normandie, Université de Caen Normandie, Normandie Univ, 14000 Caen, France

* rphotopoulos@ismans.cesi.fr, † raphael.photopoulos@unicaen.fr, ‡ aboulet@ismans.cesi.fr

## Abstract

Introducing low-energy effective Hamiltonians is usual to grasp most correlations in quantum many-body problems. For instance, such effective Hamiltonians can be treated at the mean-field level to reproduce some physical properties of interest. Employing effective Hamiltonians that contain many-body correlations renders the use of perturbative many-body techniques difficult because of the overcounting of correlations. In this work, we develop a strategy to apply an extension of the many-body perturbation theory starting from an effective interaction that contains correlations beyond the mean-field level. The goal is to re-organize the many-body calculation to avoid the overcounting of correlations originating from the introduction of correlated effective Hamiltonians in the description. For this purpose, we generalize the formulation of the Rayleigh-Schrödinger perturbation theory by including free parameters adjusted to reproduce the appropriate limits. In particular, the expansion in the bare weak-coupling regime and the strong-coupling limit serves as a valuable input to fix the value of the free parameters appearing in the resulting expression. This method avoids double counting of correlations using beyond-mean-field strategies for the description of many-body systems. The ground state energy of various systems relevant for ultracold atomic, nuclear, and condensed matter physics is reproduced qualitatively beyond the domain of validity of the standard many-body perturbation theory. Finally, our method suggests interpreting the formal results obtained as an effective field theory using the proposed reorganization of the many-body calculation. The results, like ground state energies, are improved systematically by considering higher orders in the extended many-body perturbation theory while maintaining a straightforward polynomial expansion.

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

Copyright R. Photopoulos and A. Boulet.
This work is licensed under the Creative Commons Attribution 4.0 International License.
Published by the SciPost Foundation.

Received 2024-03-05
Accepted 2024-11-18
Published 2025-01-06
doi:10.21468/SciPostPhys.18.1.003

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

Check for
updates

## Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Presentation of the problem</td><td>3</td></tr><tr><td>2.1</td><td>Hartree-Fock or MBPT at first order</td><td>4</td></tr><tr><td>2.2</td><td>Direct application to ultracold atom systems</td><td>5</td></tr></table>
<table><tr><td>3</td><td>Extension of the MBPT with strongly correlated effective Hamiltonians</td><td>5</td></tr><tr><td>3.1</td><td>Many-body perturbation theory with a bare Hamiltonian</td><td>6</td></tr><tr><td>3.1.1</td><td>Rayleigh-Schrödinger perturbation theory</td><td>6</td></tr><tr><td>3.1.2</td><td>Alternative choice for the resolvent parameter</td><td>7</td></tr><tr><td>3.2</td><td>Extended many-body perturbation theory with an effective Hamiltonian</td><td>8</td></tr><tr><td>3.2.1</td><td>First order</td><td>8</td></tr><tr><td>3.2.2</td><td>Second order</td><td>8</td></tr><tr><td>3.2.3</td><td>Third order</td><td>9</td></tr><tr><td>3.2.4</td><td>Avoiding the double counting of correlations</td><td>10</td></tr><tr><td>3.3</td><td>Effective field theory interpretation of the extended MBPT</td><td>10</td></tr><tr><td>3.3.1</td><td>Choice of the high-scale solution</td><td>11</td></tr><tr><td>3.3.2</td><td>Remove the imaginary part</td><td>11</td></tr><tr><td>4</td><td>Quantum gases context: Application to ultracold atom systems</td><td>11</td></tr><tr><td>5</td><td>Nuclear physics context: Application to the Richardson pairing Hamiltonian</td><td>12</td></tr><tr><td>5.1</td><td>Exact GS energy</td><td>14</td></tr><tr><td>5.2</td><td>Hartree-Fock and MBPT solutions</td><td>14</td></tr><tr><td>5.2.1</td><td>Hartree-Fock</td><td>14</td></tr><tr><td>5.2.2</td><td>MBPT results</td><td>15</td></tr><tr><td>6</td><td>Condensed matter physics context: Application to the Hubbard Hamiltonian</td><td>16</td></tr><tr><td>6.1</td><td>Application to the 1D Hubbard chain model</td><td>17</td></tr><tr><td>6.2</td><td>Application to the four-sites 2D Hubbard model</td><td>17</td></tr><tr><td>7</td><td>Conclusion: Discussion and outlook</td><td>19</td></tr><tr><td colspan="2">References</td><td>20</td></tr></table>

# 1 Introduction

Quantum many-body problems are, most of the time, unsolvable exactly, and several approximations have been developed to grasp the properties of such complex systems. One can mention two types of approaches that aim to solve many-body problems under some well-controlled approximations: the non-perturbative methods such as Quantum Monte Carlo (QMC) [1–6] or Density Matrix Renormalization Group (DMRG) [7–10] and the perturbative techniques such as Many-Body Perturbation Theory (MBPT) [11–17] or Coupled-Cluster (CC) theory [18–22]. On the one hand, starting from the bare Hamiltonian of the theory, non-perturbative methods can express the problem with multidimensional integrals in the Feynman’s path integral formalism, which are manageable by Monte Carlo algorithms [23–26] or, alternatively, the problem can be tackled by optimizing a matrix product state tensor network using an iterative eigensolver such as the Lanczos algorithm [27, 28]. On the other hand, perturbative techniques rely on a truncation scheme up to an arbitrary order of the exact solution expressed as an infinite series [29, 30]. In particular, truncation at the leading order, i.e. at the first order in perturbation, is the so-called mean-field approximation.

In both perturbative and non-perturbative approaches, the description includes Beyond-Mean-Field (BMF) correlations. Another approach, namely the Density Functional Theory,<fnref n="1" /> is a third class of methods used to solve quantum many-body problems involving BMF correlations in an effective Hamiltonian formally treated with the mean-field approximation [31–34].

<footnote n="1">¹Also called Energy Density Functional (EDF) theory in the nuclear physics context.</footnote>
Difficulties arise if we apply perturbative methods beyond the mean-field level from an effective Hamiltonian that contains such BMF correlations by construction. Considering an effective Hamiltonian and employing many-body methods will include BMF correlations already accounted for in the effective Hamiltonian. Several methods have been developed to avoid this overcounting in many-body techniques to deal with effective interactions at the BMF level [35]. It seems like the effective interaction to consider depends on the order of perturbation considered. One possibility to avoid overcounting of correlations is to adjust the parameter of the effective Hamiltonian according to the truncation order considered to reproduce physical properties [36]. Another strategy consists of adding counter terms in the effective interaction to avoid double counting of correlations [37].

Starting from this observation, namely the dependence of the effective interaction as a function of the perturbative order considered in the many-body technique employed, the goal of this paper is to illustrate:

1. How to incorporate properly beyond mean-field correlations into an effective interaction, e.g. using limiting constraints, symmetries, etc.?

2. How to go beyond mean-field calculations using effective interactions that already contain such correlations?

We restrict our development to the case of the MBPT to clarify the formal link between the effective interaction and the truncation order of the perturbative method. The purpose of this work is then to develop a systematic MBPT for strongly correlated effective Hamiltonians.

The article is organized as follows. In section 2, the framework of our development is briefly presented and motivated to generalize the Rayleigh-Schrödinger perturbation theory. In section 3, we develop a reformulation of the MBPT exhibiting a new extra parameter that allows the use of an effective Hamiltonian that is the core of the presented work. We also discuss the interpretation of our extended MBPT in the Effective Field Theory (EFT) framework. Then, in sections 4, 5 and 6, we apply our strategy for several many-body problems encountered, respectively, in the context of ultracold atomic quantum gases, nuclear and condensed matter physics, as a proof of principle. Finally, in section 7, we conclude and discuss the implications and further developments of our approach for the description of quantum many-body systems.

## 2 Presentation of the problem

We consider a many-particle system of N fermions interacting through a Hamiltonian parameterized by a tunable coupling constant $\lambda$ and decomposed as:

$$ \hat{H}_{\lambda} = \hat{H}_{0} + \lambda \hat{h}. \qquad (1) $$

The reference state $|\Phi_0\rangle$ is chosen as the solution of the static many-body Schrödinger equation involving $\hat{H}_0$, that is, $\hat{H}_0|\Phi_0\rangle = E_0|\Phi_0\rangle$ with the unperturbed energy $E_0$. In this work, we only consider the normal component, making the reference state $|\Phi_0\rangle$ the Hartree-Fock Slater determinant associated with $\hat{H}_0$.

Then, we assume that we know the Taylor expansion of the Ground State (GS) energy around $\lambda = 0$ and its value at the limit $\lambda \to \infty$. These two limits, respectively $\lambda \ll 1$ and $\lambda \gg 1$, will be named respectively low- and high-scale limits in the following. For simplicity, we use dimensionless GS energy. We then write:

$$ \frac{E_{\lambda}}{E_{0}} = 1 + \sum_{n=1}^{\infty} \gamma_{n} \lambda^{n}, \qquad (2) $$
where we have introduced the set {$γ_n$} for the low-scale expansion of the dimensionless GS energy. In the high-scale limit, we define the parameter $ξ_0$ as:

$$ \lim_{\lambda \to \infty} \frac{E_\lambda}{E_0} = \xi_0. \qquad (3) $$

At the standard Hartree-Fock (HF) level or first-order MBPT, i.e. when $λ\hat{h}$ can be considered in perturbation according to $\hat{H}_0$ (low-scale limit $λ \ll 1$), the energy of the system is given by the parameter $γ_1$ only: $E_λ/E_0 = 1 + γ_1λ + O(λ^2)$. In general, the result of any order will give a polynomial expansion in $λ$ that diverges in the strong-coupling regime (high-scale limit $λ \gg 1$). Thus, in this sense, the standard MBPT does not provide a finite value of the observable in this limit. The goal of this paper is to propose a method based on the MBPT at a given order and relying on the use of an effective parameterization of the Hamiltonian such that:

(i) we recover the low-scale limit up to the order of the MBPT considered;

(ii) we impose the high-scale results $E_∞/E_0 = ξ_0$, where the parameter $ξ_0$ is assumed to be known.

In the following, we will arbitrarily assume the particular form of the effective Hamiltonian as follows:

$$ \tilde{H}_{\lambda} = \hat{H}_{0} + \frac{\lambda \hat{h}}{1 + \lambda \gamma_{1}/a}, \qquad (4) $$

which mimic, in a simple manner,<fnref n="2" /> the bare Hamiltonian (1) up to first order in power of $λ$, i.e. $\tilde{H}_λ = \hat{H}_λ + O(λ^2)$, but having a finite limit in the high-scale regime tuned by the parameter $a$. In other words, the bare coupling constant $λ$ is replaced by $λ/(1 + λγ_1/a)$ which can be interpreted as a renormalized coupling constant due to the medium effect [38]. This extra parameterization allows one to include many-body BMF correlations, e.g. the high-scale constraint (3), directly in the effective interaction. Using this particular form of effective Hamiltonians, we will propose a systematic method to perform MBPT calculations capable of reproducing the low-scale expansion of the GS energy (2) up to a given order in perturbation and the high-scale limit (3).

In the following, we will first illustrate our motivations with a recent example of the GS energy of diluted many-body Fermi systems at the HF level using an effective Hamiltonian equivalent to (4). Then, we extend the method to higher orders in perturbation using a reformulation of the MBPT that allows flexibility to adjust some parameters on the proper low- and high-scale constraints.

## 2.1 Hartree-Fock or MBPT at first order

The HF energy using the effective Hamiltonian (4) is given by the Leading Order (LO) in perturbation:

$$ E_{\lambda}^{LO} = E_{0} + \frac{\langle \Phi_{0} | \lambda \hat{h} | \Phi_{0} \rangle}{1 + \lambda \gamma_{1} / a}. \qquad (5) $$

Identifying with the low-scale expansion (2), we have $\langle \Phi_0 | \lambda \hat{h} | \Phi_0 \rangle = γ_1λE_0$. Then, imposing the exact high-scale limit constraint (3), we get $ξ_0 = 1 + a$. This expression of the high-scale parameter is valid only at the HF level and depends on the form chosen to define the effective Hamiltonian (4).

<footnote n="2">²In general, $\tilde{H}_λ = \hat{H}_0 + F_a(λ)\hat{h}$ where the function $F_a$ has the following properties: $F_a(λ) = λ + O(λ^2)$ and $F_a(∞) = a/γ_1$.</footnote>
## 2.2 Direct application to ultracold atom systems

In this section, as a proof of principle, we consider an infinite spin-saturated diluted system of fermions at zero temperature. In that case, the Hamiltonian is the sum of the kinetic Hamiltonian and a two-body contact interaction of the form: $\langle k'|\lambda\hat{h}|k \rangle = 4\pi a_s/m$ [39] where $a_s$ is the s-wave scattering length and plays the role of the parameter $\lambda$. Here, we use the convention $\hbar = 1$. We know that within the framework of the MBPT applied to dilute Fermi gases with an attractive interaction ($a_s < 0$), this corresponds to the Bardeen-Cooper-Schrieffer (BCS) regime [40, 41]:

$$ \frac{E}{E_0} = 1 + \frac{10}{9\pi}(a_s k_F) + \mathcal{O}[(a_s k_F)^2], \qquad (6a) $$

$$ \lim_{|a_s| \to \infty} \frac{E}{E_0} = \xi, \qquad (6b) $$

where $\xi = 0.376(4)$ [42] is the so-called Bertsch parameter, $k_F$ is the Fermi momentum of the system directly linked to the density $\rho = k_F^3/3\pi^2$, and $E_0/N = 3k_F^2/10m$ is the free Fermi gas energy.

Then we can directly identify and make the correspondence $\lambda = a_s$, $\gamma_1 = 10k_F/9\pi$ and $\xi_0 = \xi$ that leads to the HF energy (5):

$$ \frac{E_{0}^{LO}}{E_{0}} = 1 + \frac{\frac{10}{9\pi}(a_{s}k_{F})}{1 + \frac{10}{9\pi}(\xi - 1)^{-1}(a_{s}k_{F})} \qquad (7) $$

This expression was proposed several times during the last decade in the context of ultracold atoms using a Padé[1/1] form of the energy and the unitary limit as a constraint [43, 44]. In figure 1, we display the GS energy of the spin-saturated Fermi gas at the HF level, denoted as $E^{LO}$, as a function of ($a_s k_F$) and obtained with (7) using the correspondences of parameters for the dilute Fermi gas. As discussed in [43, 44], this functional qualitatively reproduces well the GS energy and thermodynamical properties of diluted ultracold atom systems across the weak to the strong coupling regimes in a very compact form, explicitly in terms of the density and the low-energy constant of the bare two-body interaction.

Including higher-order MBPT contributions (or applying more high-scale constraints such as the Taylor expansion in the limit $\lambda \to \infty$) at the HF level is not feasible in a systematic manner. For instance, various parameterizations of the effective Hamiltonian can lead to divergences in certain cases. Such developments are not the purpose of this work. Below, we suggest pursuing BMF calculations keeping the idea and the form of the Hamiltonian discussed up to now. In the following, we develop a method that allows us to extend MBPT for an effective Hamiltonian similar to (4) and keep the low- and high-scale limits valid. The aim is to elaborate a method that avoids naturally the double counting of correlations arising from the association of perturbative techniques and effective Hamiltonians.

## 3 Extension of the MBPT with strongly correlated effective Hamiltonians

In this section, we propose an extension of the MBPT allowing the use of effective Hamiltonians (4). We first recall the equation of MBPT in Rayleigh-Schrödinger formalism and then propose another truncation scheme by introducing a free parameter. This new formulation of the theory allows for some flexibility to perform the MBPT up to a given order using the effective Hamiltonian (4) and keep the proper low-scale expansion up to the order considered, as well as the high-scale constraint.
## 3.1 Many-body perturbation theory with a bare Hamiltonian

Assuming that the GS of the Hamiltonian $\hat{H}_{\lambda} = \hat{H}_{0} + \lambda \hat{h}$ for the many-body system of interest is nondegenerate [12, 45], we decompose the exact many-body GS $|\Psi_{0}\rangle$ as a series of reference states, e.g. Slater determinants, as:

$$|\Psi_0\rangle = |\Phi_0\rangle + \sum_{i \neq 0} C_i |\Phi_i\rangle, \qquad (8)$$

where $|\Phi_0\rangle$ is the GS of the unperturbed system, i.e. $\hat{H}_0 |\Phi_0\rangle = E_0 |\Phi_0\rangle$, $E$ is the exact GS energy of the system, the set $\{|\Phi_i\rangle\}$ denotes the eigenvector basis of the unperturbed Hamiltonian $\hat{H}_0$ (note that we use the intermediate normalization $\langle\Psi_0|\Phi_0\rangle = 1$), and the $\{C_i\}$ are the associated amplitudes of the state in this basis. The static many-body Schrödinger equation reads:

$$\hat{H}_{\lambda} |\Psi_{0}\rangle = E |\Psi_{0}\rangle, \qquad (9)$$

and using the hermiticity of the Hamiltonian, we define the energy shift as:

$$\Delta E = E - E_0 = \langle \Phi_0 | \lambda \hat{h} | \Psi_0 \rangle. \qquad (10)$$

Then, we define the projectors $\hat{P} = |\Phi_0\rangle\langle\Phi_0|$ and $\hat{Q} = \sum_{i \neq 0} |\Phi_i\rangle\langle\Phi_i|$ such that $|\Psi_0\rangle = (\hat{P} + \hat{Q})|\Psi_0\rangle = |\Phi_0\rangle + \hat{Q}|\Psi_0\rangle$. Now, starting from the Schrödinger equation (9), we introduce an energy operator $\hat{\omega}$ and we project the result into the sub-Hilbert space defined by $\hat{Q}$ to get:

$$|\Psi_0\rangle = |\Phi_0\rangle + \hat{R}(\hat{\omega}) \cdot (\hat{\omega} - E + \lambda \hat{h}) |\Psi_0\rangle, \qquad (11)$$

where we denote the resolvent operator (that is assumed to exist) by:

$$\hat{R}(\hat{\omega}) = \frac{\hat{Q}}{\hat{\omega} - \hat{H}_0}. \qquad (12)$$

By consecutive iterations, we can write the energy shift (10) and the exact many-body state (11) as a perturbative series:

$$\Delta E = \sum_{n=1}^{\infty} \langle \Phi_0 | \lambda \hat{h} \cdot [\hat{R}(\hat{\omega}) \cdot (\hat{\omega} - E + \lambda \hat{h})]^{n-1} | \Phi_0 \rangle, \qquad (13a)$$

$$|\Psi_0\rangle = \sum_{n=1}^{\infty} [\hat{R}(\hat{\omega}) \cdot (\hat{\omega} - E + \lambda \hat{h})]^{n-1} |\Phi_0\rangle. \qquad (13b)$$

### 3.1.1 Rayleigh-Schrödinger perturbation theory

In Rayleigh-Schrödinger Perturbation Theory (RSPT), we chose $\hat{\omega} = E_0$. Using that particular choice, the energy shift (13a) and the exact many-body state (13b) can be written as a perturbative series:

$$\Delta E = \sum_{n=1}^{\infty} \langle \Phi_0 | \lambda \hat{h} \cdot [\hat{R}(E_0) \cdot (\lambda \hat{h} - \Delta E)]^{n-1} | \Phi_0 \rangle, \qquad (14a)$$

$$|\Psi_0\rangle = \sum_{n=1}^{\infty} [\hat{R}(E_0) \cdot (\lambda \hat{h} - \Delta E)]^{n-1} |\Phi_0\rangle. \qquad (14b)$$
Then, expanding the energy shift as a series $\Delta E = \sum_n \Delta E^{(n)}$ in powers of $\lambda$ and using the fact that $\hat{Q}\Delta E|\Phi_0\rangle = 0$, we can get the expression of $\Delta E$ as a series in terms of the interaction strength $\lambda$. For instance, up to the third order in power of $\lambda$:

$$ \Delta E^{(1)} = \langle \Phi_0 | \lambda \hat{h} | \Phi_0 \rangle, \qquad (15a) $$

$$ \Delta E^{(2)} = \langle \Phi_0 | \lambda \hat{h} \cdot \hat{R}_0 \cdot \lambda \hat{h} | \Phi_0 \rangle, \qquad (15b) $$

$$ \Delta E^{(3)} = \langle \Phi_0 | \lambda \hat{h} \cdot \hat{R}_0 \cdot \lambda \hat{w} \cdot \hat{R}_0 \cdot \lambda \hat{h} | \Phi_0 \rangle, \qquad (15c) $$

where we have defined $\hat{R}_0 = \hat{R}(E_0)$ and $\lambda \hat{w} = \lambda \hat{h} - \Delta E^{(1)}$. We can remark that by matching the low-scale expansion for the energy in series of $\lambda^n$ given by (2), we have, by definition, $\Delta E^{(n)} = \gamma_n \lambda^n E_0$.

### 3.1.2 Alternative choice for the resolvent parameter

The choice (12) for the resolvent operator is not unique. For example, the Brillouin-Wigner perturbation theory consists of setting $\hat{\omega} = E$.

In this work, we propose to set $\hat{\omega} = \hat{H}_0 + \beta(E_0 - \hat{H}_0) = \hat{\omega}_{\beta}$, such that:

$$ \hat{R}(\hat{\omega}_{\beta}) = \frac{1}{\beta} \frac{\hat{Q}}{E_0 - \hat{H}_0} = \frac{1}{\beta} \hat{R}_0, \qquad (16) $$

where $\beta$ is a number to be determined. Note that the case $\beta = 1$ is equivalent to RSPT. The energy shift (13a) and the exact many-body state (13b) are now given by:

$$ \Delta E(\beta) = \sum_{n=1}^{\infty} \langle \Phi_0 | \lambda \hat{h} \cdot \hat{B}_n(\beta) | \Phi_0 \rangle, \qquad (17a) $$

$$ |\Psi_0(\beta)\rangle = \sum_{n=1}^{\infty} \hat{B}_n(\beta) |\Phi_0\rangle, \qquad (17b) $$

where we have defined $\hat{B}_1(\beta) = 1$ and for $n \ge 2$:

$$ \hat{B}_n(\beta) = \frac{1}{\beta^{n-1}} [(\beta - 1) + \hat{R}_0 \cdot (\lambda \hat{h} - \Delta E)]^{n-2} \cdot \hat{R}_0 \cdot \lambda \hat{h}. \qquad (17c) $$

To obtain these expressions, we use $(\hat{H}_0 - E_0)|\Phi_0\rangle = 0$, $\hat{Q}\Delta E|\Phi_0\rangle = 0$, the relationship $\hat{R}_0 \cdot (E_0 - \hat{H}_0) \cdot \hat{R}_0 = \hat{R}_0$ since $\hat{Q}$ is idempotent ($\hat{Q}^2 = \hat{Q}$), and the commutation property of $\hat{H}_0$ with $\hat{R}_0$.

We can show that, following the same truncation scheme in terms of the power of $\lambda$ as in RSPT, each term is re-summed exactly to be $\beta$-independent at each order and consequently $\Delta E(\beta) = \Delta E$. More precisely, we have:

$$ \begin{aligned} \sum_{n=2}^{\infty} \hat{B}_n(\beta) &= \frac{1}{1 - \hat{R}_0 \cdot (\lambda \hat{h} - \Delta E)} \cdot \hat{R}_0 \cdot \lambda \hat{h} \\ &= \sum_{n=2}^{\infty} [\hat{R}_0 \cdot (\lambda \hat{h} - \Delta E)]^{n-2} \cdot \hat{R}_0 \cdot \lambda \hat{h}, \end{aligned} $$

and we then recover the RSPT results (14), i.e. the $\beta$-independent results.

Considering the bare interaction Hamiltonian $\lambda \hat{h}$, the only (or more natural) possibility is to use the RSPT formalism, that is to set $\beta = 1$ or to use the simplified $\beta$-independent result above, which appears naturally in calculations. However, using an effective Hamiltonian like (4), the flexibility of the $\beta$ parameter can be used to match the low-scale expansion at a given order in MBPT. In the following, we illustrate this aspect at the first, second, and third orders using the effective Hamiltonian (4).
## 3.2 Extended many-body perturbation theory with an effective Hamiltonian

We consider the effective Hamiltonian (4) where the parameter $a$ is to be determined by containing the GS energy<fnref n="3" /> at N¹LO on the high-scale limit constraint (3). Using the preliminary result on the MBPT scheme proposed in section 3.1.2, the energy shift (13a) and the exact many-body state (13b) are now given by:

$$ \Delta E_{\lambda}(\beta) = \sum_{n=1}^{\infty} \frac{\langle \Phi_0 | \lambda \hat{h} \cdot \tilde{B}_n(\beta) | \Phi_0 \rangle}{1 + \lambda \gamma_1 / a}, \qquad (18a) $$

$$ |\Psi_0(\beta)\rangle = \sum_{n=1}^{\infty} \tilde{B}_n(\beta) |\Phi_0\rangle, \qquad (18b) $$

where we have defined $\tilde{B}_1(\beta) = 1$ and for $n \ge 2$:

$$ \tilde{B}_n(\beta) = \frac{1}{\beta^{n-1}} \left[ (\beta - 1) + \hat{R}_0 \cdot \left( \frac{\lambda \hat{h}}{1 + \lambda \gamma_1 / a} - \Delta E_{\lambda}(\beta) \right) \right]^{n-2} \cdot \hat{R}_0 \cdot \frac{\lambda \hat{h}}{1 + \lambda \gamma_1 / a} \qquad (18c) $$

The basic idea of performing MBPT calculations using the particular choice (4) for an effective Hamiltonian is to reorganize the standard truncation scheme and adjust the $\beta$ parameter to match the low-scale MBPT expansion. More precisely, we will take advantage of the reformulation of the RSPT to express the GS energy and the many-body state at N¹LO as series of the form:

$$ \frac{E_{\lambda}^{N^1 LO}}{E_0} = 1 + \sum_{n=1}^{l+1} \left( \frac{\lambda}{f_{\lambda}^{N^1 LO}} \right)^n \Gamma_n^{N^1 LO}, \qquad (19a) $$

$$ |\Psi_0^{N^1 LO}\rangle = \sum_{n=1}^{l+1} \left( \frac{\lambda}{f_{\lambda}^{N^1 LO}} \right)^{n-1} \hat{G}_n^{N^1 LO} |\Phi_0\rangle, \qquad (19b) $$

where $f_{\lambda}^{N^1 LO} = 1 + \lambda \gamma_1 / a$. In the following, we illustrate the strategy in the first ($l=0$), second ($l=1$), and third ($l=2$) order, and the parameter $a$ will now be denoted $a_{l+1}$ to avoid confusion.

### 3.2.1 First order

In that case, the GS energy is given by the energy shift (18) truncated at $l=0$, that is:

$$ \frac{E_{\lambda}^{LO}}{E_0} = 1 + \gamma_1 \frac{\lambda}{f_{\lambda}^{LO}}, \qquad (20) $$

where we can identify $\Gamma_1^{LO} = \gamma_1$ using the fact that $\langle \Phi_0 | \lambda \hat{h} | \Phi_0 \rangle = \gamma_1 \lambda E_0$. It remains the determination of the high-scale parameter using (3). As above, this gives $\xi_0 = 1 + a_1$.

### 3.2.2 Second order

Following the same strategy, we truncate (18) to $l=1$. We have first to determine:

$$ \frac{\langle \Phi_0 | \lambda \hat{h} \cdot \tilde{B}_2(\beta) | \Phi_0 \rangle}{1 + \lambda \gamma_1 / a_2} = \frac{1}{\beta} \frac{\langle \Phi_0 | \lambda \hat{h} \cdot \hat{R}_0 \cdot \lambda \hat{h} | \Phi_0 \rangle}{(1 + \lambda \gamma_1 / a_2)^2} \\ = \frac{\gamma_2}{\beta} \left( \frac{\lambda}{1 + \lambda \gamma_1 / a_2} \right)^2 E_0. \qquad (21) $$

<footnote n="3">³The notation N¹LO refer to the next-to-next-to-next-to-next leading order. For example, the next-to-leading order (NLO) corresponds to the second order in perturbation, the next-to-next-to-leading-order (N²LO) corresponds to the third order in perturbation, etc.</footnote>
Then, the GS energy at NLO reads:

$$ \frac{E_{\lambda}^{NLO}}{E_0} = 1 + \gamma_1 \frac{\lambda}{f_{\lambda}^{NLO}} + \frac{\gamma_2}{\beta} \left( \frac{\lambda}{f_{\lambda}^{NLO}} \right)^2. \qquad (22) $$

**Low-scale expansion** Now, we determine the $\beta$ parameter to match the low-scale expansion up to the second order, i.e. such that:

$$ \frac{E_{\lambda}^{NLO}}{E_0} = 1 + \gamma_1 \lambda + \gamma_2 \lambda^2 + \mathcal{O}(\lambda^3), \qquad (23) $$

leading to:

$$ \frac{1}{\beta} = 1 + \frac{\gamma_1^2}{a_2 \gamma_2}. \qquad (24) $$

**High-scale constraint** Finally, we fix the high-scale parameter $a_2$ using (3) leading to the quadratic equation:

$$ \xi_0 = 1 + 2a_2 + \frac{\gamma_2}{\gamma_1^2} a_2^2. \qquad (25) $$

This quadratic equation admits two independent solutions for $a_2$. The choice of solutions to consider is discussed in the latter.

**Many-body state** Just above, we have fixed all the parameters used in the method proposed by matching on the low-scale expansion (2) and on the high-scale constraint (3). Thus, we are now able to express the many-body state in that new formulation of the MBPT using (18). To be more explicit, $\Gamma_n^{NLO}$ and $\hat{G}_n^{NLO}$ appearing in (19) are given by $\Gamma_1^{NLO} = \gamma_1$, $\hat{G}_1^{NLO} = 1$ and:

$$ \Gamma_2^{NLO} = \left[ 1 + \frac{\gamma_1^2}{a_2 \gamma_2} \right] \gamma_2, $$

$$ \hat{G}_2^{NLO} = \left[ 1 + \frac{\gamma_1^2}{a_2 \gamma_2} \right] \hat{R}_0 \cdot \hat{h}. $$

### 3.2.3 Third order

In this part, we illustrate the truncation method that we propose at the third order ($l = 2$) in perturbation. For that, we use the flexibility to adjust the $\beta$ parameter $a$ *posteriori* on the low-scale expansion and the property that, as discussed before, the energy shift (17) is, in fact, $\beta$-independent. At the second order in perturbation (see above), only one $\beta$ parameter is required to obtain the proper low-scale expansion. But, in the third order, we need two different $\beta$ parameters to match the low-scale expansion up to the third order. Here, the strategy consists of including an additional parameter to (18) as:

$$ \Delta E_{\lambda}(\beta) \rightarrow \frac{\Delta E_{\lambda}(\beta_1) + \Delta E_{\lambda}(\beta_2)}{2}. \qquad (26) $$

We note that this choice is arbitrary in the sense that any weighted arithmetic mean of the form $\Delta E(\beta) \rightarrow r\Delta E(\beta_1) + (1-r)\Delta E(\beta_2)$ can be used, where $0 < r < 1$ is a parameter to adjust. For simplicity and because the strategy developed in this work implies to adjust only three parameters at the third order in perturbation on the low- and high-scale limits, we chose
the arithmetic mean $r = 1/2$. Then, making a truncation of (18) up to $l = 2$, we have the following:

$$ \frac{\langle \Phi_0 | \lambda \hat{h} \cdot \tilde{B}_3(\beta) | \Phi_0 \rangle}{1 + \lambda \gamma_1 / a_3} = \frac{\beta - 1}{\beta^2} \frac{\langle \Phi_0 | \lambda \hat{h} \cdot \hat{R}_0 \cdot \lambda \hat{h} | \Phi_0 \rangle}{(1 + \lambda \gamma_1 / a_3)^2} + \frac{1}{\beta^2} \frac{\langle \Phi_0 | \lambda \hat{h} \cdot \hat{R}_0 \cdot \lambda \hat{w} \cdot \hat{R}_0 \cdot \lambda \hat{h} | \Phi_0 \rangle}{(1 + \lambda \gamma_1 / a_3)^3} \\ = \frac{(\beta - 1) \gamma_2}{\beta^2} \left( \frac{\lambda}{1 + \lambda \gamma_1 / a_3} \right)^2 E_0 + \frac{\gamma_3}{\beta^2} \left( \frac{\lambda}{1 + \lambda \gamma_1 / a_3} \right)^3 E_0, \quad (27) $$

where we have moreover truncated the expression of $\tilde{B}_3(\beta)$ up to the third order in power of $\lambda$ (as in the standard RSPT formulation). Finally, using (26), the N$^2$LO GS energy is given by:

$$ \frac{E_{\lambda}^{N^2LO}}{E_0} = 1 + \gamma_1 \frac{\lambda}{f_{\lambda}^{N^2LO}} + \frac{\gamma_2}{2} \left[ \frac{2\beta_1 - 1}{\beta_1^2} + \frac{2\beta_2 - 1}{\beta_2^2} \right] \left( \frac{\lambda}{f_{\lambda}^{N^2LO}} \right)^2 + \frac{\gamma_3}{2} \left[ \frac{1}{\beta_1^2} + \frac{1}{\beta_2^2} \right] \left( \frac{\lambda}{f_{\lambda}^{N^2LO}} \right)^3 . \quad (28) $$

As before, the parameters $\{\beta_1, \beta_2\}$ are to be determined in terms of $\{\gamma_1, \gamma_2, \gamma_3\}$ and $a_3$ to match the low-scale expansion up to the third order, and then the high-scale parameter $a_3$ can be determined by imposing the high-scale constraint for the GS energy.

### 3.2.4 Avoiding the double counting of correlations

We can observe that our method does not overcount the many-body correlations. This is due to the fact that we impose the low-scale expansion up to a given order in perturbation. Thus, even if many-body correlations are included in an effective Hamiltonian and potentially taken into account within the MBPT framework, the adjustment of the $\beta$ parameters in order to reproduce the low-scale expansion (2) in the limit $\lambda \ll 1$ avoid this eventual double counting. To be more precise, if we reduce our strategy to the RSPT, that is to say setting $\beta = 1$, and keeping an effective Hamiltonian similar to (4), the low-scale expansion (2) cannot be recovered in the limit $\lambda \ll 1$. Therefore, our strategy consists in mimic the infinite series (2) with effective parameters $\{\tilde{\gamma}_n\}$ for which the first parameters match the physical parameters $\{\gamma_n\}$ up to the order $m$ that is considered to truncate the MBPT calculation, i.e. $\tilde{\gamma}_1 = \gamma_1, \dots, \tilde{\gamma}_m = \gamma_m$, and $\tilde{\gamma}_n \neq \gamma_n$ for $n > m$.

## 3.3 Effective field theory interpretation of the extended MBPT

Considering the new method proposed and described above using a correlated effective Hamiltonian, a general and comprehensive Effective Field Theory (EFT) expansion emerges. Indeed, in the GS energy expansion (19) obtained, $\lambda$ plays the role of the low-scale and $f_{\lambda}^{N^1LO}$ the role of the high-scale. This expansion is similar to the low-scale expansion (2) except that the low-scale parameter $\lambda$ is renormalized by the high-scale in-medium function $f_{\lambda}^{N^1LO}$. The proposed method is a systematic procedure to:

* get a better convergence of the result compared to a standard MBPT approach, e.g. finite limit in the strong coupling regime (except when $f_{\lambda}^{N^1LO} = 0$);

* include correlations in an effective way much BMF correlations and even much beyond the correlations accessible using the standard MBPT approach;

* allow the use of the MBPT systematically using an effective correlated Hamiltonian (4);

* reorganize the MBPT expansion to be valid in the low- and high-scale limit;

* remove automatically the double counting of correlations included in the effective Hamiltonian (4) and resulting of the MBPT approach;
* obtain an EFT expansion of the observables in term of $(\lambda/f_{\lambda}^{N^1LO})^n$ where $\lambda$ and $f_{\lambda}^{N^1LO}$ play the role of the low- and high-scales, respectively;

* have analytical, simple, and compact expressions for the observable as a function of well-determined fixed quantities, e.g. {$\gamma_n$} and $\xi_0$, without fitting procedure (eventually the used coefficients can be provided by ab initio calculations or experiments).

But two problems remain:

(i) The polynomial equation similar to (25) admits $l$ different complex solutions for the determination of the parameter $a_{l+1}$ defining the effective Hamiltonian (4).

(ii) The GS has an imaginary non-physical part due to the fact that we allow a complex value for the $a_{l+1}$ parameter.

Before discussing the results for some applications of the extended MBPT method developed and applied to physical systems, we give a simple argument and an explanation to (i) make a pragmatic choice for the parameter $a_{l+1}$ to use and (ii) remove the pathology and recover a real GS energy.

### 3.3.1 Choice of the high-scale solution

The $N^1LO$, the GS energy in the limit $\lambda \to \infty$ is given by a polynomial equation similar to (25). The main purpose of the proposed method is to grasp the high-scale correlations from the HF level. Thus, this motivates the choice of the minimum solution, i.e. set $a_{l+1} = \min\{|x_{l+1}|\}$, where the set $\{x_{l+1}\}$ denotes the $l$ solutions of the high-scale polynomial equation similar to (25), such that the main contribution to the energy in the high-scale limit is given by the leading order<fnref n="4" /> of the expansion (19).

### 3.3.2 Remove the imaginary part

The appearance of an imaginary part in the GS energy comes from the resolution of the polyno- mial equation similar to (25) having $l+1$ distinct solutions. However, the low-scale expansion of the GS energy (19) matches the standard MBPT expansion up to the $l+1$ order and the high-scale constraint (3). Consequently, we have the following properties:

$$ \mathrm{Im} \left( \frac{E_{\lambda}^{N^1LO}}{E_0} \right) = \mathcal{O}(\lambda^{l+2}). $$

This flexibility allows the redefinition of $\Gamma_{l+1}^{N^1LO}$ to get a real GS energy since this change does not affect the systematic procedure requiring only the first $l+1$ order of the low-scale expansion. Thus, this minimal subtraction of the imaginary part consists of keeping the real part of the GS energy.

## 4 Quantum gases context: Application to ultracold atom systems

As in section 2.2, the correspondence for spin-saturated dilute ultracold atom systems with the $N^1LO$ GS energy is given by the high-scale constraint, i.e. the Bertsch parameter $\xi_0 = 0.376$, and the low-scale expansion (2) where $\gamma_1 = 10k_F/9\pi$, $\gamma_2 = 4(11-2\ln 2)k_F^2/21\pi^2$, $\gamma_3 = 0.032k_F^3$, and $\gamma_4 = 0.451k_F^4$ [46]. These results are summarized in table 1.

<footnote n="4">⁴Valid in the case $|a_{l+1}| < 1$, but the same conclusion on this choice occurs to minimize the correction of the next order.</footnote>
<table><caption>Table 1: Low-scale coefficients and high-scale limit (Bertsch parameter), for ($a_s k_F$) > 0 and ($a_s k_F$) < 0, of the ultracold Fermi gas. The dashes mean that the coefficient is not defined.</caption><thead><tr><th></th><th>(a<sub>s</sub>k<sub>F</sub>) &gt; 0</th><th>(a<sub>s</sub>k<sub>F</sub>) &lt; 0</th></tr></thead><tbody><tr><td>&gamma;<sub>1</sub></td><td>k<sub>F</sub>/6&pi;</td><td>10k<sub>F</sub>/9&pi;</td></tr><tr><td>&gamma;<sub>2</sub></td><td>-</td><td>4k<sub>F</sub><sup>2</sup>(11 &minus; 2ln2)/21&pi;<sup>2</sup></td></tr><tr><td>&gamma;<sub>3</sub></td><td>-</td><td>0.032k<sub>F</sub><sup>3</sup></td></tr><tr><td>&xi;<sub>0</sub></td><td>0.376</td><td>0.376</td></tr></tbody></table>

We show in figure 1 the GS energy of the spin-saturated Fermi gas at N'LO, $E^{N'LO}$ given by (19), as a function of ($a_s k_F$). On the BCS side of the crossover ($a_s < 0$), the results obtained with our method go much beyond the standard MBPT approach and reproduce very well the *ab initio* calculations from the weak to the strong coupling regime without substantial improvement when we increase the order of perturbation using the effective Hamiltonian (4). On the Bose-Einstein Condensate (BEC) side of the crossover ($a_s > 0$), due to the appearance of bound dimers in the gas, the GS energy of the system (containing the binding energies of the dimers) takes the following expansion [43,47]:

$$ \frac{E}{E_0} = \frac{1}{6\pi}(a_s k_F) + \frac{32}{375\sqrt{\pi^5/10}}(a_s k_F)^{5/2} + \dots \qquad (29) $$

This expansion does not exhibit a simple polynomial form in ($a_s k_F$) and, thus, the proposed method, in particular the choice of the effective Hamiltonian (4), is not adapted for such a low-scale expansion (except up to the first order). That is why we do not discuss this case further. For reference, the leading order is given in figure 1(c,d) and shows a large error in the non-perturbative regime ($a_s k_F > 1$).

# 5 Nuclear physics context: Application to the Richardson pairing Hamiltonian

We consider a quantum many-body system composed of N doubly degenerate in spin $\sigma = \uparrow, \downarrow$ and equally spaced single-particle levels $n = 1, 2, \dots$. The many-body Hamiltonian expressed in the second quantization formalism, with at most two-body interaction, is given by the so-called Richardson Pairing Hamiltonian<fnref n="5" /> (RPH) [50]:

$$ \begin{aligned} \hat{H}_{\text{RPH}} &= \hat{H}_0 + \hat{V} \\ &= e \sum_{n\sigma} (n-1) \hat{a}_{n\sigma}^{\dagger} \hat{a}_{n\sigma} - \frac{g}{2} \sum_{nm} \hat{a}_{n\uparrow}^{\dagger} \hat{a}_{n\downarrow}^{\dagger} \hat{a}_{m\downarrow} \hat{a}_{m\uparrow}. \end{aligned} \qquad (30) $$

For simplicity and without loss of generality, we consider $e = 1$. The coupling $g$ corresponds to the pairing contribution of the two-body interaction $\hat{V}$. Note that the RPH can be rewritten as (in the total spin $S = 0$ channel only):

$$ \hat{H}_{\text{RPH}} = e \sum_{n\sigma} (n-1) \hat{a}_{n\sigma}^{\dagger} \hat{a}_{n\sigma} - \frac{g}{2} \sum_{nm} \hat{P}_{n}^{+} \hat{P}_{m}^{-}, \qquad (31a) $$

in terms of the pair creation/annihilation operators defined as:

$$ \hat{P}_{n}^{+} = \hat{a}_{n\uparrow}^{\dagger} \hat{a}_{n\downarrow}^{\dagger}, \quad \text{and} \quad \hat{P}_{m}^{-} = \hat{a}_{m\uparrow}^{\dagger} \hat{a}_{m\downarrow}^{\dagger}. \qquad (31b) $$

<footnote n="5">⁵This is a simplified version of the RPH: in a more realistic case, $e$ and $g$ are not constants.</footnote>
<figure><img src="image_3.png" /></figure>

Figure 1: GS energy of the spin-saturated Fermi gas $E_{\lambda}^{N^1LO}$ discussed in section 4 as a function of $(a_s k_F)$ [panels (a) for $a_s < 0$ and (c) for $a_s > 0$] and $(a_s k_F)^{-1}$ [panels (b) for $a_s < 0$, and (d) for $a_s > 0$] obtained with (19) using the dilute Fermi gas correspondence (see section 4) for $l=0$, $l=1$, and $l=2$. $E_0$ corresponds to the exact GS energy of the non-interacting Fermi gas. For reference, using thin gray lines, the first, second, and third orders of MBPT for the dilute Fermi gas [46] are shown, and, the *ab initio* calculations of [48, 49] are represented in each panel by black circles.

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

Figure 2: Configurations (Slater determinants) of the RPH discussed in the text with $N=4$ and $p=4$ ($N_p = 2p$). On the left is displayed the GS and zero-particle-zero-hole excitations ($0p-0h$), the four central diagrams correspond to the $2p-2h$ excitations, and the right diagrams correspond to the $4p-4h$ excitations.

This Hamiltonian commutes with the product of the pair creation and annihilation operators and thus corresponds to a system with no broken pairs, i.e. the RPH links two-particle states in spin-reversed states.

In the following, we consider a system with $N=4$ and no broken pairs (total spin $S=0$) in the GS. We only consider $N_p = 8$ single-particle states, that is, $p=4$ levels. The schematic configurations of this model are shown in figure 2.
<table><caption>Table 2: Energy of the single-particle states for the RPH discussed in the text labeled by the quantum number *n* and the spin-projection on the quantification axis *σ*.</caption><thead><tr><th>n</th><th>σ</th><th>E<sub>nσ</sub></th></tr></thead><tbody><tr><td>1</td><td>±1/2</td><td>-g/2</td></tr><tr><td>2</td><td>±1/2</td><td>e-g/2</td></tr><tr><td>3</td><td>±1/2</td><td>2e</td></tr><tr><td>4</td><td>±1/2</td><td>3e</td></tr></tbody></table>

Using table 2 and figure 2, we deduce that the Hamiltonian has the following matrix form:

$$ \hat{H}_{\text{RPH}} = 2e \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 3 & 0 & 0 & 0 \\ 0 & 0 & 0 & 3 & 0 & 0 \\ 0 & 0 & 0 & 0 & 4 & 0 \\ 0 & 0 & 0 & 0 & 0 & 5 \end{bmatrix} - \frac{g}{2} \begin{bmatrix} 2 & 1 & 1 & 1 & 1 & 0 \\ 1 & 2 & 1 & 1 & 0 & 1 \\ 1 & 1 & 2 & 0 & 1 & 1 \\ 1 & 1 & 0 & 2 & 1 & 1 \\ 1 & 0 & 1 & 1 & 2 & 1 \\ 0 & 1 & 1 & 1 & 1 & 2 \end{bmatrix}. \quad (32) $$

## 5.1 Exact GS energy

In figure 3, we show the GS energy as a function of *g*. In the following, we introduce the reduced GS energy $\bar{E}$ defined as:

$$ \frac{\bar{E}}{E_0} = \frac{E}{E_0} + \frac{3g}{2} \Theta(g), \quad (33) $$

where *E* is the GS energy of the system and *Θ* is the Heaviside step function. This definition of the reduced GS energy allows for a finite limit in the high-scale limit, that is, $\lim_{g \to \pm\infty} \bar{E}/E_0 = 3$, which gives the high-scale constraint $\xi_0 = 3$.

## 5.2 Hartree-Fock and MBPT solutions

### 5.2.1 Hartree-Fock

Starting in the single-particle basis $|i\rangle = |n_i \sigma_i\rangle$ given by table 2 (eigenvectors of $\hat{H}_0$), we define the single-particle HF basis as:

$$ |i\rangle_{HF} = \sum_j C_{ij} |j\rangle, \quad (34) $$

and then the HF equation reads (eigenvalues problem):

$$ \sum_j h_{kj}^{HF} C_{ij} = \epsilon_i^{HF} C_{ik}, \quad (35a) $$

where the HF matrix is defined as:

$$ h_{ij}^{HF} = \langle i | \hat{h}_0 | j \rangle + \sum_{k=1}^{N} \sum_{pq} \langle ip | \hat{v} | jq \rangle, \quad (35b) $$

with $\langle i | \hat{h}_0 | j \rangle = e \delta_{ij}$ and:

$$ \langle ip | \hat{v} | jq \rangle = -\frac{g}{2} (-1)^{\sigma_i + \sigma_j} \delta_{n_i n_p} \delta_{n_j n_q} (1 - \delta_{\sigma_i \sigma_p}) (1 - \delta_{\sigma_i \sigma_q}), \quad (35c) $$
<figure><img src="image_5.png" /><figcaption>Figure 3: GS energy $E_{\lambda}^{N^1LO}$ or reduced GS energy $\bar{E}_{\lambda}^{N^1LO}$ discussed in the text for the RPH as a function of $g$ [panels (a) for $g < 0$ and (c) for $g > 0$] and $g^{-1}$ [panels (b) for $g < 0$, and (d) for $g > 0$] obtained with (19) using the low-scale parameter correspondence of the model for $l = 0, l = 1$, and $l = 2$. $E_0$ corresponds to the reference energy of the GS. For reference, using thin gray lines, the first, second, and third orders of MBPT are shown, and the exact result is represented in each panel by black circles.</figcaption></figure>

is antisymmetrized under the exchange of particles.

We can show that the solution to the HF equation is given by $C_{ij} = \delta_{ij}$ and we deduce the HF GS energy:

$$E_{HF} = \sum_{i=1}^{N} \sum_{j=1}^{N} \left[ \langle i | \hat{h}_0 | j \rangle - \frac{1}{2} \langle ij | \hat{v} | ij \rangle \right] = 2e - g. \quad (36)$$

We note that the HF GS energy is equal to the reference energy obtained considering only the $0p-0h$ excitation, i.e. $E_{HF} = E_{0p-0h}$.

### 5.2.2 MBPT results

Calculations of the many-body perturbation theory provide the low-scale expansion of the GS energy for the pairing model discussed in this section, i.e. the coefficients $\{\gamma_n\}$ appearing in (2) up to fifth order are given by: $\gamma_1 = -1/2$, $\gamma_2 = -7/48$, $\gamma_3 = -1/24$, $\gamma_4 = -77/27648$, $\gamma_5 = 5/864$, etc. For the reduced GS energy (33), only the first order is affected and $\bar{\gamma}_1 = \gamma_1 + 3\Theta(g)/2$ and $\bar{\gamma}_n = \gamma_n$ where $\{\bar{\gamma}_n\}$ are the coefficients of the low-scale expansion (2) for the reduced GS energy $\bar{E}$. These results are summarized in table 3.

In practice, we perform the method proposed in this paper for the reduced GS energy $\bar{E}$ using the low-scale coefficients $\bar{\gamma}_n$ and the high-scale constraint $\xi_0 = 3$, then the energy is obtained simply using (33), i.e. $E^{N^1LO} = \bar{E}^{N^1LO} - 3gE_0\Theta(g)/2$.
<table><caption>Table 3: Low-scale coefficients and high-scale limit, for $g > 0$ and $g < 0$, of the RPH model introduced in the text. In parentheses, we give the values of the regularized coefficients to remove divergences of the GS energy when it is suitable.</caption><thead><tr><th></th><th>$g > 0$</th><th>$g < 0$</th></tr></thead><tbody><tr><td>$\gamma_1$</td><td>$-1/2$ (1)</td><td>$-1/2$</td></tr><tr><td>$\gamma_2$</td><td>$-7/48$</td><td>$-7/48$</td></tr><tr><td>$\gamma_3$</td><td>$-1/24$</td><td>$-1/24$</td></tr><tr><td>$\xi_0$</td><td>$-\infty$ (3)</td><td>3</td></tr></tbody></table>

The final result for the GS energy as a function of the coupling constant $g$ is shown in figure 3. We observe that, even at leading order, the GS energy is reproduced in good approximation in a wide range of interaction strengths. In particular, contrary to the standard MBPT results that diverge for $|g| > 1$ (thin gray lines), our approaches grasp most of the BMF correlations in the strong coupling regime.

## 6 Condensed matter physics context: Application to the Hubbard Hamiltonian

In this section, we apply our strategy developed above for the description of a model that describes low-energy physical properties of strongly correlated fermions navigating into a lattice, namely the Hubbard model.

The Hubbard model [51–53] stands as one of the simplest and most frequently used effective models in theoretical condensed matter physics [54,55]. Its purpose is to capture the general properties of spin-1/2 electrons moving through a lattice by hopping between neighboring sites and subject to a local two-body interaction with strength $U$. In its second-quantized form, the correlated tight-binding Hamiltonian is expressed as:

$$ \hat{H}_{\text{Hubbard}} = -t \sum_{\langle i,j \rangle, \sigma} (\hat{c}_{i\sigma}^{\dagger} \hat{c}_{j\sigma} + \hat{c}_{j\sigma}^{\dagger} \hat{c}_{i\sigma}) + U \sum_{i} \hat{n}_{i\uparrow} \hat{n}_{i\downarrow}, \quad (37) $$

where $t$ is the hopping integral, $\hat{c}_{i\sigma}^{\dagger}$ and $\hat{c}_{i\sigma}$ are, respectively, the creation and annihilation operators of electron with spin $\sigma = \uparrow, \downarrow$ on site $i$, the $\hat{n}_{i\sigma} = \hat{c}_{i\sigma}^{\dagger} \hat{c}_{i\sigma}$ is the occupation number operator, and $U$ is the on-site nearest-neighbor Coulomb repulsion of spin-up and spin-down electrons occupying the same lattice site. We note that in the positive $U$ regime, the model received a considerable renewed interest in two-dimensional (2D) geometry after P. W. Anderson’s proposal in connection to high-$T_c$ superconducting cuprates [56,57]. Furthermore, the one-dimensional (1D) Hubbard model [58] has previously also been proposed as a minimal model to describe the low-energy physical properties of 1D conductors [59] and quasi-1D copper oxides compounds [60–62]. In the following, we will put $t = 1$ to set our unit of energy without losing generality. Note that the summation appearing in the kinetic term in (37) is done on the nearest neighbor sites $j$ for all sites $i$ designed by $\langle i, j \rangle$.

For the sake of simplicity of the discussion and to be able to compare our results with the exact solution, we restrict the model to the so-called half-filled 1D Hubbard chain and to the half-filled 2D four-sites Hubbard Hamiltonian for which exact formulas of the GS energy from exact diagonalization can be found in the literature [63,65].
## 6.1 Application to the 1D Hubbard chain model

As mentioned above, we apply our strategy to the 1D Hubbard chain model described by the Hubbard Hamiltonian (37) at half-filling and considering only nearest-neighbor hopping amplitudes with a positive on-site repulsion ($U > 0$).

The exact GS energy is then obtained using the Bethe anzats and is given in [63]:

$$ \frac{E}{E_0} = \pi \int_0^\infty \frac{d\omega}{\omega} \frac{J_0(\omega)J_1(\omega)}{1 + \exp((U/t)\omega/2)}, \quad (38) $$

where the $J_n$ are the Bessel functions of the first kind. This expression can then be expanded to the low-scale limit $U/t \ll 1$ leading to the low-scale coefficients of expansion (2) and presented in table 4 [64]. Taking the high-scale limit provides the high-scale parameter $\xi_0 = 0$ for $U/t > 0$ and $\xi_0 = 2$ for $U/t < 0$.

The results of the strategy developed in this article and applied in the previous sections, are given in figure 4. Again, we observe a rather good approximation of the GS energy along the wide range of ($U/t$) values. We mention that the strategy is not applicable at the third order because of the vanishing of the third order term in the low-scale expansion ($\gamma_3 = 0$). However, we argue that careful decomposition of the fourth-order MBPT expansion will lead to a substantial improvement. We stick here to the third order to be consistent with the level of approximation employed in the discussion above.

## 6.2 Application to the four-sites 2D Hubbard model

Finally, we follow our strategy for the four-sites 2D Hubbard model (square $2 \times 2$ cluster) that also takes the prerequisites decomposition (37) for which the exact GS energy can be found in [65] and takes the expression:

$$ \frac{E}{E_0} = -\frac{U}{4t} + \frac{1}{2} \sqrt{\frac{16 + (U/t)^2}{3}} \cos\left(\frac{\theta}{3}\right), \quad (39) $$

where:

$$ \theta = \arccos\left(\frac{4U}{t}\left(\frac{3}{16 + (U/t)^2}\right)^{3/2}\right), $$

leading to the low- and high-scale parameters given in table 5. The only subtlety is that, in this case, we have a linear divergence of the GS energy for ($U/t$) < 0 that can be treated as the one appearing in the study of the Richardson Pairing Hamiltonian in section 5. The results are compared to the exact calculations and MBPT results in figure 5 for the first, second, and third order expansion. We observe that the GS energy is reproduced accurately at the third order in perturbation. Also, convergence according to the level of approximation is obtained, as expected by the formal aspect of our expansion.

<table><caption>Table 4: Low-scale coefficients and high-scale limit, for ($U/t$) > 0 and ($U/t$) < 0, of the 1D Hubbard chain model introduced in the text.</caption><thead><tr><th></th><th>(U/t) > 0</th><th>(U/t) < 0</th></tr></thead><tbody><tr><td>γ<sub>1</sub></td><td>-π/16</td><td>-π/16</td></tr><tr><td>γ<sub>2</sub></td><td>7ζ(3)/64π<sup>2</sup></td><td>7ζ(3)/64π<sup>2</sup></td></tr><tr><td>γ<sub>3</sub></td><td>0</td><td>0</td></tr><tr><td>ξ<sub>0</sub></td><td>0</td><td>2</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>Figure 4: GS energy $E_{\lambda}^{N^l LO}$ discussed in the text for the 1D Hubbard chain model as a function of $U/t$ [panels (a) for $U/t < 0$ and (c) for $U/t > 0$] and $(U/t)^{-1}$ [panels (b) for $U/t < 0$ and (d) for $U/t > 0$] obtained with (19) using the low-scale parameter correspondence of the model for $l = 0$ and $l = 1$. $E_0$ corresponds to the reference GS energy. For reference, using thin gray lines, the first and second orders of MBPT are shown, and the exact result is represented in each panel by the black circles.</figcaption></figure>

<table><caption>Table 5: Low-scale coefficients and high-scale limit, for $(U/t) > 0$ and $(U/t) < 0$, of the four-sites 2D Hubbard model introduced in the text. In parentheses, we give the values of the regularized coefficients to remove divergences of the GS energy when it is suitable.</caption><thead><tr><th></th><th>(U/t) > 0</th><th>(U/t) < 0</th></tr></thead><tbody><tr><td>γ<sub>1</sub></td><td>-3/16</td><td>-3/16 (5/16)</td></tr><tr><td>γ<sub>2</sub></td><td>13/512</td><td>13/512</td></tr><tr><td>γ<sub>3</sub></td><td>-3/1024</td><td>-3/1024</td></tr><tr><td>ξ<sub>0</sub></td><td>0</td><td>∞ (0)</td></tr></tbody></table>
<figure><img src="image_7.png" /></figure>

Figure 5: GS energy $E_{\lambda}^{N^1LO}$ or reduced GS energy $\bar{E}_{\lambda}^{N^1LO}$ discussed in the text for the four-sites Hubbard model as a function of $U/t$ [panels (a) for $U/t < 0$, and (c) for $U/t > 0$] and $(U/t)^{-1}$ [panels (b) for $U/t < 0$, and (d) for $U/t > 0$] obtained with (19) using the low-scale parameter correspondence of the model for $l = 0, l = 1$, and $l = 2$. $E_0$ corresponds to the reference GS energy. For reference, using thin gray lines, the first, second, and third orders of the MBPT are shown, and the exact result is represented in each panel by the black circles.

# 7 Conclusion: Discussion and outlook

In this work, we emphasized the possibility of exhibiting free parameters in the MBPT expansion that allow for the use of effective Hamiltonians that contain BMF correlations. In particular, we show that the problem of double counting of BMF correlations, arising from the use of many-body techniques on a renormalized theory, can be automatically solved due to the flexibility of the expansion without fitting procedure or further adjustments. The main advantage of the strategy used in this paper is twofold:

* A systematic convergence of the results, faster than the standard MBPT, is observed at each order of truncation without inducing an increasing computational complexity.

* An explicit parameterization of the GS energy (and of the many-body states) is obtained in terms of the MBPT amplitude and high-scale limit of the bare theory, that is to say, a finite and restricted number of physical constants, e.g. {$\gamma_n$}, accessible from standard many-body theory or analysis of experimental data.

We can mention some restrictions since the applicability of our strategy depends strongly on the state-of-the-art standard MBPT, ab initio methods, and experiments to provide the low-scale constants and the high-scale limit of the bare Hamiltonian.
A compelling aspect resides in the fact that the formalism is independent of the parameterization of the correlated effective Hamiltonian chosen. As shown in this work, our extended MBPT is valid for various correlated systems, from ultracold atoms and finite nuclei, to condensed matter systems described by the Hubbard Hamiltonian. The simple cases for which we tested our strategy have highlighted some restrictions that require further investigations to grasp most of the complexity of realistic models [66–78]. We contend that our approach to characterize quantum many-body problems possesses a degree of universality, as it allows the integration of parameters that can be adjusted to match specific physical properties.

This novel approach could finally be a valuable guide providing new insights and ideas to predict the collective behavior of quantum many-body systems, e.g. via the linear response theory [79–82] that requires the knowledge of the ground state. Another aspect of the new method developed in this paper is the possibility of designing a new parameterization of BMF DFT that extends the domain of applicability for a wide range of systems. This is also a promising way to link DFT and EFT [39, 83]. Finally, due to the link between MBPT and Coupled-Cluster theory, it can be interesting to investigate possibilities to apply our method in such a formalism.

## Acknowledgments

We warmly thank Philippe Cirot, the former director of our institute, who unconditionally and systematically encouraged our work.

**Conflict of interest** The authors declare that they have no conflict of interest.

## References

[1] J. Carlson, S. Gandolfi, F. Pederiva, S. C. Pieper, R. Schiavilla, K. E. Schmidt and R. B. Wiringa, *Quantum Monte Carlo methods for nuclear physics*, Rev. Mod. Phys. **87**, 1067 (2015), doi:10.1103/RevModPhys.87.1067.

[2] T. Pang, *An introduction to quantum Monte Carlo methods*, Morgan & Claypool, Williston, USA, ISBN 9781681741093 (2016), doi:10.1088/978-1-6817-4109-3.

[3] F. Becca and S. Sorella, *Quantum Monte Carlo approaches for correlated systems*, Cambridge University Press, Cambridge, UK, ISBN 9781107129931 (2017), doi:10.1017/9781316417041.

[4] W. M. C. Foulkes, L. Mitas, R. J. Needs and G. Rajagopal, *Quantum Monte Carlo simulations of solids*, Rev. Mod. Phys. **73**, 33 (2001), doi:10.1103/RevModPhys.73.33.

[5] J. E. Lynn, I. Tews, J. Carlson, S. Gandolfi, A. Gezerlis, K. E. Schmidt and A. Schwenk, *Quantum Monte Carlo calculations of light nuclei with local chiral two- and three-nucleon interactions*, Phys. Rev. C **96**, 054007 (2017), doi:10.1103/PhysRevC.96.054007.

[6] J. E. Lynn, I. Tews, S. Gandolfi and A. Lovato, *Quantum Monte Carlo methods in nuclear physics: Recent advances*, Annu. Rev. Nucl. Part. Sci. **69**, 279 (2019), doi:10.1146/annurev-nucl-101918-023600.

[7] U. Schollwöck, *The density-matrix renormalization group*, Rev. Mod. Phys. **77**, 259 (2005), doi:10.1103/RevModPhys.77.259.
[8] T. Papenbrock and D. J. Dean, *Density matrix renormalization group and wavefunction factorization for nuclei*, J. Phys. G: Nucl. Part. Phys. **31**, S1377 (2005), doi:10.1088/0954-3899/31/8/016.

[9] J. Rotureau, N. Michel, W. Nazarewicz, M. Płoszajczak and J. Dukelsky, *Density matrix renormalization group approach for many-body open quantum systems*, Phys. Rev. Lett. **97**, 110603 (2006), doi:10.1103/PhysRevLett.97.110603.

[10] J. Rotureau, N. Michel, W. Nazarewicz, M. Płoszajczak and J. Dukelsky, *Density matrix renormalization group approach to two-fluid open many-fermion systems*, Phys. Rev. C **79**, 014304 (2009), doi:10.1103/PhysRevC.79.014304.

[11] R. J. Bartlett and M. Musiał, *Coupled-cluster theory in quantum chemistry*, Rev. Mod. Phys. **79**, 291 (2007), doi:10.1103/RevModPhys.79.291.

[12] I. Shavitt and R. J. Bartlett, *Many-body methods in chemistry and physics*, Cambridge University Press, Cambridge, UK, ISBN 9780521818322 (2009), doi:10.1017/CBO9780511596834.

[13] R. Roth and J. Langhammer, *Padé-resummed high-order perturbation theory for nuclear structure calculations*, Phys. Lett. B **683**, 272 (2010), doi:10.1016/j.physletb.2009.12.046.

[14] J. Langhammer, R. Roth and C. Stumpf, *Spectra of open-shell nuclei with Padé-resummed degenerate perturbation theory*, Phys. Rev. C **86**, 054315 (2012), doi:10.1103/PhysRevC.86.054315.

[15] A. Tichai, J. Langhammer, S. Binder and R. Roth, *Hartree-Fock many-body perturbation theory for nuclear ground-states*, Phys. Lett. B **756**, 283 (2016), doi:10.1016/j.physletb.2016.03.029.

[16] B. S. Hu, F. R. Xu, Z. H. Sun, J. P. Vary and T. Li, *Ab initio nuclear many-body perturbation calculations in the Hartree-Fock basis*, Phys. Rev. C **94**, 014303 (2016), doi:10.1103/PhysRevC.94.014303.

[17] A. Tichai, R. Roth and T. Duguet, *Many-body perturbation theories for finite nuclei*, Front. Phys. **8**, 00164 (2020), doi:10.3389/fphy.2020.00164.

[18] K. Kowalski, D. J. Dean, M. Hjorth-Jensen, T. Papenbrock and P. Piecuch, *Coupled cluster calculations of ground and excited states of nuclei*, Phys. Rev. Lett. **92**, 132501 (2004), doi:10.1103/PhysRevLett.92.132501.

[19] G. Hagen, T. Papenbrock, D. J. Dean and M. Hjorth-Jensen, *Ab initio coupled-cluster approach to nuclear structure with modern nucleon-nucleon interactions*, Phys. Rev. C **82**, 034330 (2010), doi:10.1103/PhysRevC.82.034330.

[20] T. M. Henderson, J. Dukelsky, G. E. Scuseria, A. Signoracci and T. Duguet, *Quasiparticle coupled cluster theory for pairing interactions*, Phys. Rev. C **89**, 054305 (2014), doi:10.1103/PhysRevC.89.054305.

[21] G. R. Jansen, J. Engel, G. Hagen, P. Navratil and A. Signoracci, *Ab initio coupled-cluster effective interactions for the shell model: Application to neutron-rich oxygen and carbon isotopes*, Phys. Rev. Lett. **113**, 142502 (2014), doi:10.1103/PhysRevLett.113.142502.
[22] G. Hagen, T. Papenbrock, M. Hjorth-Jensen and D. J. Dean, *Coupled-cluster computations of atomic nuclei*, Rep. Prog. Phys. **77**, 096302 (2014), doi:10.1088/0034-4885/77/9/096302.

[23] J. A. Barker, *A quantum-statistical Monte Carlo method; Path integrals with boundary conditions*, J. Chem. Phys. **70**, 2914 (1979), doi:10.1063/1.437829.

[24] J. W. Negele and H. Orland, *Quantum many-particle systems*, Addison-Wesley, Boston, USA, ISBN 9780201125931 (1988).

[25] J. Shumway and D. M. Ceperley, *Path integral Monte Carlo simulations for fermion systems: Pairing in the electron-hole plasma*, J. Phys. IV Fr. **10**, Pr5 (2000), doi:10.1051/jp4:2000501.

[26] C. Cazorla and J. Boronat, *Simulation and understanding of atomic and molecular quantum crystals*, Rev. Mod. Phys. **89**, 035003 (2017), doi:10.1103/RevModPhys.89.035003.

[27] S. R. White, *Density-matrix algorithms for quantum renormalization groups*, Phys. Rev. B **48**, 10345 (1993), doi:10.1103/PhysRevB.48.10345.

[28] L. Liu, H. Yao, E. Berg, S. R. White and S. A. Kivelson, *Phases of the infinite U Hubbard model on square lattices*, Phys. Rev. Lett. **108**, 126406 (2012), doi:10.1103/PhysRevLett.108.126406.

[29] B. Farid, *A note on the many-body perturbation theory*, Philos. Mag. Lett. **79**, 581 (1999), doi:10.1080/095008399176968.

[30] V. I. Yukalov and E. P. Yukalova, *From asymptotic series to self-similar approximants*, Physics **3**, 829 (2021), doi:10.3390/physics3040053.

[31] E. Engel and R. M. Dreizler, *Density functional theory: An advanced course*, Springer, Berlin, Heidelberg, Germany, ISBN 9783642140891 (2011), doi:10.1007/978-3-642-14090-7.

[32] R. M. Dreizler and E. K. U. Gross, *Density functional theory: An approach to the quantum many-body problem*, Springer, Berlin, Heidelberg, Germany, ISBN 9783642861079 (1990), doi:10.1007/978-3-642-86105-5.

[33] N. Schunck, *Energy density functional methods for atomic nuclei*, IOP Publishing, Bristol, UK, ISBN 9780750314237 (2019), doi:10.1088/2053-2563/aae0ed.

[34] G. Colò, *Nuclear density functional theory*, Adv. Phys.: X **5**, 1740061 (2020), doi:10.1080/23746149.2020.1740061.

[35] M. Grasso, *Effective density functionals beyond mean field*, Prog. Part. Nucl. Phys. **106**, 256 (2019), doi:10.1016/j.ppnp.2019.02.002.

[36] K. Moghrabi, M. Grasso, X. Roca-Maza and G. Colò, *Second-order equation of state with the full Skyrme interaction: Toward new effective interactions for beyond-mean-field models*, Phys. Rev. C **85**, 044323 (2012), doi:10.1103/PhysRevC.85.044323.

[37] C.-J. Yang, M. Grasso and D. Lacroix, *Toward a systematic strategy for defining power counting in the construction of the energy density functional*, Phys. Rev. C **96**, 034318 (2017), doi:10.1103/PhysRevC.96.034318.
[38] D. Lacroix, A. Boulet, M. Grasso and C.-J. Yang, *From bare interactions, low-energy constants, and unitary gas to nuclear density functionals without free parameters: Application to neutron matter*, Phys. Rev. C **95**, 054306 (2017), doi:10.1103/PhysRevC.95.054306.

[39] R. J. Furnstahl, *EFT for DFT*, in *Renormalization group and effective field theory approaches to many-body systems*, Springer, Berlin, Heidelberg, Germany, ISBN 9783642273193 (2012), doi:10.1007/978-3-642-27320-9_3.

[40] A. L. Fetter and J. D. Walecka, *Quantum theory of many-particle systems*, Dover Publications, Mineola, USA, ISBN 9780486428277 (2003).

[41] W. Zwerger, *The BCS-BEC crossover and the unitary Fermi gas*, Springer, Berlin, Heidelberg, Germany, ISBN 9783642219771 (2012), doi:10.1007/978-3-642-21978-8.

[42] M. J. H. Ku, A. T. Sommer, L. W. Cheuk and M. W. Zwierlein, *Revealing the superfluid lambda transition in the universal thermodynamics of a unitary Fermi gas*, Science **335**, 563 (2012), doi:10.1126/science.1214987.

[43] S. K. Adhikari, *Nonlinear Schrödinger equation for a superfluid Fermi gas in the BCS-BEC crossover*, Phys. Rev. A **77**, 045602 (2008), doi:10.1103/PhysRevA.77.045602.

[44] D. Lacroix, *Density-functional theory for resonantly interacting fermions with effective range and neutron matter*, Phys. Rev. A **94**, 043614 (2016), doi:10.1103/PhysRevA.94.043614.

[45] A. Szabo and N. S. Ostlund, *Modern quantum chemistry: Introduction to advanced electronic structure theory*, Dover Publications, Mineola, USA, ISBN 9780486691862 (1996).

[46] C. Wellenhofer, C. Drischler and A. Schwenk, *Dilute Fermi gas at fourth order in effective field theory*, Phys. Lett. B **802**, 135247 (2020), doi:10.1016/j.physletb.2020.135247.

[47] D. S. Petrov, C. Salomon and G. V. Shlyapnikov, *Weakly bound dimers of fermionic atoms*, Phys. Rev. Lett. **93**, 090404 (2004), doi:10.1103/PhysRevLett.93.090404.

[48] J. Carlson, S. Gandolfi and A. Gezerlis, *Quantum Monte Carlo approaches to nuclear and atomic physics*, Prog. Theor. Exp. Phys. **01**, 01A209 (2012), doi:10.1093/ptep/pts031.

[49] G. E. Astrakharchik, J. Boronat, J. Casulleras and a. S. Giorgini, *Equation of state of a Fermi gas in the BEC-BCS crossover: A quantum Monte Carlo study*, Phys. Rev. Lett. **93**, 200404 (2004), doi:10.1103/PhysRevLett.93.200404.

[50] R. W. Richardson and N. Sherman, *Exact eigenstates of the pairing-force Hamiltonian*, Nucl. Phys. **52**, 221 (1964), doi:10.1016/0029-5582(64)90687-X.

[51] J. Hubbard, *Electron correlations in narrow energy bands*, Proc. R. Soc. Lond. A: Math. Phys. Eng. Sci. **276**, 238 (1963), doi:10.1098/rspa.1963.0204.

[52] M. C. Gutzwiller, *Effect of correlation on the ferromagnetism of transition metals*, Phys. Rev. Lett. **10**, 159 (1963), doi:10.1103/PhysRevLett.10.159.

[53] J. Kanamori, *Electron correlation and ferromagnetism of transition metals*, Prog. Theor. Phys. **30**, 275 (1963), doi:10.1143/PTP30.275.

[54] M. Qin, T. Schäfer, S. Andergassen, P. Corboz and E. Gull, *The Hubbard model: A computational perspective*, Annu. Rev. Condens. Matter Phys. **13**, 275 (2022), doi:10.1146/annurev-conmatphys-090921-033948.
[55] D. P. Arovas, E. Berg, S. A. Kivelson and S. Raghu, *The Hubbard model*, Annu. Rev. Condens. Matter Phys. **13**, 239 (2022), doi:10.1146/annurev-conmatphys-031620-102024.

[56] P. W. Anderson, *The resonating valence bond state in La$_2$CuO$_4$ and superconductivity*, Science **235**, 1196 (1987), doi:10.1126/science.235.4793.1196.

[57] E. Dagotto, *Correlated electrons in high-temperature superconductors*, Rev. Mod. Phys. **66**, 763 (1994), doi:10.1103/RevModPhys.66.763.

[58] E. H. Lieb and F. Y. Wu, *The one-dimensional Hubbard model: A reminiscence*, Phys. A: Stat. Mech. Appl. **321**, 1 (2003), doi:10.1016/S0378-4371(02)01785-5.

[59] J. M. P. Carmelo, N. M. R. Peres and P. D. Sacramento, *Finite-frequency optical absorption in 1D conductors and Mott-Hubbard insulators*, Phys. Rev. Lett. **84**, 4673 (2000), doi:10.1103/PhysRevLett.84.4673.

[60] R. Neudert et al., *Four-band extended Hubbard Hamiltonian for the one-dimensional cuprate Sr$_2$CuO$_3$: Distribution of oxygen holes and its relation to strong intersite Coulomb interaction*, Phys. Rev. B **62**, 10752 (2000), doi:10.1103/PhysRevB.62.10752.

[61] Z. Chen et al., *Anomalously strong near-neighbor attraction in doped 1D cuprate chains*, Science **373**, 1235 (2021), doi:10.1126/science.abf5174.

[62] H.-X. Wang, Y.-M. Wu, Y.-F. Jiang and H. Yao, *Spectral properties of a one-dimensional extended Hubbard model from bosonization and time-dependent variational principle: Applications to one-dimensional cuprates*, Phys. Rev. B **109**, 045102 (2024), doi:10.1103/PhysRevB.109.045102.

[63] E. H. Lieb and F. Y. Wu, *Absence of Mott transition in an exact solution of the short-range, one-band model in one dimension*, Phys. Rev. Lett. **20**, 1445 (1968), doi:10.1103/PhysRevLett.20.1445.

[64] W. Metzner and D. Vollhardt, *Ground-state energy of the d = 1,2,3 dimensional Hubbard model in the weak-coupling limit*, Phys. Rev. B **39**, 4462 (1989), doi:10.1103/PhysRevB.39.4462.

[65] A. Leprévost, O. Juillet and R. Frésard, *Exact ground state of strongly correlated electron systems from symmetry-entangled wave-functions*, Ann. Phys. **526**, 430 (2014), doi:10.1002/andp.201400129.

[66] M. Xu, L. H. Kendrick, A. Kale, Y. Gang, G. Ji, R. T. Scalettar, M. Lebrat and M. Greiner, *Frustration- and doping-induced magnetism in a Fermi-Hubbard simulator*, Nature **620**, 971 (2023), doi:10.1038/s41586-023-06280-5.

[67] B. Windt, M. Bello, E. Demler and J. I. Cirac, *Fermionic matter-wave quantum optics with cold-atom impurity models*, Phys. Rev. A **109**, 023306 (2024), doi:10.1103/PhysRevA.109.023306.

[68] F. Mivehvar, F. Piazza, T. Donner and H. Ritsch, *Cavity QED with quantum gases: New paradigms in many-body physics*, Adv. Phys. **70**, 1 (2021), doi:10.1080/00018732.2021.1969727.

[69] L. Coraggio, A. Covello, A. Gargano, N. Itaco and T. T. S. Kuo, *Shell-model calculations and realistic effective interactions*, Prog. Part. Nucl. Phys. **62**, 135 (2009), doi:10.1016/j.ppnp.2008.06.001.
[70] L. Coraggio, A. Covello, A. Gargano, N. Itaco and T. T. S. Kuo, *Effective shell-model Hamiltonians from realistic nucleon-nucleon potentials within a perturbative approach*, Ann. Phys. **327**, 2125 (2012), doi:10.1016/j.aop.2012.04.013.

[71] O. Juillet and R. Frésard, *Exotic spin, charge and pairing correlations of the two-dimensional doped Hubbard model: A symmetry-entangled mean-field approach*, Phys. Rev. B **87**, 115136 (2013), doi:10.1103/PhysRevB.87.115136.

[72] J.-Y. P. Delannoy, M. J. P. Gingras, P. C. W. Holdsworth and A.-M. S. Tremblay, *Low-energy theory of the $t-t'-t''-U$ Hubbard model at half-filling: Interaction strengths in cuprate superconductors and an effective spin-only description of La$_2$CuO$_4$*, Phys. Rev. B **79**, 235130 (2009), doi:10.1103/PhysRevB.79.235130.

[73] J. Ripoche, D. Lacroix, D. Gambacurta, J.-P. Ebran and T. Duguet, *Combining symmetry breaking and restoration with configuration interaction: A highly accurate many-body scheme applied to the pairing Hamiltonian*, Phys. Rev. C **95**, 014326 (2017), doi:10.1103/PhysRevC.95.014326.

[74] H.-C. Jiang and T. P. Devereaux, *Superconductivity in the doped Hubbard model and its interplay with next-nearest hopping $t'$*, Science **365**, 1424 (2019), doi:10.1126/science.aal5304.

[75] R. Photopoulos and R. Frésard, *A 3D tight-binding model for La-based cuprate superconductors*, Ann. Phys. **531**, 1900177 (2019), doi:10.1002/andp.201900177.

[76] V. Marino, F. Becca and L. Fausto Tocchio, *Stripes in the extended $t-t'$ Hubbard model: A variational Monte Carlo analysis*, SciPost Phys. **12**, 180 (2022), doi:10.21468/SciPostPhys.12.6.180.

[77] M. Kitatani, L. Si, P. Worm, J. M. Tomczak, R. Arita and K. Held, *Optimizing superconductivity: From cuprates via nickelates to palladates*, Phys. Rev. Lett. **130**, 166002 (2023), doi:10.1103/PhysRevLett.130.166002.

[78] H. Xu, C.-M. Chung, M. Qin, U. Schollwöck, S. R. White and S. Zhang, *Coexistence of superconductivity with partially filled stripes in the Hubbard model*, Science **384**, eadh7691 (2024), doi:10.1126/science.adh7691.

[79] D. Pines and P. Nozières, *The theory of quantum liquids*, W. A. Benjamin, Reading, USA, (1966).

[80] C. Kittel, *Introduction to solid state physics*, Wiley, Hoboken, USA, ISBN 9780471415268 (2004).

[81] G. D. Mahan, *Many-particle physics*, Springer, New York, USA, ISBN 9780306463389 (2000), doi:10.1007/978-1-4757-5714-9.

[82] A. Pastore, D. Davesne and J. Navarro, *Linear response of homogeneous nuclear matter with energy density functionals*, Phys. Rep. **563**, 1 (2015), doi:10.1016/j.physrep.2014.11.002.

[83] R. J. Furnstahl, *Turning the nuclear energy density functional method into a proper effective field theory: Reflections*, Eur. Phys. J. A **56**, 85 (2020), doi:10.1140/epja/s10050-020-00095-y.