{'fraction_non_alphanumeric': 0.037194426427501065, 'fraction_numerical': 0.009393661515063408, 'mean_word_length': 4.539147670961348, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Particle-number projected Bogoliubov coupled cluster theory . Application to the pairing Hamiltonian

Background: While coupled cluster theory accurately models weakly correlated quantum systems, it often fails in the presence of strong correlations where the standard mean-field picture is qualitatively incorrect. In many cases, these failures can be largely ameliorated by permitting the mean-field reference to break physical symmetries. Symmetry-broken coupled cluster, e.g. Bogoliubov coupled cluster, theory can indeed provide reasonably accurate energetic predictions, but the broken symmetry can compromise the quality of the resulting wave function and predictions of observables other than the energy.


I. INTRODUCTION
The description of weakly correlated quantum systems is straightforward. One typically starts with a mean-field approximation providing a qualitative but not quantitative characterization of the problem before incorporating many-body correlations via a method of choice. Coupled cluster (CC) theory [1][2][3][4][5] is one of the most popular and powerful ways to do so, as it accurately accounts for these correlation effects at a reasonable computational cost that scales polynomially with system size.
Strongly correlated systems are much more challenging, but mean-field methods are surprisingly resilient. Indeed, the mean-field picture can, by breaking some or even all of the symmetries of the system, correctly describe at least some of the important qualitative physics in such a situation. For example, dealing with the pairing or reduced BCS Hamiltonian via the breaking of U (1) global gauge symmetry associated with particle-number conservation, mean-field theory correctly occupies every single-particle level. The symmetry-breaking mean-field state can then serve as a reference to expand the exact ground-state wave-function and, in many strongly correlated problems, deliver accurate energies at a reasonable price.
However, methods based on a symmetry-breaking mean-field reference state have important shortcomings. They do not account for the entanglement which is a hallmark of strong correlations. Indeed, quantum fluctuations do eventually lift the fictitious degeneracy associated with the broken symmetry given that the latter is in fact only emergent in finite systems [6][7][8]. It is thus mandatory to restore good symmetry quantum numbers, which not only revises the energy (dramatically in certain situations) but also allows the proper handling of transition operators characterized by symmetry selection rules. The shortcomings are diminished but not resolved in high-orders broken-symmetry, e.g. coupled cluster, methods given that they can in fact only restore the broken symmetry in an all-order limit.
At the mean-field level, it has long been realized that the qualitative failures of breaking symmetries are largely resolved by symmetry projection [9][10][11][12][13][14], at least in finite systems. Indeed, the action of the projector retains only that portion of the symmetry-breaking mean-field state that has the correct symmetry properties. It seems natural to extend this idea to the case of symmetry-breaking coupled cluster theory. Due to the complexity of the coupled cluster wave function, there has not been any practical attempt in the past to merge full symmetry projection and coupled cluster theory. This has begun to change in the past years given that Duguet [15,16] and Scuseria [17] have proposed independent ways to achieve such a goal. Tsuchimochi and Ten-no have provided a third approach [18], albeit so far at the linearized coupled cluster level.
The first goal of this paper is to discuss the formal relationship between the approaches introduced in Refs. [15,16] and in Ref. [17]. While these two methods were originally formulated in different ways, it is presently shown that they follow closely related ideas and can be formulated within the same basic framework. We seek to standardize the language and reveal the similarities between both methods, as well as to highlight their differences.
The second goal of the present paper is to test symmetry-projected coupled cluster theory on the reduced BCS or pairing Hamiltonian. In this context, the broken and restored symmetry is U (1) global-gauge symmetry associated with the conservation of particle number. Indeed, the mean-field solutions generated in the strongly correlated regime of the pairing Hamiltonian are particle-number breaking BCS states. Although the results are restricted to the schematic pairing Hamiltonian, the formalism presented in this manuscript is directly applicable to the more general context of Bogoliubov coupled-cluster [19,20] such that the applied method is coined as particle-number projected Bogoliubov coupled cluster (PBCC) theory.
The paper is organized as follows. In Sec. II, the main features of the PBCC formalism are introduced. Additional formal and technical details are supplied in set of appendices. While the presentation remains general in Sec. II, it is further specified to the pairing Hamiltonian in App. E. Numerical results obtained for the pairing Hamiltonian are discussed in Sec. III while conclusions are given in Sec. IV.

II. MANY-BODY FORMALISM
In addition to laying down the formalism in view of applying it, the goal of the present section is to clarify the structural identity of the two versions of the PBCC formalism introduced independently in Refs. [15,16] and [17] on the one hand, and to compare the equations put forward in both cases to solve for the gauge-angledependent cluster amplitudes on the other hand.

A. Bogoliubov coupled cluster
Before coming to PBCC theory, it is first necessary to briefly outline the underlying Bogoliubov coupled cluster (BCC) formalism [19,20]. Bogoliubov coupled cluster begins with the introduction of a Bogoliubov reference state which is a vacuum for the set of quasiparticle operators β and β † obtained from particle ones c and c † via a linear Bogoliubov transformation [10] β where the unitary Bogoliubov transformation W is typically obtained by solving Hartree-Fock-Bogoliubov (HFB) equations but does not need to be. While |Φ reduces to a Hartree-Fock (HF) Slater determinant when V = 0, it is in general not an eigenstate of the particle number operator 1 A, i.e. it breaks U (1) global gauge symmetry.
Many-body correlations are incorporated via the usual exponential ansatz of coupled cluster theory where U n is a n-tuple (i.e. 2n quasiparticle) connected excitation operator with respect to |Φ such that [U p , U q ] = 0. Defining the n-tuple excitation operator 2 allows one to introduce n-tuple excitations of the Bogoliubov vacuum through With those definitions at hand, the insertion of the BCC wave function ansatz into the Schrödinger equation yields the BCC energy and amplitude equations 1 The reader is advised not to be confused between the operator A (math style) and its eigenvalue A (roman style) used throughout the paper. 2 The n-tuple (i.e. 2n quasiparticle) de-excitation operator is con- where intermediate normalization Φ|Ψ BCC = Φ|e U |Φ = 1 was used in Eq. 7a.
These equations are more usually formulated in terms of a similarity-transformed Hamiltonian, but these similarity-transformed equations are equivalent to those presented in the foregoing. The amplitudes U 2n0 µ are obtained by solving Eq. 7b, which constitutes a set of coupled nonlinear equations.

B. Particle number projection
If BCC Eq. 7b is satisfied for all µ, the BCC wave function is an exact eigenfunction of the Hamiltonian. Provided that the Hamiltonian governing the A-body quantum system is a particle-number conserving operator such that [A, H] = 0, elementary considerations dictate that in the absence of degeneracies the eigenstates of the Hamiltonian H carry good symmetry quantum number A, i.e. are particle number eigenfunctions. The BCC wave function is thus an eigenfunction of both H and A.
In practical calculations, however, the expansion of U = n U n must be truncated, and generally only a subset of the BCC amplitude equations can then be satisfied. In this case, the BCC wave function is no longer a particle number eigenfunction. The idea of PBCC is to remedy this by applying a particle number projection operator to the BCC wave function 3 The projection operator on particle number A reads as where the operator R(ϕ) ≡ e iAϕ forms the abelian compact Lie group U (1) ≡ {R(ϕ), ϕ ∈ [0, 2π]} associated with the global rotation of the A-body fermion system in gauge space. Let us note that it is helpful to use a variation after projection (VAP) scheme in which the Bogoliubov transformation W minimizes the projected mean-field energy rather than the straight mean-field one. Results presented in Sec. III B are all based on a Bogoliubov transformation obtained in this way.

C. Ground-state energy
Inserting the PBCC wave-function ansatz into the Schrödinger equation and commuting the Hamiltonian with the projection operator leads to 3 When all BCC amplitude equations are satisfied, the number projection is obviously harmless and |Ψ PBCC = |Ψ BCC .
Expanding P A according to Eq. 9 and introducing offdiagonal, i.e. left-rotated, unexcited norm and Hamiltonian kernels where the gauge-rotated Bogoliubov bra state is defined as Φ(ϕ)| ≡ Φ|R(ϕ), the PBCC energy is written as Expressing Φ(ϕ)| via a Thouless transformation of Φ| (see App. A 1 for details) where the Thouless operator Z(ϕ) is a pure de-excitation operator when acting to the right (and a pure excitation operator when acting to the left), the two kernels can be rewritten as where use was made of the identity e Z(ϕ) |Φ = |Φ . The operators entering Eq. 14 are defined through the similarity transformation which is discussed in details in App. A 2.
The angle-dependent state exp(U Z (ϕ))|Φ can be written in coupled cluster form, as Here, the gauge-angle-dependent operator W (ϕ) ≡ n=0 W n (ϕ) is made out of n-fold excitation operators and includes a normalization constant W 0 (ϕ). They can be obtained by expanding e U Z (ϕ) , normal ordering each term with respect to |Φ and retaining the non-zero terms of e U Z (ϕ) |Φ before rewriting the resulting wave function in coupled cluster form. The excitations W n (ϕ) are denoted as disentangled clusters, because this process above is equivalent to disentangling the algebra formed by excitations, de-excitations and quasiparticle numberconserving operators. Adopting the language of disentangled clusters allows one to rewrite Eq. 14 as where T (ϕ) ≡ W (ϕ) − W 0 (ϕ) and with where C denotes the normal-ordered part of n-tuple excited contributions to e U Z (ϕ) obtained by contracting strings of U k operators. It is to be noted that W n (ϕ) receives contributions from all U k and from all excitation ranks in the expansion of e U Z (ϕ) . Morever, even if U is truncated at some low order of excitation, W will generally be non-zero for all excitation orders. Practical calculations thus require that both the broken-symmetry cluster operator U and the disentangled cluster operator W (ϕ) be truncated, or at least that higher-order cluster operators be approximated in terms of lower-order ones. While Eq. 18 provides a formal definition of the gaugerotated cluster operators W n (ϕ), it is not advocated to access them in this way in the following, essentially because the summation over l in the second line of Eq. 18 converges too slowly to be of use.
Equation 17 demonstrates the typical factorization property of off-diagonal kernels that leads to introducing the correlated part of N (ϕ) and the connected part of H(ϕ) via where the index C stipulates the connected character of the kernel. With all these developments at hand, the final form of the PBCC energy reads as

D. Discussion
Equations 17, 19 and 20 stipulating the way the energy is to be computed in PBCC theory appear identically in Ref. [15,16] and [17]. Both approaches make the same use of the similarity-transformed operator H Z (ϕ), once recognized that the Thouless operator Z(ϕ) corresponds to both the operator V 1 introduced in Ref. [17] and to the operator (R −− (ϕ)) † introduced in Ref. [16]. Furthermore, n(ϕ) = exp(W 0 (ϕ)) and the disentangled clusters appear identically once recognized that T (ϕ) [15,16] is nothing but W (ϕ) − W 0 (ϕ) [17]. Eventually, both formalisms are thus structurally identical. However, they effectively differ in the way the disentangled clusters are actually determined as is illustrated in details in the next section. The latter difference is of importance in a practical, i.e. truncated, implementation of the formalism as is illustrated later via the numerical application.

E. Amplitude equations
Before coming to the disentangled clusters, let us first make a comment on the determination of the gaugeunrotated ones. In the PBCC context, the U n operators may be determined in a projection-after-optimization (PAO) scheme, i.e. they can be obtained in a first step by solving standard BCC equations (Eq. 26), or via a more involved optimization-after-projection (OAP) scheme by solving BCC equations in presence of the projection operator [21], e.g.
The distinction between PAO and OAP is similar in spirit to obtaining the Bogoliubov transformation from a variational problem in presence or in absence of P A as discussed earlier.
Let us now come to the main point concerning the disentangled clusters. As for the gauge-dependent constant W 0 (ϕ), the same first-order ordinary differential equation (ODE) was proposed in Refs. [15,16] and [17] (see Sec. II E 2 below). Contrarily, the equations used to determine the gauge-rotated cluster operator T (ϕ) = W (ϕ) − W 0 (ϕ) are different in both variants of the formalism. While amplitude equations generalizing those at play in BCC, i.e. at ϕ = 0, were put forward in Ref. [15,16], a set of first-order ODEs were introduced in Ref. [17]. In the limit where PBCC becomes exact, the two versions of the formalism are identical. When PBCC is not exact, however, solving angle-dependent BCC amplitude equations for T (ϕ) produces a constant connected hamiltonian kernel h(ϕ), as a consequence of which PBCC implemented in this manner is in fact energetically equivalent to BCC and must be discarded. See App. D for more details on this important point.
Let us now detail the two options to determine the disentangled cluster amplitudes.

Angle-dependent BCC equations
Left-multiplying the Schrödinger equation by the ntuply (i.e. 2n quasiparticle) excited rotated Bogoliubov state one obtains the equations of motion which can be written in a compact form by invoking n-tuply excited norm and Hamiltonian kernels 4 defined throughÑ which reduce to the unexcited kernels introduced in Eq. 14 for B µ = 1.
After several steps of algebraic manipulations detailed in App. B, Eq. 23 can be rewritten as a workable set of equations for the n-tuple gauge-dependent cluster amplitudes. Illustrating this for single and double excitations, Eq. 23 reduces to in which the exponential naturally terminates by virtue of the connected character of the kernel. Equations 25 are formally identical to standard BCC equations except that the Hamiltonian H is replaced by its similaritytransformed partner H Z (ϕ) at each gauge angle ϕ. As a matter of fact, Eq. 25 reduces at ϕ = 0 to plain, i.e. unrotated, BCC equations (Eq. 7b) that are to be used to determine the plain cluster amplitudes U n , e.g.
Determining the angle-dependent cluster amplitudes T 1 (ϕ), T 2 (ϕ), . . . by solving the angle-dependent BCC equations (Eq. 25) was put forward in Refs. [15,16]. Gauge-rotated BCC amplitude equations do not provide access to the norm kernel and must thus be complemented by another equation delivering W 0 (ϕ). To do so, it was proposed in Refs. [15,16] to use a first-order ordinary differential equation (ODE) that belongs in fact to a systematic set of first-order ODEs that can be used to determine all W n (ϕ).

Ordinary differential equations
At gauge angle ϕ = 0, the rotation operator R(ϕ) is the identity. Accordingly U Z (ϕ = 0) = U , and therefore That we know a set of initial conditions suggests that one might obtain W (ϕ) by some form of evolution. The differential equations governing W (ϕ) can be obtained 5 by differentiating both sides of Eq. 16 defining W (ϕ). Writing Eq. 16 as and differentiating it, one finds or, since W (ϕ) is in normal order and commutes with its derivative, It can be shown that 6 where A 02 Z (ϕ) is the de-excitation part of the similaritytransformed number operator Left-multiplying Eq. 30 by all possible Φ|B µ leaves us with a set of coupled ODEs for the amplitudes defining W n (ϕ) for n ≥ 0. For example, this gives up to W 2 (ϕ) The above ODEs stipulate that the dependence of the cluster amplitudes on the gauge angle ϕ is driven by kernels of the particle number operator (i.e. of its deexcitation part) that is nothing but the infinitesimal generator of the U (1) group. From the above examples, it is clear that the derivative of W k (ϕ) involves contributions from W k+1 (ϕ). By truncating at a certain excitation level, these equations can be decoupled, and the low order amplitudes can be obtained from an ODE solver, supplemented by the initial conditions given in Eq. 27. Note however that truncating the set of ODEs (for any given truncation of the cluster operator U ) implies an approximation on the action of the projection operator. The effect of this truncation will be gauged in the numerical application below. In Refs. [15,16], it was advocated to complement Eq. 25 with the first differential equation above to access W 0 (ϕ). In Ref. [17], it was advocated to determine angle-dependent cluster amplitudes W n (ϕ) by solving the set of differential equations.
While the derivation of the ODEs given above is straightforward, it may be useful to provide an alternative perspective providing a deeper insight. Because projected states are eigenstates of the particle number operator A, one may write which can be rearranged to yield 7 where the n-tuply excited particle-number kernel 8 involved is defined as Equation 36 provides a set of exact identities highlighting the effect of the symmetry projection when focusing on the particle number operator for which we know the exact answer, i.e. the number A. Contrarily to the gauge-rotated BCC equations (Eq. 25) deriving from the Schrödinger equation, these identities reflect the structure of the U (1) symmetry group and are independent 7 When the de-excitation operator is trivial, i.e. Bµ = 1, Eq. 34 is nothing but the left projection of the eigenequation for the particle number by the reference bra Φ|. 8 The nature of the operator kernelȂµ(ϕ) differs from the one introduced in Eq. 24b for the Hamiltonian by the relative position of the de-excitation operator Bµ and of the operator of interest, i.e. H or A. Most interestingly, while an equation similar to Eq. 23 could be written for A, the same is not true regarding Eq. 34 for H.
of the actual dynamics generated by the Hamiltonian. While Eq. 25 equates the BCC residuals to zero at each gauge angle ϕ, property 36 spans the entire volume of the symmetry group at once. As a matter of fact, fulfilling Eq. 36 for each µ can be rephrased into satisfying a set of differential equations stipulating how the evolution of the n-tuply excited norm kernel with respect to the gauge angle is driven by the n-tuply excited particlenumber kernel Indeed, satisfying Eq. 38 implies that such that Eq. 36 is fulfilled. It happens that, after several steps of algebraic manipulations detailed in App. C, Eq. 38 actually delivers the set of ODEs examplified in Eq. 33 from which the disentangled cluster amplitudes W n (ϕ), n ≥ 0, are determined. This ends demonstrating that the set of ODEs is nothing but a rephrasing of Eq. 34, which itself translates the restoration of good particle number symmetry at each n-tuply excitation level.

III. APPLICATION
We apply the PBCC formalism laid out in Sec. II to the pairing Hamiltonian problem at the single and double excitation levels (i.e. U = U 1 + U 2 ). In this context, the general Bogoliubov reference state simplifies into a simpler BCS vacuum. Consequently, the presently applied method is termed PBCS-CCSD, or PBCS-CCD in case one only includes doubles. Our main objective is to gauge the merits of the symmetry projection, which will be done by comparing PBCS-CCSD results to exact ones [22,23] and to those obtained in Ref. [19] via the BCS-CC(S)D method from which the particle-number projection is absent. The reader interested in comparing the results displayed below with those obtained with recently proposed alternatives of similar quality is referred to Refs. [24,25].
Results are only provided for the ODE-based approach since, as mentioned earlier and as discussed 9 in details in App. D, solving angle-dependent BCC equations for T (ϕ) makes PBCC energetically equivalent to BCC and must thus be discarded. In addition to the truncation of the cluster operator U to singles and doubles, the set of ODEs is truncated at the doubles level, i.e. they are solved for W 0 (ϕ), W 1 (ϕ) and W 2 (ϕ) while assuming that W n (ϕ) = 0 for n ≥ 3.

A. Pairing Hamiltonian
The pairing or reduced BCS Hamiltonian is defined as Here, λ is the chemical potential, p denotes singleparticle energy levels, and G is the interaction strength.
The pair operators can be mapped to fermions using and satisfy an SU (2) algebra Single-particle states p andp are said to be pair conjugated and are degenerate such that p = p . For simplicity, the single-particle levels are taken to be equally spaced, so that p = p ∆ where ∆ is the level spacing. Because the interaction is attractive, the mean-field solution spontaneously breaks U (1) global-gauge symmetry for all G greater than some critical value denoted as G c and takes the form of a BCS reference state. For G < G c , the mean-field solution reduces to the symmetry conserving HF Slater determinant.
Having introduced the pairing Hamiltonian, the explicit expression of all the quantities necessary to implement PBCS-CCSD theory are provided in App. E. In the weakly correlated regime (G G c ) the symmetry adapted coupled cluster (i.e. coupled cluster based upon a Hartree-Fock reference, here called HF-CCSD though the effects of singles are zero by symmetry) is fairly accurate. It is entirely equivalent to BCS-CCSD in this regime given that the BCS reference state reduces to the HF Slater determinant. Increasing G, HF-CCSD fails spectacularly when approaching G c . Contrarily, BCS-CCSD continues to deliver reasonably accurate energies, albeit with a discontinuous derivative at G = G c reflecting the change in character of its mean-field reference. Number projected BCS (PBCS), happens to be wonderfully accurate for larger G but, while exact at G = 0, does not perform as well for small G. Its error is largest for G slightly larger than G c . It is just a remarkable feature that PBCS is better for large G and coupled cluster for small G as well as that PBCS is better for smaller numbers of levels and coupled cluster for larger as the energetic benefit of number projection dwindles as the number of correlated particles increases.

Energetics
While both PBCS and BCS-CCSD perform well for some values of G and not as well for others, PBCS-CC gives almost exact energies for all G, even the region near G c where neither PBCS nor BCS-CCSD performs particularly well. As anticipated from the well-celebrated PBCS theory, the strongest impact of the projection is obtained in the regime where the symmetry is weakly broken in the first place, i.e. after the phase transition for intermediate coupling strengthes. Eventually, merging symmetry projection and coupled cluster theory leads thus to a wave function that is significantly better than either of its components.
The fractional correlation picture can be misleading in certain cases. The right panel shows that for small G, all methods are close to the right answer. For larger G, the error in PBCS decreases (indeed, PBCS is an eigenfunction of the Hamiltonian in the G → ∞ limit) while the total errors of BCS and BCS-CCD are still increasing at G = 3 G c . Adding number projection to BCS-CCD yields an absolute error in the total energy too small to be distinguished on the plot.
The benefits of symmetry projection are most pronounced for smaller systems, while as Ref. [19] shows, BCS-CCD seems to be on average more accurate for larger system sizes. It is thus worth checking how PBCS-CCD performs as a function of system size. This is shown in Fig. 2 where the fraction of correlation energy recovered at half filling is displayed as a function of the number of levels for two different values of G/G c . In keeping with our expectations, PBCS tends to perform worse as system size increases, while BCS-CCD does better. While the quality of PBCS-CCD does deteriorate ever so slightly with increasing system size, the method delivers nearly exact energies across the board.
The situation away from half filling is broadly similar, as can be seen in Fig. 3 where the fraction of correlation energy recovered for 100 levels as a function of the filling fraction is displayed. The interaction strength G = 1.5 G c is employed, where G c denotes here the critical G at half filling. This value of G is sufficiently large that number symmetry is broken for all filling fractions under consideration. BCS-CCD is most accurate near half filling, while PBCS is in contrast most accurate for small or large filling fractions. PBCS-CCD seems to combine these two effects and is exceptionally accurate for all filling fractions. On the other hand, the strongest benefit from the projection is obtained for filling fractions close to zero or one, i.e. where the symmetry breaking is the weakiest 10 . Note that single excitations have been presently excluded from BCS-CC and PBCS-CC calculations simply because including them is more cumbersome computationally. Because single excitations have the effect of adjusting the reference determinant, BCS-CCSD or PBCS-CCSD results can only be expected to be better than those displayed in Fig 3.

One-body observables
Thus far, the focus has been on the energy. It is important to check the quality of the predictions for other properties as well. Generally, properties can be computed by differentiating the total energy with respect to some parameter in the Hamiltonian. For example, single-particle occupation probabilities n p can be expressed as dE/d p . From them, any one-body ground-state quantity can be obtained.
The errors on single-particle occupation probabilities are shown in Fig. 4. As one might expect BCS displays noticeable errors that are to some extent mitigated by PBCS. It was already noted in Ref. 19 that BCS-CCD gives surprisingly accurate occupancies. Adding number 10 In the nuclear physics context, this typically corresponds to semimagic nuclei near shell closures for which we expect PBCC to be the most superior to BCC. This however remains to be seen in realistic nuclear applications. projection yields essentially no error such that PBCS-CCSD delivers basically exact results for any one-body property.

Particle-number restoration
In addition to the great performance of PBCS-CC for ground-state energy and one-body observables, the actual restoration of good particle number at each BCC truncation order must be proven. To do so, Fig. 5 displays the dispersion of particle number with and without projection as well as with and without BCC corrections, i.e. at the (P)BCS and (P)BCS-CCSD levels.
As expected, one first observes that the particlenumber dispersion in zero in all cases for G < G c given that the methods are symmetry-conserving to begin with over this interval. Beyond G c , the dispersion grows monotonically in the BCS calculation, reaching large values in the strongly correlated regime. While including dynamical correlations at the BCS-CCSD level does effectively reduce the particle-number dispersion while improving the energetics, the dispersion remains significant 11 . On the other hand, including the projection on top of BCS does restore the symmetry exactly for all values of G c via the inclusion of non-dynamical correlations but misses dynamical correlations that are of im- 11 It remains of interest to test in which fashion the particle number is effectively restored as one includes higher cluster amplitudes toward the exact BCC wave function. portance for high-precision energetics as illustrated previously. Eventually, going to PBCS-CCSD brings the best of both worlds, i.e. both dynamical and non-dynamical correlations are captured in a consistent way such that the particle-number symmetry is indeed restored for all values of G c while obtaining nearly perfect energetics.
While being convincingly close to zero for any practical purpose, one observes that the dispersion is not strictly null near G c in the PBCS-CCSD calculation. This relates to the fact that the set of ODEs is truncated (i.e. W 0 (ϕ), W 1 (ϕ) and W 2 (ϕ) are presently retained), which implies an approximation on the action of the projection operator. As a matter of fact, the present result demonstrates that the truncation of the ODEs at the double level does provide a faithful account of the projection operator for all values of the coupling strength. If needed, one can anyway envision to include W 3 (ϕ) to reach very high precision. PBCS-CCD results are displayed in Fig. 6 while varying the chemical potential around its BCS value. The quantity reported is the grand potential H − λ(A − A) , where A = 100 denotes the targeted number of particles. Clearly, BCS and BCS-CCD produce essentially identical deviations, which means that the BCS-CCD correlation energy depends only weakly on the chemical potential.
Projecting the corresponding BCS states, the PBCS energy depends (mildly) on the chemical potential, i.e. on the average particle number defining the underlying BCS state. Doing the same on top of the BCS-CCD states, the PBCS-CCD energy is independent of the underlying chemical potential. This demonstrate that the consistent combination of dynamical correlations through the inclusion of singles and doubles and of non-dynamical correlations via the projection erase the memory of the underlying BCS state at rather low CC order.

IV. CONCLUSIONS
Coupled cluster theory very efficiently captures weak correlations but, when based upon a symmetry-adapted reference determinant, tends to break down when the mean-field solution is strongly unstable toward symmetry breaking. While this failure can be avoided by using a symmetry-breaking reference state, one pays the price in lost quantum numbers and potential inaccuracies in properties sensitive to the broken symmetry. The choice between accurate energetics on the one-hand and good symmetries on the other is resolved at the mean-field level by symmmetry projection, and the same appears to be true at the correlated level. While symmetry-projected mean-field methods do not typically describe weak correlations and symmetry-adapted coupled cluster methods do not generally describe strong correlations, the combi-nation of symmetry projection and coupled cluster theory yields a wave function that can accurately describe both regimes. Indeed, for the pairing Hamiltonian PBCS-CC is energetically exact for both small and large G and even gives essentially exact results near G c where neither PBCS nor BCS-CC alone is adequate. Moreover, the wave function of PBCS-CC appears to yield essentially exact one-body properties. All of this can be obtained with a computational cost not too much larger than that of the underlying BCS-CC. While we have not yet applied our techniques to more physical Hamiltonians, these early applications suggest that symmetry-projected coupled cluster theory has the potential to be an important part of the computational arsenal.
The gauge-rotated Bogoliubov state can be expressed via a non-unitary Thouless transformation of the unrotated one according to where the one-body Thouless operator solely contains a pure de-excitation part. The skewsymmetric Thouless matrix is expressed via the transition Bogoliubov transformation made out of the two matrices where (U, V ) define the unrotated Bogoliubov vacuum |Φ .

Similarity transformation
Similarity-transformed quasiparticle operators are defined through

Given the elementary commutators
one finds which can be compacted under the form of a non-unitary linear gauge-dependent Bogoliubov transformation One notes that a pure de-excitation operator is invariant under such a transformation, i.e.
Fermionic anti-commutation rules of the transformed quasiparticle operators are maintained given that Given a normal-ordered operator defined through its similarity-transformed partner reads as 13 13 It is important to understand the notation used in the present document, i.e. W (ij) Z (ϕ) denotes the transformed operator of W ij such that the notation (ij) is a sole reminder of the normalordered nature of the original operator but does not characterize the normal-ordered nature of the transformed operator. On the other hand, W mn Z (ϕ) does denote the normal-ordered part of the transformed operator W Z (ϕ) containing m (n) quasiparticle The last line defines a sum of normal-ordered terms obtained by applying Wick's theorem with respect to |Φ to creation (anihilation) operators. Consistently, one may use the notation W mn(ij) Z (ϕ) to further represent the part of W mn Z (ϕ) originating from W ij ; see Ref [16] for more details. the transformed operator. Given the nature of the transformation defined in Eq. A8, the resulting terms have at least as many annihilation operators as the original operator (j) and possibly up to the total number of original quasiparticle operators (i + j). The number of creation operators ranges from 0 to the original one (i) such that the overall number of quasiparticle operators is bound to remain between j and i + j in each term.

Appendix B: Angle-dependent BCC equations
Let us derive the first set of equations providing access to all amplitudes but W 0 (ϕ), i.e. to all the T n (ϕ). One starts from Eq. 23 obtained by left-multiplying the Schrödinger equation with n-tuply excited rotated Bogoliubov states and written in a compact form as is obviously valid for B µ = 1. Inserting Eq. 19b into it, one obtains the connected form of the energy eigenequation where the connected Hamiltonian kernel, denoted as h(ϕ) in Eq. 19b, naturally terminates such that the above equation formally reads as the standard, i.e., diagonal or unrotated, BCC energy equation [20]. These connected residuals can be trivially expressed in terms of the double similarity transformed Hamiltonian according to

Singly-excited equation
Following the same steps as in Sec. II C, singly-excited kernels introduced in Eq. 24 can be rewritten as in terms of connected singly-excited kernels defined through where the connected character stipulates that the cluster operator is necessarily contracted with the operator H Z (ϕ). Inserting Eq. B5 into Eq. B1 and exploiting Eq. B2, one obtains the connected single equation under the form where h k1k2 (ϕ) naturally terminates such that Eq. B7 formally reads exactly as the standard, i.e., diagonal or unrotated, single BCC amplitude equation [20].

Doubly-excited equation
To unreveil the hierarchy of angle-dependent BCC equations, the doubly-excited one is now worked out. Following the same steps as before, doubly-excited kernels introduced in Eq. 24 can be rewritten as where connected doubly-excited kernels are defined as Inserting Eq. B8 into Eq. B1 and exploiting Eqs. B2 and B7, one obtains the connected double equation under the form which naturally terminates such that Eq. B10 formally reads exactly as the standard, i.e., diagonal or unrotated, double BCC amplitude equation [20].

n-tuply excited equations
As for connected singly-and doubly-excited equations, the derivation of the n-tuple equation invokes all equations of lower rank. Reasoning by recurrence, one can prove that where the naturally terminating connected n-tuply excited kernel of the Hamiltonian is defined as

Appendix C: Differential equations
We derive now the second set of equations providing access to all angle-dependent amplitudes W n (ϕ), i.e. to W 0 (ϕ) and T n (ϕ) for n ≥ 1.

Master equation
Taking the derivative of the n-tuply excited norm kernel with respect to the gauge angle gives 2. Particle number kernel The particle number operator is defined in an arbitrary single-particle basis through with a 11 pq = δ pq . In the quasiparticle basis of |Φ , the normal-ordered form of A reads as where the matrix elements of the various normal-ordered parts are given in terms of the Bogoliubov transformation defining |Φ via In agreement with Eq. A12, the similarity-transformed particle-number operator is written as where the expression of each angle-dependent matrix element can be found in Ref [16]. With the above definitions at hand, the particlenumber kernel introduced in Eq. 37 can be re-expressed as where the kernelÃ 02 µ (ϕ) is defined similarly to the Hamiltonian kernel introduced in Eq. 24.

Unexcited amplitude equation
Applying Eqs. C1-C6 for B µ = 1, the unexcited norm kernel satisfies Introducing the connected unexcited kernel of A 02 and further noticing that Eq. C7 can be rewritten as a differential equation for the unexcited norm kernel Given that N (ϕ) = Φ(ϕ)|Φ n(ϕ) and that one obtains the differential equation for the correlated part of the norm kernel d dϕ n(ϕ) = i a 02 (ϕ)n(ϕ) .

Singly-excited amplitude equation
Applying Eqs. C1-C6 for B µ = B k1k2 , the singlyexcited norm kernel satisfies where the connected singly-excited kernel of A 02 reads as and where it was exploited that Using Eqs. B5a and C10, along with the fact that W 1 (ϕ) is eventually shown to satisfy the first-order differential equation which connects it to W 2 (ϕ) through Eq. C15.

Doubly-excited amplitude equation
Applying Eqs. C1-C6 for B µ = B k1k2k3k4 , the doublyexcited norm kernel satisfies Using Eq. B8a and exploiting that along with the differential equations associated with unexcited and singly-excited norm kernels derived above, W 2 (ϕ) is shown to satisfy the first-order differential equation where the connected doubly-excited kernel of A 02 connects it to W 3 (ϕ).

n-tuple amplitude equation
As for single and double amplitudes, the derivation of the differential equation satisfied by the n-tuple amplitude invokes all those lower ranks. Reasoning by recurrence, one can prove where the connected n-tuply excited kernel of A 02 connects it to W n+1 (ϕ).

Initial conditions
The ODEs need to be complemented with initial conditions. Working up to the n-tuple level, one most naturally employs the initial conditions for k ≤ n and W k (0) = 0 for k > n, where the gaugeunrotated amplitudes U k are to be determined by solving standard BCC equations (Eq. 26) at the same n-tuple level.
Appendix D: Off-diagonal kernels

Results
The connected part h(ϕ) of the off-diagonal hamiltonian kernel and the norm kernel N (ϕ) are displayed as a Before discussing the results, it is worth noting that in an exact 14 limit of PBCC, the connected part of the off-diagonal hamiltonian or particle number kernel is independent of the gauge angle [15,16] while the norm kernel is nothing but N (ϕ) = e iAϕ . Contrarily, the nontrivial impact of the particle-number projection precisely relates to the gauge-angle dependency acquired by h(ϕ) and a(ϕ) as well as to the departure of N (ϕ) from a single IRREP of the U (1) group as a result of the underlying symmetry breaking.
The results obtained for at the well-celebrated PBCS ("MF") level provides a reference. The norm of h(ϕ) displays a typical bell-shaped curve that is maximum at ϕ = 0 15 . The norm kernel displays a similar qualitative behavior with N (0) = 1 as a result of Φ|Φ . The phase of both h(ϕ) and N (ϕ) are strictly linear in ϕ but while it is zero for h(ϕ) that is real, the slope is equal to A for N (ϕ).
Moving now to PBCS-CCSD, let us first focus on the case where T 1 (ϕ) and T 2 (ϕ), from which h(ϕ) is computed according to Eq. 19b, are obtained by solving the ODEs introduced in Sec. II E 2. The behavior of h(ϕ) is qualitatively very similar to what is observed at the PBCS level, i.e. when U = W (ϕ) = 0. While the phase remains strictly null, the norm h(ϕ) is more spread out with a maximum still located at ϕ = 0 16 . By going from PBCS to PBCS-CCSD, one starts observing the flattening of h(ϕ) and the converging towards the exact solution that must occur as one includes higher-order cluster operators 17 . As for h(ϕ), the norm of N (ϕ) resembles the one at the PBCS level but is flatter as a result of included many-body correlations 18 . The fact that N (0) = 1 remains true testifies of the intermediate normalization satisfied by the BCC wave function at any truncation order.
Let us now focus on the case where T 1 (ϕ) and T 2 (ϕ) are obtained by solving the gauge-dependent BCC amplitude equations introduced in Sec. II E 1. One observes that h(ϕ) is constant and real 19 , i.e. h(ϕ) = h(0) for all ϕ. This behavior is thus entirely different from what is obtained at the lowest order, i.e. PBCS level, which is both surprising and at odd with the results obtained from the ODEs. An analytical proof of this unexpected result is given in Sec. D 2 below. By virtue of Eq. C10, the norm kernel behaves trivially as N (ϕ) = e ia(0)ϕ = e iAϕ is this case. As mentioned above, the actual independence of h(ϕ) on the gauge angle renders the particle number projection trivial, i.e. inactive, even though the cluster operator U is truncated such that the symmetry is indeed broken at the BCC level. While requiring that T (ϕ) satisfies angle-dependent BCC equations is correct in the exact limit, it is not pertinent as soon as U is trun- cated. Indeed, this demand happens to break completely the link to the original definition of T (ϕ) in Eq. 18 in which its dependence on the gauge angle is driven by a gauge rotation. In conclusion, the determination of T (ϕ) via angle-dependent BCC amplitude equations is inappropriate and must be discarded in actual calculations where the symmetry is broken and needs to be effectively restored.

Constant connected hamiltonian kernel
Using the disentanglement properties of the Lie algebra [26], the symmetry rotation operator R(ϕ) can always be decomposed exactly as R(ϕ) = e X(ϕ) e Y (ϕ) e Z(ϕ) where In other words, Z(ϕ) destroys a pair of quasiparticles, X(ϕ) creates a pair of quasiparticles, and Y (ϕ) preserves the number of quasiparticles. The precise details of the matrix elements need not concern us here, but note that Z(ϕ) is nothing but the Thouless operator introduced in App. A 1.
Because the Hamiltonian is invariant under gauge rotation, so that To simplify the Hamiltonian kernel, we can definẽ BecauseT (ϕ) contains only quasiparticle creation operators, it commutes with X(ϕ), and it will prove convenient to also introducê T (ϕ) = X(ϕ) + e Y (ϕ) T (ϕ) e −Y (ϕ) .
Note thatT (ϕ) contains all the same excitation levels as does T (ϕ), and also contains single excitations even if T (ϕ) does not, simply due to the presence of X(ϕ). Now, with these definitions in hand the Hamiltonian kernel is H(ϕ) = Φ|H e X(ϕ) eT (ϕ) e Y (ϕ) |Φ e W0(ϕ) (D9a) where in the second line we have noted that e Y (ϕ) |Φ = |Φ e Y0(ϕ) , (D10) and the connected Hamiltonian kernel h(ϕ) defined in Eq. 19b is just This identity is true regardless of the truncation of the broken symmetry cluster operator U = T (ϕ = 0). Thus far all we have done is to note that, regardless of the truncation of the broken symmetry cluster operator, the reduced Hamiltonian kernel can be written in terms of the untransformed Hamiltonian and an alternative angle-dependent cluster operatorT (ϕ). Our next step is to show that if T (ϕ) satisfies the angledependent BCC equations, thenT (ϕ) satisfies the angleindependent BCC equations, and the reduced Hamiltonian kernel is therefore constant.
which is nothing but the BCC energy associated with T (0).
To restate our result, regardless of the truncation of U = T (0), so long as is contains single excitations (be-causeT (ϕ) always contains single excitations) and so long as all BCC equations for the relevant excitation levels are satisfied, one solution of the angle-dependent BCC equations is such that the reduced Hamiltonian kernel is constant and equal to its value at ϕ = 0, which is the unprojected BCC energy.