Löwdin's symmetry dilemma within Green functions theory for the one‐dimensional Hubbard model

The energy gap of correlated Hubbard clusters is well studied for one‐dimensional systems using analytical methods and density‐matrix‐renormalization‐group (DMRG) simulations. Beyond 1D, however, exact results are available only for small systems by quantum Monte Carlo. For this reason and, due to the problems of DMRG in simulating 2D and 3D systems, alternative methods such as Green functions combined with many‐body approximations (GFMBA), that do not have this restriction, are highly important. However, it has remained open whether the approximate character of GFMBA simulations prevents the computation of the Hubbard gap. Here we present new GFMBA results that demonstrate that GFMBA simulations are capable of producing reliable data for the gap which agrees well with the DMRG benchmarks in 1D. An interesting observation is that the accuracy of the gap can be significantly increased when the simulations give up certain symmetry restriction of the exact system, such as spin symmetry and spatial homogeneity. This is seen as manifestation and generalization of the “symmetry dilemma” introduced by Löwdin for Hartree–Fock wave function calculations.

flexibility of the HF wavefunction, giving an approximate energy value closer to the exact one, but typically does not preserve the system's good quantum numbers and symmetry properties. Well-known examples are HF simulations of the ground state of the uniform electron gas, e.g., Refs. [2][3][4]. From Löwdin's original insight, extensive research has spurred, to develop approaches where symmetries in a single-determinant wavefunction are deliberately broken, and subsequently reintroduced via symmetry projection operators, to attain a variationally improved multi-determinant state (see, e.g., Ref. [5] for a recent discussion).
Ambivalence in the use of symmetry is in fact of very general occurrence, and concerns both finite and infinite systems. An interesting example is provided by spontaneous symmetry breaking (SSB). [6] In a rigorous sense, SSB only takes place in the thermodynamic limit. However, exact numerical evidence from finite systems shows that the "disjointness" (typical of macroscopic sizes) of phases or different values/orientations of the order parameter is replaced by a crossover behaviour across finite barriers in Hilbert space, whose sharpness and strength increases on enlarging the system's size. An example is Wigner crystallization in finite electron systems, such as quantum dots, e.g., Refs. [7][8][9], for which also SSB in HF calculations was investigated; [10] other examples are found in ultracold bosons in traps, nuclear matter, and quantum chemistry; for an overview, see Ref. [11]. It can thus be methodologically expedient to artificially break the symmetry in a finite system, to gain insight about the system behaviour in the thermodynamic limit. An often used prescription is the addition of small external sources lowering the symmetry, [12,13] but under the stipulation that it is understood that true SSB occurs only asymptotically.
Another central element to consider in addressing the symmetry-related behaviour of a system is electronic or inter-particle correlations. These have deep influence in various situations, e.g. condensed-matter systems and materials, plasmas, nuclear matter, and cold atoms, to mention a few. [14][15][16][17] Clearly, the interplay of electronic interactions and symmetry constraints affects the system's properties in a way that is not accountable for within a free-particle or mean-field picture. It should be noted, though, that already within a wavefunction framework, some theories going beyond the mean-field picture can mimic the effect of strong electronic correlations with wavefunctions that do not respect the expected symmetry (see, e.g., Ref. [18]).
In this work, we take a different route from wavefunctions, and we study the effect of lifting/breaking symmetry in the presence of significant electronic correlations within many-body perturbative Green functions theory. As the system of choice with strong correlations, we consider the Hubbard model [14,19,20] which via a minimum-complexity Hamiltonian, describes the key trends in the behaviour of interacting electrons in the energy bands of a solid. For this reason, Hubbard-like models have been applied in many contexts and to a wide typology of systems, both in and out of equilibrium, see, e.g., Refs. [21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37].
Even though the Hubbard Hamiltonian is considerably simpler than that of a realistic material, exact solutions for the Hubbard model are only known in special cases: in one dimension, an exact analytical treatment is possible via Bethe-ansatz techniques, [38] and exact numerical solutions for finite samples can be obtained via the density-matrix-renormalization-group (DMRG) method [39][40][41] (including spin-charge separation effects [42][43][44] ). Using configuration interaction (CI) [45] or quantum Monte Carlo methods, [46][47][48][49] exact solutions can be obtained for any dimensionality, but only for small clusters. Finally, an exact description is possible in the limit of infinite dimensionality via dynamical-mean-field theory (DMFT). [50][51][52] In all other cases (notably, in two and three dimensions), some level of approximation must be introduced. Yet, high accuracy can be attained, for example, via the diagrammatic Monte Carlo technique, [53] or by extensions of DMFT via cluster [54] or diagrammatic approaches. [55,56] A premier method traditionally applied to the Hubbard model is the Green functions formalism combined with many-body perturbation theory (GFMBA). [57,58] The GFMBA method is a general framework that can be used in any dimension (i.e. also for 2D and 3D Hubbard models), scales not too unfavourably with system's size, can deal with both static and time-dependent regimes, [17,59] and is also practically viable for implementation for realistic systems. Furthermore, a very recent reformulation of the method in terms of coupled one-and two-particle propagators [60,61] has considerably increased its scope and range of applicability. In the case of the GFMBA method, correlation effects are included via selected classes (possibly infinite sums) of diagrams in the selfenergy, or via truncated iterative functional-derivative schemes. This leads to different perturbative treatments, [57,58] e.g. HF, second-order Born approximation (SOA), third-order approximation, [62,63] GW, and the particle-particle and particle-hole T-matrix approximation. [17,63] Comparisons against exact benchmarks for finite systems have shown that the GFMBA method works well for not too strong interactions. [62][63][64][65] The GFMBA approach has recently been under extensive scrutiny in relation to the existence of multiple solutions in the ground state and a potential convergence to nonphysical ones [66][67][68] (for an example from out-of-equilibrium systems, see e.g. Ref. [69]). However, in the discussion the multiplicity of GFMBA solutions, little or no attention has been given so far in the literature to the explicit role of symmetry, particularly in relation to the Hubbard model (for some discussion of broken symmetry in the case of an isolated Hubbard dimer, see Refs. [70,71]). For the exact solution of this model, important ground-state properties and symmetries are well known [72,73] and in common practice, these symmetries are granted by default also to the approximate solutions.
There is scarce knowledge on how these symmetries affect approximate treatments beyond mean-field theory. For example, for Hartree-Fock treatments of the Hubbard model, it is well known that a phase transition from a paramagnetic ground state to an antiferromagnetic one of unphysical nature occurs at a critical interaction, U c , where the specific value depends on the system size and geometry. [47,74] Yet, at the same time, the (unphysical) broken spin-symmetry solutions result in a ground-state energy closer to the exact one, as well as in the emergence of a band gap (the latter is absent within spin-symmetric/spin-restricted mean-field schemes). This raises the question whether selfenergy approximations beyond Hartree-Fock can be found that violate the symmetry properties of the exact solution as well, and yet provide "improved" values of relevant observables, such as the ground-state energy or the Hubbard gap. In other words, 1. Is there a Löwdin symmetry dilemma for the Hubbard model within many-body perturbation theory? 2. And in case, how is this related to solution multiplicity?
Our answer to the first question is in the affirmative: By considering the Hubbard model in the one-dimensional case, and comparing GFMBA, CI, and DMRG results, we find that, lifting the symmetry constraints artificially, simulates the effect of having more correlation effects in the system, and leads to a significant improvement of observables like the system's ground-state energy or the spectral functions, even at fairly large interactions. In particular, our discussion will be especially focussed on the charge energy gap Δ (also known as Mott-Hubbard gap). For a system containing N particles, it depends on the ground-state energies of the N, N + 1 and N − 1 systems, This is a central quantity of the Hubbard model, related to the metal-insulator transition at a characteristic interaction strength. 2 For the second question, we find that the occurrence of multiple (metastable) solutions is central to the connection between symmetry lifting and improved values of certain observables in the Hubbard model. Besides looking at what happens when lifting symmetry and why, we will also consider the possible implications in physical terms. However, due to the explorative nature of this initial study, we will not explore/discuss strategies to restore symmetry. This is left to future work.
The paper is organized as follows: in Section 2, we introduce the Hubbard Hamiltonian and the Green functions formalism for HF and SOA treatments within the three self-consistency protocols with different levels of symmetry constraints. Section 3 presents results from a mean-field treatment, specifically for the ground-state energy, the Hubbard gap, and the magnetic moment as function of the interaction strength U. In Section 4, we consider the SOA, discussing at the same time the ground-state density matrix and equilibrium spectral function for all the three self-consistency prescriptions. The multiple solutions in SOA are further analysed in Section 5 in terms of ground-state energy values and self-consistency convergence errors. Finally, Section 6 focuses on a characterization of the Hubbard gap, where HF and SOA results are compared to DMRG ones, and the dependence on system size, L, is taken into account. Our conclusions and a brief outlook are presented in Section 7.

GREEN FUNCTIONS THEORY
We consider the Fermi-Hubbard model which is described by the Hamiltonian 2 The characteristics of the MIT (and thus of the gap) depend on dimensionality (see, e.g., Refs. [14,75]). For D = 1, the exact solution shows that the insulating phase (and the charge gap) exist for any U > 0 (i.e. no MIT occurs). [38] For infinite dimensions, where the model is exactly solvable via DMFT, the MIT occurs at finite interaction values [52] (however, there is coexistence of metallic and insulating solutions in an interval U c 2 < U < U c 1 ). This picture qualitatively remains in three dimensions, with however quantitative modifications, due to antiferromagnetic fluctuations. [48] However, antiferromagnetic correlations appear to play an even more dramatic role in two dimensions, where it has been recently also proposed in Ref.
[76] that the system is gapped for any U > 0. For a recent review of the situation in 2D see e.g. Ref. [77].
where J is the hopping amplitude between adjacent lattice sites, and U is the on-site interaction strength. The operatorŝ c † i, andĉ j, create and annihilate an electron with spin projection at site i and j, respectively, and the density operator is given byn i, =ĉ † i,ĉ i, . The one-body nonequilibrium Green function of the system is defined by the canonical operators for complex times z on the Keldysh contour , [78,79] where   is the time-ordering operator on the contour, and the averaging is performed with the correlated unperturbed density operator of the system. In the rest of this work, we specialize to the equilibrium regime, i.e. when the system is not acted upon by external fields. In that case, the real-time components [17,59] of the Green function depend only on the difference of the two time arguments. Further, the retarded (R) component for the spin projection obeys the Dyson equation where all quantities are matrices in the orthonormal single-particle basis {|i⟩}. For spin-compensated situations, the noninteracting retarded Green function G R 0 ( ) is independent of the spin projection and given by depending only on the single-particle contribution to the Hamiltonian where the sum over ⟨i, j⟩ ensures only nearest neighbour hopping. For numerical reasons, a finite value of = 10 −2 is used throughout this work. If the exact selfenergy R ( ) of the system was known, Equation (4) would provide the exact single-particle Green function. However, in practice many-body approximations to the selfenergy have to be used. In this work, we employ the time-diagonal Hartree-Fock as well as the time-non-local second-order Born approximation. The retarded component of the HF selfenergy is defined as 3 where is the Dirac delta function for the relative time t in equilibrium and diag(⋅) represents a diagonal matrix with the given arguments as diagonal entries. Further, denotes the spin-projection opposite to , L the number of lattice sites, and The density matrix is given by the less component of the Green function, For the correlated Green function, the less (<) and greater (>) components can be determined by with the Fermi function f For the SOA selfenergy, the retarded component is given by with the Heaviside step function Θ(t) and the greater and less components of the selfenergy Here, • denotes the Hadamard product between matrices, and the G ≷ (t) are determined by the inverse Fourier transform, Finally, using the Fourier transform, , the selfenergies, Equations (7) and (11), can be included in Equation (4).
Since the selfenergy, in general, is a functional of the single-particle Green function, Equation (4) has to be solved iteratively for both spin components until a self-consistent solution is found. The choice of a suitable initial value is crucial and can affect the final result of the iteration. Here, if not mentioned otherwise, G R 0 ( ) is chosen as the starting point. In summary, for the self-consistent solution of Equation (4) the following scheme is iterated until convergence is achieved: 1. Diagonalize the single-particle Hamiltonianĥ 0 , cf. Equation (6), set G R 0 ( ) via Equation (5) and choose an initial value to start the iteration, e.g. (9) and (10) 3. Perform the inverse Fourier transform G , cf. Equation (13) 4. Calculate ≷ (t) and R (t) using Equations (7), (11) and (12) 5. Perform the Fourier 6. Solve the Dyson equation for G R ( ), Equation (4), using the new R ( ) 7. If G R ( ) is not yet converged start again at (1) To improve the convergence of the above scheme, the input Green function at iteration k, G R k,in (the spin index is neglected here), is determined by mixing the solutions of the two previous iterations where a mixing parameter of = 0.05 is used. The error at iteration k is given by where D( ) is the density of states (DOS) of the system, During the calculations, it is verified that the sum rules for the spectral function hold, i.e. that the DOS remains positive and integrates to the number of lattice sites L within numerical accuracy. To ensure the convergence of the iteration scheme an error threshold of thr = 10 −12 is used in this work. In addition to spectral properties, the single-particle Green function gives access to the total energy of the system. [17,57] which combines the kinetic part and the Galitskii-Migdal interaction energy where, in the GFMBA simulations of the present paper, the chemical potential is set to = 0. In this paper, we consider three distinct cases for which different symmetry restrictions are imposed on the Green function during the solution of the Dyson equation: 1. "uniform" (uni): the system is required to be translationally invariant. In this case, the iteration scheme is solved in momentum space where the Green function and selfenergy are diagonal, i.e. G R ij, ( ) → G R p, ( ). 2. "restricted spin" (rs): the system is required to be spin-symmetric. In this case, the iteration scheme is solved for one spin projection only since the Green function and selfenergy are spin-independent, i.e. G R ↑ ( ) = G R ↓ ( ). 3. "unrestricted spin" (us): no restrictions regarding both, the translation and spin symmetry are imposed. In this case, an antiferromagnetic state is chosen to start the iteration. 4

SPIN SYMMETRY IN THE MEAN-FIELD APPROXIMATION
Since the exact ground state of the half-filled 1D Hubbard Hamiltonian is known to be spin-symmetric (i.e. paramagnetic) for systems with an even number of particles, [72,73] a logical prescription is to introduce spin symmetry also for approximate solutions like HF and SOA. However, it is well known that, beyond a critical interaction strength U c , unrestricted-spin HF (usHF) spontaneously breaks spin symmetry resulting in an antiferromagnetic ground state. In the following, we quantify the influence of this artificial phase transition on important ground-state properties by comparing the performance of restricted-spin (rs) and unrestricted-spin (us) HF for finite one-dimensional Hubbard chains with open (hard-wall) boundary conditions. In Figure 1a the ground-state energy of three Hubbard clusters containing L = 2, 4, 6 sites is plotted vs. the interaction strength U for rsHF, usHF, and the result obtained by exact diagonalization of the Hamiltonian. The qualitative observations are similar for all three systems. In the limit of vanishing on-site interaction all three methods agree perfectly and show a linear increase of the ground-state energy with U. For interactions beyond U ≳ 1J, the exact energy is reduced due to increasing correlations giving rise to mounting differences compared to the two HF solutions. In the case of rsHF, which by design fulfils the exact spin symmetry of the system, the linear increase of the ground-state energy is present for all values of U resulting in a strong deviation from the exact result, for U ≳ 1J. Most notably, since correlations are not included in rsHF, no Mott regime is observed in the presence of the on-site interaction. By contrast, removing the requirement of spin-symmetry (usHF) results in a lower ground-state energy for interactions beyond U c ≈ 2 J which approaches the exact value for U → ∞. Additionally, the usHF density of states, shown in Figure 1c for a Hubbard chain of 10 sites, is indicative that a correlation gap in the spectrum (corresponding to a Mott transition) emerges for a critical interaction U c . However, since the usHF selfenergy accounts only for mean-field effects, the improved results for the Hubbard gap and the ground-state energy cannot be attributed to the effects of correlations. Instead, they are connected to the emergence of the antiferromagnetic state. This becomes apparent when looking at the local magnetic moment on the outermost site of finite Hubbard chains depicted in Figure 1b. The local magnetic moment is defined as with n i = n i,↑ + n i,↓ and the local double occupancy d i which contains a mean-field and a correlation part, 4 When choosing the spin-symmetric G R 0 ( ) as the initial value, the iteration will not break spin symmetry. In the exact case, where the spin densities are homogeneous, an increasing magnetic moment at high interaction strengths is caused by an increase of electronic correlations leading to a negative d i corr and, thus, a decrease in the double occupancy. However, for usHF, where d i corr ≡ 0, the inhomogeneous spin-density distribution of the antiferromagnetic spin state mimicks the effect of additional correlations. 5 To summarize, removing the requirement of a homogeneous spin-symmetric ground state, as observed in the exact solution, allows HF simulations to achieve results for ground-state energies closer to the exact ones and gives rise to a Hubbard gap of reasonable magnitude.

SYMMETRIES IN SECOND-BORN APPROXIMATION
We next analyse the effect of imposing symmetry restrictions when including selfenergies with correlation effects (beyond HF). The first correction beyond HF that takes into account interparticle scattering is the second-order Born approximation (SOA). We thus consider, within SOA, the behaviour of finite Hubbard chains with periodic boundary conditions for which the exact ground state is known to be both, spin symmetric and invariant under space translations. [72,73,80] The effect of relaxing the aforementioned symmetry constraints is quantified by comparing results from the three SOA approaches I, II, and III, (i.e. for the uniform, restricted-spin, and unrestricted-spin treatments) introduced in Section 2; we will refer to these three treatments as uniSOA, rsSOA, and usSOA, respectively. For all three cases, we compare in Figure 2 the spin-up density matrix (a-d), the DOS (e-g) and the ground-state energy (h) to the exact results for an 8-site Hubbard chain with periodic boundary conditions at U = 4 J. In the uniform case, the density matrix by design exhibits perfect 5 Note that the critical interaction U c for which the symmetry-broken state emerges decreases with increasing length of the chain. An in-depth analysis on the exact value of U c goes beyond the scope of this work. Relaxing the requirement of translational symmetry (rsSOA) leads to unphysical inhomogeneities in the odd minor diagonals of the spin-up density matrix which is in contrast to the exact solution. 6 However, the ground-state energy is improved, to −3.84 J, closer to the exact result. Additionally, in the DOS, a correlation-induced gap emerges at the Fermi energy in conjunction with an, in general, better qualitative agreement with the exact spectrum. Still, the rsSOA gap of 0.77 J is less than half the size of the exact result of 2.01 J.

F I G U R E 2
As a next step, we no longer enforce spin symmetry (usSOA) which results in an antiferromagnetic ground state indicated by the spin-density wave on the diagonal and the inhomogeneities on the even minor diagonals of the spin-up density matrix, shown in Figure 2c. This Néel state has an energy of −4.08 J which is in much better agreement with the exact value. In this case, the DOS nicely reproduces the position of the main peaks and the Hubbard gap of 1.77 J is much closer to the exact value. 7 Similar to the findings for the HF selfenergy presented in Section 3, removing the requirement of translation and spin symmetry for the SOA leads to a significant improvement for the ground-state energy and the DOS. This way the exact Mott gap can be remarkably well reproduced, even for relatively large interactions such as U = 4 J. 6 Note that, while the translational symmetry is broken, the spin symmetry is still fulfilled. 7 The remaining differences to the exact DOS, namely, the missing high-energy satellites and the degenerate peak at ∼±1.8 J, can be attributed to the shortcoming of the SOA at the large interaction U = 4 J.

MULTIPLE SOLUTIONS OF THE DYSON EQUATION IN SECOND-BORN APPROXIMATION
The possibility of multiple solutions is a well-known feature of self-consistent treatments of the Dyson equation. [67,68,81] Usually, a self-consistency requirement is employed in approximate treatments via e.g. perturbation theory. The existence of a self-consistent solution to the Dyson equation was shown to depend on the choice of the selfenergy approximation, in particular its positive semi-definiteness (PSD). The SOA selfenergy used in this work is known to be conserving and PSD, [82] and is thus suitable to study the emergence of multiple solutions. With a specific choice of the translational invariance, and without addressing the role of broken symmetry, this was done before in a specific case. [83] Here, the effect of artificially breaking specific symmetries is brought to the foreground, and explicitly connected to the insurgence of multiple approximate self-consistent solutions. To this end, we solve Equation (4) with no symmetries enforced (usSOA), for a half-filled Hubbard chain of length L = 8 with periodic boundary conditions and U = 4 J, i.e., for the same system as in Section 4. The initial state of the iteration scheme is chosen to be the homogeneous, spin-restricted HF ground state where the density is modified by a small random perturbation of the order 10 −5 . 8 In Figure 3, the DOS (a,b), ground-state energy (c), and iteration error (d) are shown during the iterative procedure. The calculation starts from the slightly disturbed rsHF ground state with energy −1.52 J (not shown in Figure 3c). Within the first 170 iterations, the system converges into the homogeneous state (cf. the DOS in Figure 2e), as the relative iteration error drops to 10 −3 . However, at around 300 iterations the error increases and the system transitions into the inhomogeneous but spin-symmetric state (cf. the DOS in Figure 2f). This state appears to be even more stable, with the relative iteration error temporarily dropping to 10 −5 . However, after 3000 iterations, a final transition sets in, and the system arrives in the spin-asymmetric state (cf. the DOS in Figure 2g). The system remained in this state for the remainder of the iteration process, and the relative error of the calculation eventually reached the order of machine precision. 9 While, during the iteration, the DOS evolves through the three states shown in Figure 2, the total energy passes through the values of the respective states that are depicted in Figure 2h. This is shown in detail, as a function of the iteration number, in Figure 3c, cf. the coloured lines. This shows that the three states, corresponding to the different symmetry restrictions discussed in Section 4, can be reached by a single iterative solution of the Dyson equation when no symmetries are enforced. During the iteration, in the vicinity of each of the three states, the iteration error drops significantly (reaching a high degree of self-consistency). This suggests that all of them are solutions of the same Dyson equation, with only the unrestricted-spin state being absolutely numerically stable. A well-conditioned iteration scheme will, ultimately, reach this minimum-energy state. Of course, the sequence of states reached on the way to the minimum depends on the choice of the initial state of the iteration and on details of the numerical procedure.
The present example is a direct illustration of Löwdin's symmetry dilemma discussed above and shows that a successive reduction of the symmetry of explored states may allow one to improve certain target properties of a system, such as the ground-state energy, also in a many-body calculation with the SOA selfenergy. The behaviour just discussed here is expected to be a specific facet of a general scenario underlying the search for multiple solutions of the Dyson equation that are a consequence of the nonlinear dependence of the collision integrals on the Green functions which is a general property of selfenergies beyond Hartree-Fock. These observations should give useful hints how to improve iterative solutions of the Dyson equation or similar nonlinear equations.

BENCHMARKING AGAINST DMRG
The existence of a Mott gap for large on-site interactions is one of the most important features of the Hubbard model. After analysing the general properties of GFMBA simulations with HF and SOA selfenergies for the three methods (I)-(III), we now focus on their performance regarding the Hubbard gap, Equation (1), in particular. Since we are considering finite systems we are interested in the correlation part of the band gap which we define as 8 Recall that, without this small initial inhomogeneity, no broken spin symmetry would occur, even in spin-resolved calculations. 9 It should be noted that there is a little dip/kink in the relative error at around iteration 3500. This could hint to a possible forth viable state that we did not reach in our calculation. where Δ rsHF is the band gap obtained from a rsHF calculation that contains only the finite-size contribution. In the thermodynamic limit Δ rsHF vanishes, and Δ corr = Δ. In Figure 4, we compare the correlation gap of finite Hubbard chains of varying length at U = 4 J to the (exact) result obtained by DMRG (we employed the size-increasing scheme as in Refs. [84,85]). Additionally, we extrapolate the data to L → ∞ where the DMRG result agrees with the Bethe-ansatz solution. [86] In the case of the restricted-spin HF and the homogeneous SOA (uniSOA) state the Hubbard gap vanishes, cf. Figure 2d for the latter.

F I G U R E 3
In contrast, starting with open boundary conditions, rsSOA shows a finite gap and correctly predicts its qualitative dependence on the length of the system. However, since SOA captures only part of the correlation effects, the correct band gap is underestimated by ∼0.7 J for all system sizes. As discussed for the local magnetic moment in Section 3, the antiferromagnetic ground state of the unrestricted-spin methods can compensate shortcomings of the selfenergy approximations in treating correlations. In the case of usHF, this results in the opening of a correlation gap on the mean-field level. However, for the large interaction of U = 4 J, the size of the gap is severely overestimated, especially, in the thermodynamic case, where usHF predicts a size of 3.073 J, as opposed to the exact value of 1.287 J. In contrast to the exact case, the correlated gap is monotonically increasing with the system length for the unrestricted-spin methods. Nevertheless, usSOA shows the best agreement with the exact gap out of all selfenergies and symmetry restrictions considered here. Especially for large system sizes, including the case of the infinite Hubbard chain, the deviation does not exceed ∼0.25 J.
The results for periodic boundary conditions (shown as dashed lines in Figure 4) differ only slightly from the above observations. While, for small systems (L < 40), there is a noticeable difference between both cases, the size of the correlated gap converges to the same value for larger systems. The speed of this convergence is considerably faster for the unrestricted-spin methods which possess an antiferromagnetic ground state.
Independent of the type of boundary conditions, giving up on the symmetries of the exact solution, the description of the Hubbard gap in finite systems, especially for the selfenergy in SOA, improves dramatically, even for large interaction strengths, again confirming Löwdin's symmetry dilemma. with crosses) and periodic (dashed lines with circles) boundary conditions. Density-matrix-renormalization-group results are compared to different restricted-spin (rs) and unrestricted-spin (us) selfenergy approximations. The extrapolated results for the limit of the infinite Hubbard chain are shown as dotted lines on the right To further assess the quality of not only the correlation gap but the total DOS, in Figure 5 we compare the results of the rsSOA and usSOA methods to the spectral function obtained by Nocera et al. [87] using the time-dependent DMRG (tDMRG) method for an open Hubbard chain of 40 sites for U = 4 J. The finite real-time window of the tDMRG approach results in an artificial broadening of the peaks in the spectral function. For a better comparison, we emulate this effect by performing on our much sharper DOS a convolution with a Gaussian that corresponds to a time-window size of 15 J −1 . The general trends observed in Figure 2f,g for a short chain of eight sites are confirmed also for this larger system of length L = 40. The main characteristics of the tDMRG spectrum are the major peaks at ±1 J and ± 2 J, the satellites at ±4.5 J and a slowly descending slope between those features. According to the observations of Figure 4, rsSOA severely underestimates the size of the Hubbard gap since the spectrum in general is shifted toward the Fermi energy resulting in a poor overall agreement with tDMRG. On the contrary, usSOA reproduces the position of the major peaks and the satellites remarkably well, whereas only the slope between them is not described correctly. The spectral weight that is missing in the slope is, instead, transferred to the oversized major peaks. Nevertheless, the agreement of usSOA with the tDMRG results is striking, considering the high on-site interaction of U = 4 J.

CONCLUSION
Currently, there is great interest in the theoretical condensed-matter community in devising approaches for strongly correlated systems, i.e. for those systems where a description based on an independent-particle picture is qualitatively inadequate.
In this work, we have looked at a certain aspect of the problem, namely, the interplay of symmetry constraints and electronic correlations. Specifically, using finite Hubbard chains at half filling as a case study for strongly correlated systems, we investigated the role of symmetry requirements in many-body perturbation Green function theory (GFMBA), where approximations of increasing complexity can be systematically devised for the many-body selfenergy.  [87]. The rsSOA and usSOA spectra are broadened using a Gaussian corresponding to a real-time window of size 15 J −1 that matches the width in Ref. [87] In the literature, GFMBA is often seen as an ill-suited conceptual paradigm to describe strong correlations, because of, e.g., possible multiple self-consistent approximate solutions, or uncertainty about the convergence radius of the perturbation expansion. While not calling into question this point of view, in this study we have used GFMBA to address the interplay of symmetry and correlation effects. To our knowledge this is a point that, irrespective of the methodology used, has received little systematic attention so far. For our GFMBA description of Hubbard systems, we used the second-Born approximation which accounts for electronic correlations at lowest perturbative order in a "skeleton-diagram" sense. However, for the sake of comparison we also presented results from an Hartree-Fock treatment.
Already at the HF level, our comparison of spin-symmetry-restricted and unrestricted self-consistent solutions indicates that one is facing a so-called Löwdin symmetry dilemma for the one-dimensional Hubbard model: namely, the violation of spin-symmetry can lead to a spectral function and ground-state energy that are closer to the exact ones than those obtained when spin symmetry is imposed.
Our results suggest that this behaviour is robust against the inclusion of electronic correlations within GFMBA. For our self-consistent SOA treatment of finite periodic Hubbard chains at half-filling, in addition to spin symmetry we also considered translation symmetry. While the exact solution corresponds to densities which are both spin-projection symmetric and translationally invariant, the violation of these properties in SOA leads to results for the ground-state energy and the local spectral function which are remarkably close to the exact ones. This is particularly so when using the unrestricted-spin symmetry SOA: in this case, the exact value of the Hubbard energy gap is surprisingly well produced, within ∼15%, also for strong on-site interaction U = 4 J. Finally, an interesting overall trait of our results is that the symmetry-restricted states are in fact solutions (albeit metastable) to the unrestricted Dyson equation: starting from a specifically crafted initial state, the self-consistency iteration dynamics passes through the symmetry-constrained states before reaching the unrestricted ground state.
Thus, altogether our work illustrates a direct connection between symmetry constraints and solution multiplicity in GFMBA, adding to the already available body of knowledge on the behaviour of multiple solutions for the Dyson equation. As possible future directions, an obvious point to address is if more complex selfenergy approximations, such as the T matrix and GW approximations, would give improved results in symmetry-lifted treatments as well. Other straightforward extensions would be the exploration of the case of higher dimensions (where, however, an increased complexity is expected due to a richer structure of the phase diagram) and electron occupancies different from half filling.
More in general (and very much in the spirit of Löwdin's original lines of thought), it could be of some interest to see if it is possible (and what happens) when concretely profiting from the symmetry dilemma within GFMBA, by restoring symmetry at the end via projection techniques. This procedure, in different variants and extensions, has been already used in the literature for wavefunctions, typically starting from HF symmetry-unrestricted states, e.g., Ref. [11]. Pursuing the same strategy within GFMBA would allow one to explore the feasibility of convenient ways to access some important physical quantities in strongly correlated systems.

ACKNOWLEDGMENTS
Open Access funding enabled and organized by ProjektDEAL.

DATA AVAILABILITY STATEMENT
Research data are not shared