An extended Agassi model: algebraic structure, phase diagram, and large size limit

Symmetry plays a major role in physics to handle problems that because of its complexity or size cannot be treated by brute force diagonalizations or with standard many-body approximations. On the other hand, algebraic boson or fermion models with few-degrees of freedom are easily solved and have provided for years with benchmarks to test many-body approximations. These models are characterized by Lie algebras [1] and for particular situations (called dynamical symmetries) can be solved analytically. In addition, in more general situations they can be solved numerically even for very large system sizes. In different areas of physics these families of models can be found. For instance, in quantum optics popular models for radiation-matter interaction are the Jaynes–Cummings [2] or the Dicke [3] models. In nuclear physics, the interacting boson model (IBM) [4] is based on the su(6) and has been used extensively in the last years for the study of medium and heavy nuclei. Also in nuclear physics the Elliot SU(3) model [5] was proposed to study rotational structures in nuclei. Other simple two-level models as the Lipkin–Meshkov–Glick (LMG) model [6] and the two-level pairing model [7] were originally proposed in the context of nuclear physics but then their simplicity and success allow them to be used in different physics branches: condense matter, molecular physics or cold atom physics.

The Agassi model was proposed in the late 1960s by Agassi as a model that could mimic the much more complex pairing-plus-quadrupole model of nuclear physics [8]. To this end, he tailored a model through the interplay of a two-level pairing plus a monopole term. In this sense, the Agassi model is an extension of the LMG model which, in addition, contains a pairing interaction and, therefore, is a much more rich model. The model was proposed originally to test the goodness of two popular many-body approximations, namely, the Bardeen–Cooper–Schrieffer (BCS) and the Hartree–Fock–Bogoliubov (HFB) approximations. The Agassi model allows its exact diagonalization, that can be carried out easily owing the algebraic structure of the model [8]. These exact results can be confronted to those obtained with approximation techniques to validate under which conditions these approximations are sufficiently precise. In spite of its unquestionable attractive, the Agassi model has been seldom used for years regardless of its great flexibility and its simplicity to be solved for large systems. There have been few works published in which the random phase approximation, HFB [8–10] or perturbation theory [11] were applied to the Agassi model. Moreover, modern many-body approximations that have been used intensively in nuclear physics did not take advantage of the model to check their applicability, and to establish their limitations and accuracy. Only very recently, an application to explore the coupled cluster theory in relation to symmetry breaking and restoration [12] has been published. This work opens the possibility of future applications to quantum chemistry and nuclear physics.

In addition, the Agassi model can play a relevant role in the field of quantum phase transitions (QPTs). The study of QPTs and critical points in algebraic models have received much attention in the last two decades. Many publications on the subject have been done (see, e.g. [13] and [14]). In particular, QPTs have been studied within the IBM [15, 16], the LMG [17], the vibron model [18], and the Jaynes–Cummings [2] and Dicke [3] models. Therefore, it seems quite natural to use the Agassi model, which encompasses the LMG and the two-level paring models, with an underlying so(5) algebra, to study the QPTs that arise in it. In [10] the Agassi model phase diagram was established. This diagram presents, apart from the symmetric phase, two broken symmetry phases: one connected to parity breaking related to the monopole interaction and another, called superconducting phase, linked to the pairing interaction. Recently, an extended Agassi model, in which a more general monopole interaction is included, was proposed in [19]. The corresponding richer phase diagram was reported and it includes several QPTs of different character. The model was studied using the mean-field HFB theory and the results checked against exact diagonalizations.

In this paper, the extended Agassi model is revised and its boson image using the Schwinger representation is worked out. Then, the phase diagram of the model is discussed and a technique that allows to obtain the mean field results as the large size limit of exact finite number of particle calculations is introduced. The paper is organized as follows: in section 2, the algebraic structure of the proposed extended Agassi model is discussed. Also in section 2, the Schwinger representation of the model is worked out. In section 3, a mean-field approach based on a condensate of bosons plus the HFB approach is applied to obtain the corresponding mean-field energy surfaces and to analyze their stability. In this section, the phase diagram structure and the nature of the QPTs are presented. In section 4, a 1/j analysis of several quantities is proposed to obtain their thermodynamic limit. Finally, section 5 is devoted to the summary and conclusions of this work.

The Agassi model is a two-level system, each of the levels has Ω (even) degeneracy. The labeling for the single particle states is the following: σ = 1 for the upper level and σ = −1 for the lower one, plus an additional quantum number m = ±1, ±2, ..., ±Ω/2, which labels the states within a given level. Consequently, one can use a spin image taking j as an integer number and Ω = 2j. Moreover, σ is used as the level parity, σ = +1 and σ = −1 imply positive and negative parity, respectively. The extended Agassi model Hamiltonian is

$$\begin{eqnarray}\begin{array}{rcl}H & = & \varepsilon {J}^{0}-g\displaystyle \sum _{\sigma ,{\sigma }^{{\prime} }=-1,1}{A}_{\sigma }^{\dagger }{A}_{{\sigma }^{{\prime} }}\\ & & -\displaystyle \frac{V}{2}\left[{\left({J}^{+}\right)}^{2}+{\left({J}^{-}\right)}^{2}\right]-2{{hA}}_{0}^{\dagger }{A}_{0}.\end{array}\end{eqnarray} \tag{ 1 }$$

All operators in the Hamiltonian (1) are defined in terms of the fermion creation and annihilation operator as

$$\begin{eqnarray}&&{J}^{+}=\displaystyle \sum _{m=-j}^{j}{c}_{1,m}^{\dagger }{c}_{-1,m}={\left({J}^{-}\right)}^{\dagger },\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{J}^{0}=\displaystyle \frac{1}{2}\displaystyle \sum _{m=-j}^{j}\left({c}_{1,m}^{\dagger }{c}_{1,m}-{c}_{-1,m}^{\dagger }{c}_{-1,m}\right),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{A}_{1}^{\dagger }=\displaystyle \sum _{m=1}^{j}{c}_{1,m}^{\dagger }{c}_{1,-m}^{\dagger },\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{A}_{-1}^{\dagger }=\displaystyle \sum _{m=1}^{j}{c}_{-1,m}^{\dagger }{c}_{-1,-m}^{\dagger },\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{A}_{0}^{\dagger }=\displaystyle \sum _{m=1}^{j}\left({c}_{-1,m}^{\dagger }{c}_{1,-m}^{\dagger }-{c}_{-1,-m}^{\dagger }{c}_{1,m}^{\dagger }\right),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{A}_{1}=\displaystyle \sum _{m=1}^{j}{c}_{1,-m}{c}_{1,m},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{A}_{-1}=\displaystyle \sum _{m=1}^{j}{c}_{-1,-m}{c}_{-1,m},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{A}_{0}=\displaystyle \sum _{m=1}^{j}\left({c}_{1,-m}{c}_{-1,m}-{c}_{1,m}{c}_{-1,-m}\right),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{N}_{\sigma }=\displaystyle \sum _{m=-j}^{j}{c}_{\sigma ,m}^{\dagger }{c}_{\sigma ,m},\qquad N={N}_{1}+{N}_{-1},\end{eqnarray} \tag{ 10 }$$

where ${c}_{\sigma ,m}^{\dagger },{c}_{\sigma ,m}$ are single fermion creation and annihilation operators in the state $| \sigma ,m\rangle$, respectively. The last term in (1), $-2{{hA}}_{0}^{\dagger }{A}_{0}$, has been added to the original Agassi model. This term does not correspond to the standard pairing since it involves correlated pairs with opposite angular momentum projection, but sitting in different levels, while the standard one treats with correlated pairs sitting in the same level, either the lower or the upper one. The addition of this extra term gives rise to new relevant effects not present in the original model formulation. It is convenient to redefine the Hamiltonian parameters to have them rescaled with the shell size. Thus, new parameters χ, Σ, and Λ (see [10]) are introduced as

$$\begin{eqnarray}&&V=\displaystyle \frac{\varepsilon \chi }{2j-1},\quad g=\displaystyle \frac{\varepsilon {\rm{\Sigma }}}{2j-1},\quad h=\displaystyle \frac{\varepsilon {\rm{\Lambda }}}{2j-1}.\end{eqnarray} \tag{ 11 }$$

These three model parameters have to be positive in order to do not produce unphysical situations. With this parametrization, the extended Agassi Hamiltonian is

$$\begin{eqnarray}\begin{array}{rcl}H & = & \varepsilon \left[{J}^{0}-\displaystyle \frac{{\rm{\Sigma }}}{2j-1}\,\displaystyle \sum _{\sigma {\sigma }^{{\prime} }}{A}_{\sigma }^{\dagger }{A}_{{\sigma }^{{\prime} }}\right.\\ & & -\displaystyle \frac{\chi }{2(2j-1)}\,\left[{\left({J}^{+}\right)}^{2}+{\left({J}^{-}\right)}^{2}\right]\\ & & -\left.2\displaystyle \frac{{\rm{\Lambda }}}{2j-1}\,{A}_{0}^{\dagger }{A}_{0}\right].\end{array}\end{eqnarray} \tag{ 12 }$$

The model algebraic structure is o(5), with 10 independent generators: six A's, plus three J's, plus the particle number, N. It is worth noting that N₁ and N₋₁ are linear combinations of J⁰ and N. In a pictorial way the operators of the model are represented in figure 1. It is of interest to write down the commutation relations of the operators (2)–(10) that underlie an o(5) algebra

$$\begin{eqnarray}&&\left[{J}^{+},{J}^{-}\right]=2{J}^{0},\quad \left[{J}^{0},{J}^{\pm }\right]=\pm {J}^{\pm },\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\left[{A}_{1},{A}_{1}^{+}\right]=j-{N}_{1},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\left[{A}_{-1},{A}_{-1}^{+}\right]=j-{N}_{-1},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\left[{A}_{0},{A}_{0}^{+}\right]=2j-N,\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\left[{A}_{0},{A}_{1}^{+}\right]=-{J}^{+},\qquad \left[{A}_{0},{A}_{-1}^{+}\right]=-{J}^{-},\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\left[{A}_{1},{A}_{0}^{+}\right]=-{J}^{-},\qquad \left[{A}_{-1},{A}_{0}^{+}\right]=-{J}^{+},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\left[{J}^{+},{A}_{1}^{+}\right]=0,\,\left[{J}^{+},{A}_{-1}^{+}\right]={A}_{0}^{+},\,\left[{J}^{+},{A}_{0}^{+}\right]=2{A}_{1}^{+},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\left[{J}^{-},{A}_{1}\right]=0,\,\left[{J}^{-},{A}_{-1}\right]=-{A}_{0},\,\left[{J}^{-},{A}_{0}\right]=-2{A}_{1},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&\left[{J}^{+},{A}_{1}\right]=-{A}_{0},\,\left[{J}^{+},{A}_{0}\right]=-2{A}_{-1},\,\left[{J}^{+},{A}_{-1}\right]=0,\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&\left[{J}^{-},{A}_{1}^{+}\right]={A}_{0}^{+},\,\left[{J}^{-},{A}_{0}^{+}\right]=2{A}_{-1}^{+},\,\left[{J}^{-},{A}_{-1}^{+}\right]=0,\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&\left[{J}^{0},{A}_{1}^{+}\right]={A}_{1}^{+},\,\left[{J}^{0},{A}_{-1}^{+}\right]=-{A}_{-1}^{+},\,\left[{J}^{0},{A}_{0}^{+}\right]=0,\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&\left[{J}^{0},{A}_{1}\right]=-{A}_{1},\,\left[{J}^{0},{A}_{-1}\right]={A}_{-1},\,\left[{J}^{0},{A}_{0}\right]=0.\end{eqnarray} \tag{ 24 }$$

Therefore, the Hamiltonian (12) can be diagonalized within a o(5) basis [8, 20]. On the other hand, the system eigenstates have a well defined parity since the parity operator ${{\rm{e}}}^{-{\imath }\pi {J}^{0}}$ commutes with the Hamiltonian (12). The complete, but non-orthogonal, basis of the model is defined in terms of three quantum numbers: ${n}_{+}$ for the number of pairs in the upper orbit (σ = 1), ${n}_{-}$ for the number of pairs in the lower orbit (σ = −1), and n₀ for the number of A₀ pairs

$$\begin{eqnarray}&&\left|{n}_{+},{n}_{-},{n}_{0}\right\rangle ={\left({A}_{1}^{+}\right)}^{{n}_{+}}{\left({A}_{-1}^{+}\right)}^{{n}_{-}}{\left({A}_{0}^{+}\right)}^{{n}_{0}}\left|0\right\rangle .\end{eqnarray} \tag{ 25 }$$

In order to perform calculations within the Agassi model, one has to select the j quantum number, that fixes the system size, and the number of fermions, N in the space. Here, the ratio $N/4j$ is fixed to the value 1/2 ($N=2j$) which means half-filling or, equivalently, j fermion pairs are considered. In this half-filling situation, the quantum numbers fulfill ${n}_{-}+{n}_{+}+{n}_{0}=j$.

Zoom In Zoom Out Reset image size

2.1. Schwinger representation

It is of interest to build the boson image of the extended Agassi model through the Schwinger representation. To this end one introduces five bosons: ${p}_{-1}^{\dagger }$ which creates a correlated fermion pair in the lower level (σ = −1), ${p}_{1}^{\dagger }$ which creates a correlated fermion pair in the upper level (σ = 1), ${p}_{0}^{\dagger }$ which creates a correlated fermion pair with a given m (one fermion in the upper level and the other in the lower one), ${v}^{+}$ which creates a quartet of fermions (one pair $(m,-m)$ above and other pair $(m,-m)$ below), and ${u}^{+}$ that corresponds to the identity (vacuum). These five operators (${b}_{i}^{\dagger }\equiv {p}_{-1}^{\dagger },{p}_{0}^{\dagger },{p}_{1}^{\dagger },{u}^{\dagger },{v}^{\dagger }$) fulfill the usual boson commutation relations

$$\begin{eqnarray}&&[{b}_{i},{b}_{j}^{\dagger }]={\delta }_{{ij}},\,[{b}_{i}^{\dagger },{b}_{j}^{\dagger }]=0,\,[{b}_{i},{b}_{j}]=0.\end{eqnarray} \tag{ 26 }$$

The O(5) generators in terms of Schwinger representation are written as

$$\begin{eqnarray}&&{J}^{0}={p}_{1}^{\dagger }{p}_{1}-{p}_{-1}^{\dagger }{p}_{-1},\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{J}^{+}=\sqrt{2}\left({p}_{1}^{\dagger }{p}_{0}+{p}_{0}^{\dagger }{p}_{-1}\right),\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{J}^{-}=\sqrt{2}\left({p}_{-1}^{\dagger }{p}_{0}+{p}_{0}^{\dagger }{p}_{1}\right),\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{A}_{1}^{\dagger }={p}_{1}^{\dagger }u+{v}^{\dagger }{p}_{-1},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{A}_{0}^{\dagger }=\sqrt{2}\left({p}_{0}^{\dagger }u-{v}^{\dagger }{p}_{0}\right),\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{A}_{-1}^{\dagger }={p}_{-1}^{\dagger }u+{v}^{\dagger }{p}_{1}\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{A}_{1}={p}_{-1}^{\dagger }v+{u}^{\dagger }{p}_{1},\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{A}_{0}=\sqrt{2}\left({u}^{\dagger }{p}_{0}-{p}_{0}^{\dagger }v\right),\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{A}_{-1}={p}_{1}^{\dagger }v+{u}^{\dagger }{p}_{-1},\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{N}_{1}={p}_{0}^{\dagger }{p}_{0}+2{p}_{1}^{\dagger }{p}_{1}+2{v}^{\dagger }v,\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&{N}_{-1}={p}_{0}^{\dagger }{p}_{0}+2{p}_{-1}^{\dagger }{p}_{-1}+2{v}^{\dagger }v,\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&N={N}_{-1}+{N}_{1}.\end{eqnarray} \tag{ 38 }$$

The Hamiltonian (1) in terms of Schwinger bosons is written as

$$\begin{eqnarray}&&\begin{array}{l}{H}_{B}=\left(\varepsilon -g\right){p}_{1}^{\dagger }{p}_{1}-\left(\varepsilon +g\right){p}_{-1}^{\dagger }{p}_{-1}\\ \ -\,\left(g+V\right)\left({p}_{1}^{\dagger }{p}_{-1}+{p}_{-1}^{\dagger }{p}_{1}\right)-2{{gv}}^{\dagger }v-2h\left({p}_{0}^{\dagger }{p}_{0}+{v}^{\dagger }v\right)\\ \ -\,V\left({p}_{1}^{\dagger }{p}_{1}^{\dagger }{p}_{0}{p}_{0}+{p}_{-1}^{\dagger }{p}_{-1}^{\dagger }{p}_{0}{p}_{0}\right.\\ \ +\,\left.{p}_{0}^{\dagger }{p}_{0}^{\dagger }{p}_{1}{p}_{1}+{p}_{0}^{\dagger }{p}_{0}^{\dagger }{p}_{-1}{p}_{-1}+2{p}_{1}^{\dagger }{p}_{0}^{\dagger }{p}_{0}{p}_{-1}+2{p}_{-1}^{\dagger }{p}_{0}^{\dagger }{p}_{0}{p}_{1}\right)\\ \ -\,g\left[2{p}_{1}^{\dagger }{p}_{-1}^{\dagger }{uv}+2{u}^{\dagger }{v}^{\dagger }{p}_{-1}{p}_{1}\right.\\ \ +\,{p}_{1}^{\dagger }{u}^{\dagger }{{up}}_{1}+{p}_{1}^{\dagger }{v}^{\dagger }{{vp}}_{1}+{p}_{-1}^{\dagger }{u}^{\dagger }{{up}}_{-1}+{p}_{-1}^{\dagger }{v}^{\dagger }{{vp}}_{-1}\\ \ +\,{p}_{1}^{\dagger }{p}_{1}^{\dagger }{uv}+{u}^{\dagger }{v}^{\dagger }{p}_{1}{p}_{1}+{p}_{-1}^{\dagger }{p}_{-1}^{\dagger }{uv}+{u}^{\dagger }{v}^{\dagger }{p}_{-1}{p}_{-1}\\ \ +\,\left.{p}_{1}^{\dagger }{u}^{\dagger }{{up}}_{-1}+{p}_{1}^{\dagger }{v}^{\dagger }{{vp}}_{-1}+{p}_{-1}^{\dagger }{u}^{\dagger }{{up}}_{1}+{p}_{-1}^{\dagger }{v}^{\dagger }{{vp}}_{1}\right]\\ \ -\,2h\left({p}_{0}^{\dagger }{u}^{\dagger }{{up}}_{0}-{p}_{0}^{\dagger }{p}_{0}^{\dagger }{uv}-{u}^{\dagger }{v}^{\dagger }{p}_{0}{p}_{0}+{v}^{\dagger }{p}_{0}^{\dagger }{p}_{0}v\right).\end{array}\end{eqnarray} \tag{ 39 }$$

In order to diagonalize the above Hamiltonian the following basis can be used

$$\begin{eqnarray}&&\left|{n}_{-1},{n}_{0},{n}_{1},n\right\rangle =\displaystyle \frac{{\left({p}_{-1}^{\dagger }\right)}^{{n}_{-1}}{\left({p}_{0}^{\dagger }\right)}^{{n}_{0}}{\left({p}_{1}^{\dagger }\right)}^{{n}_{1}}{\left({u}^{\dagger }\right)}^{n}{\left({v}^{\dagger }\right)}^{n}}{\sqrt{{n}_{-1}!{n}_{0}!{n}_{1}!{\left(n!\right)}^{2}}}| 0\rangle .\end{eqnarray} \tag{ 40 }$$

The states should fulfill the number constraint, which for the half-filled case implies

$$\begin{eqnarray}&&{n}_{-1}+{n}_{0}+{n}_{1}+2\,n=j.\end{eqnarray} \tag{ 41 }$$

Because n is an integer, $j-{n}_{0}-{n}_{1}-{n}_{-1}$ should be an even number. In addition, the states own a given parity, positive for even values of n₀ and negative for odd values.

The space spanned by the basis states (40) corresponds to a representation of the su(5) algebra. However, we are interested in the so(5) subspace that has a one to one correspondence with the physical states of the fermionic Agassi model. This physical subspace is characterized by a null value of the quadratic Casimir operator of so(5) that can be written as a boson pairing interaction

$$\begin{eqnarray}&&{P}^{\dagger }P=\left(2{p}_{1}^{\dagger }{p}_{-1}^{\dagger }-{p}_{0}^{\dagger }{p}_{0}^{\dagger }-2{u}^{\dagger }{v}^{\dagger }\right)\left(2{p}_{1}{p}_{-1}-{p}_{0}{p}_{0}-2{uv}\right),\end{eqnarray} \tag{ 42 }$$

meaning that these states have a generalized seniority zero. Alternatively, this condition can be used to select the physical states after a diagonalization in the enlarged basis (40) as those eigenstates that fulfill

$$\begin{eqnarray}&&P\left|{\rm{\Phi }}\right\rangle =\left(2{p}_{1}{p}_{-1}-{p}_{0}{p}_{0}-2{uv}\right)\left|{\rm{\Phi }}\right\rangle =0.\end{eqnarray} \tag{ 43 }$$

In this section, the extended Agassi model is studied by mean-field techniques. Two different approaches that lead to the same results are elaborated. First, a boson condensate depending on some variational parameters is used to obtain the energy surface associated to the Hamiltonian within its Schwinger image and, second, the HFB approach with the fermionic image of the Hamiltonian is worked out.

3.1. Mean field with a boson condensate in the Schwinger representation

The initial point of this variational method is a boson condensate, ${{\rm{\Gamma }}}^{\dagger }$, that is proposed to be a linear combination of the five bosons appearing in the Schwinger image and five variational parameters, namely, ${\eta }_{-1},{\eta }_{0},{\eta }_{1},\delta$, and γ

$$\begin{eqnarray}&&{{\rm{\Gamma }}}^{\dagger }=\displaystyle \frac{{\eta }_{-1}{p}_{-1}^{\dagger }+{\eta }_{0}{p}_{0}^{\dagger }+{\eta }_{1}{p}_{1}^{\dagger }+\delta {u}^{\dagger }+\gamma {v}^{\dagger }}{\sqrt{{\eta }_{-1}^{2}+{\eta }_{0}^{2}+{\eta }_{1}^{2}+{\delta }^{2}+{\gamma }^{2}}}.\end{eqnarray} \tag{ 44 }$$

Then, the trial wave function of the system ground state is proposed to be a condensate of j boson pairs as

$$\begin{eqnarray}&&| j,{\eta }_{-1},{\eta }_{0},{\eta }_{1},\delta ,\gamma \rangle =\displaystyle \frac{{\left({{\rm{\Gamma }}}^{\dagger }\right)}^{j}}{\sqrt{j!}}| 0\rangle .\end{eqnarray} \tag{ 45 }$$

Note that it is assumed that the number of fermion pairs is j and, therefore, the system is half filled.

Next step is to calculate the expectation value of the Hamiltonian (39) with the trial wave function (45)

$$\begin{eqnarray}&&\begin{array}{l}E(j,{\eta }_{-1},{\eta }_{0},{\eta }_{1},\delta ,\gamma )\\ \quad =\,\langle j,{\eta }_{-1},{\eta }_{0},{\eta }_{1},\delta ,\gamma | H| j,{\eta }_{-1},{\eta }_{0},{\eta }_{1},\delta ,\gamma \rangle ,\end{array}\end{eqnarray} \tag{ 46 }$$

to obtain the corresponding energy surface. In the $j\to \infty$ limit this energy surface reads as

$$\begin{eqnarray}\begin{array}{rcl}E & = & \varepsilon j\displaystyle \frac{{\eta }_{-1}^{2}+{\eta }_{1}^{2}}{{\eta }_{-1}^{2}+{\eta }_{0}^{2}+{\eta }_{1}^{2}+{\delta }^{2}+{\gamma }^{2}}\\ & & -{{gj}}^{2}\displaystyle \frac{({\eta }_{-1}^{2}+{\eta }_{1}^{2}){\left(\delta +\gamma \right)}^{2}}{{\left({\eta }_{-1}^{2}+{\eta }_{0}^{2}+{\eta }_{1}^{2}+{\delta }^{2}+{\gamma }^{2}\right)}^{2}}\\ & & -{{Vj}}^{2}\displaystyle \frac{2{\left({\eta }_{-1}+{\eta }_{1}\right)}^{2}{\eta }_{0}^{2}}{{\left({\eta }_{-1}^{2}+{\eta }_{0}^{2}+{\eta }_{1}^{2}+{\delta }^{2}+{\gamma }^{2}\right)}^{2}}\\ & & -2{{hj}}^{2}\displaystyle \frac{{\eta }_{0}^{2}{\left(\delta -\gamma \right)}^{2}}{{\left({\eta }_{-1}^{2}+{\eta }_{0}^{2}+{\eta }_{1}^{2}+{\delta }^{2}+{\gamma }^{2}\right)}^{2}},\end{array}\end{eqnarray} \tag{ 47 }$$

where the contributions from the two-body terms to the one-body terms have been neglected, since they vanish assuming that the two-body coefficients are scaled with an extra 1/j factor. As already commented before, the use of the parameters (11) that take into account the system size j is convenient. Using them, the energy functional is

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{E}{\varepsilon j} & = & \displaystyle \frac{{\eta }_{-1}^{2}+{\eta }_{1}^{2}}{{\eta }_{-1}^{2}+{\eta }_{0}^{2}+{\eta }_{1}^{2}+{\delta }^{2}+{\gamma }^{2}}\\ & & -\displaystyle \frac{{\rm{\Sigma }}}{2}\displaystyle \frac{({\eta }_{-1}^{2}+{\eta }_{1}^{2}){\left(\delta +\gamma \right)}^{2}}{{\left({\eta }_{-1}^{2}+{\eta }_{0}^{2}+{\eta }_{1}^{2}+{\delta }^{2}+{\gamma }^{2}\right)}^{2}}\\ & & -\chi \displaystyle \frac{{\left({\eta }_{-1}+{\eta }_{1}\right)}^{2}{\eta }_{0}^{2}}{{\left({\eta }_{-1}^{2}+{\eta }_{0}^{2}+{\eta }_{1}^{2}+{\delta }^{2}+{\gamma }^{2}\right)}^{2}}\\ & & -{\rm{\Lambda }}\displaystyle \frac{{\eta }_{0}^{2}{\left(\delta -\gamma \right)}^{2}}{{\left({\eta }_{-1}^{2}+{\eta }_{0}^{2}+{\eta }_{1}^{2}+{\delta }^{2}+{\gamma }^{2}\right)}^{2}}.\end{array}\end{eqnarray} \tag{ 48 }$$

Note that ε is an overall energy constant, that in practice can be considered as ε = 1. Please note, that in mean-field a large system size j is assumed, consequently, the term −1 in the denominator $(2j-1)$ of equations (11) is negligible.

In principle, the energy surface depends on five variational parameters, but in practice, this number is smaller. First, for positive values of ε, which corresponds to the physical cases, ${\eta }_{-1}$ can be fixed to ${\eta }_{-1}=1$. Moreover, because of the shape of the function (48), either δ = γ or δ = −γ for all possible minima. This fact has very important consequences on the way the energy function behaves with respect to the control parameters. As a matter of fact, parameters Σ and Λ never affect the value of the energy in its minimum at the same time. In general, for small enough values of Λ the energy only depends on χ and Σ, while for large enough values of Λ the energy only depends on χ and Λ. Therefore, the system behaves as if two separated energy surfaces were competing. One has to look for the minima of the two competing surfaces by minimizing both energy surfaces with respect to the variational parameters and look for the global minimum. The detailed study of the energy surfaces in the Schwinger representation can be carried out both numerical and analytically, however, the results are fully equivalent to those obtained with the HFB approach. Since the HFB derivation is more straightforward, only this approach is worked out in detail in the next subsection.

3.2. The HFB approach

Our mean-field treatment for the extended Agassi model follows closely the approach presented in [10, 19]. First, a Hartree–Fock transformation is applied, and then it is complemented with a Bogoliubov one. Different energy surfaces, called A and B, are obtained with two alternative phase selections for the Bogoliubov transformation (see [19] for details). The variational mean-field state is defined via the HFB formalism in terms of the parameters φ and β. Then, the expectation value of the Hamiltonian (12) produces an energy surface

$$\begin{eqnarray}&&E(\varphi ,\beta )=\displaystyle \frac{\langle {\rm{HFB}}(\varphi ,\beta )| H| {\rm{HFB}}(\varphi ,\beta )\rangle }{\langle {\rm{HFB}}(\varphi ,\beta )| {\rm{HFB}}(\varphi ,\beta )\rangle }.\end{eqnarray} \tag{ 49 }$$

The minimization with respect to φ and β of the energy functional E(φ, β) allows one to find the extrema. Their characterization (maximum, minimum or saddle point) is done by constructing the Hessian matrix. The complete analysis for the two alternative phase selections is done in [19] and a short summary is presented here for completeness. Depending on the phase selection two energy surfaces, A and B, are obtained with scaled energies

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{A}}{j\varepsilon }=-\cos \varphi \cos \beta -\displaystyle \frac{{\rm{\Sigma }}}{2}{\sin }^{2}\beta -\displaystyle \frac{\chi }{2}{\sin }^{2}\varphi {\cos }^{2}\beta \end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{B}}{j\varepsilon }=-\cos \varphi \cos \beta -{\rm{\Lambda }}{\sin }^{2}\beta {\sin }^{2}\varphi -\displaystyle \frac{\chi }{2}{\sin }^{2}\varphi {\cos }^{2}\beta .\end{eqnarray} \tag{ 51 }$$

In order to determine the ground state of the system one has to look for the absolute minima of the two surfaces. Because of the existence of two surfaces there are regions in which different phases coexist. The detailed discussion on these points was presented in [19]. Five different phases are established in the phase diagram (figure 2)

• Symmetric or spherical solution: it is the region with χ < 1, Σ < 1 and Λ < 1.
• HF deformed solution: it is the region χ > 1, χ > Σ and ${\rm{\Lambda }}\lt \tfrac{1+{\chi }^{2}}{2\chi }$.
• BCS deformed solutions: it is the region Σ > 1, Σ > χ and ${\rm{\Lambda }}\lt \tfrac{1+{{\rm{\Sigma }}}^{2}}{2{\rm{\Sigma }}}$.
• Combined HF-BCS deformed solution: it is the region with ${\rm{\Lambda }}\gt 1,{\rm{f}}{\rm{o}}{\rm{r}}\,{\rm{\Sigma }}\lt 1\,{\rm{a}}{\rm{n}}{\rm{d}}\,\chi \lt 1,\,{\rm{\Lambda }}\gt {\textstyle {\textstyle \tfrac{1+{{\rm{\Sigma }}}^{2}}{2{\rm{\Sigma }}}}}\,{\rm{f}}{\rm{o}}{\rm{r}}\,{\rm{\Sigma }}\,\gt \chi \,{\rm{a}}{\rm{n}}{\rm{d}}\,{\rm{\Sigma }}\gt 1$ and ${\rm{\Lambda }}\gt {\textstyle \tfrac{1+{\chi }^{2}}{2\chi }}\,{\rm{f}}{\rm{o}}{\rm{r}}\,{\rm{\Sigma }}\lt \chi \,{\rm{a}}{\rm{n}}{\rm{d}}\,\chi \gt 1$.
• Closed valley solution: it is the plane χ = Σ with ${\rm{\Lambda }}\lt \tfrac{1+{\chi }^{2}}{2\chi }$.

Zoom In Zoom Out Reset image size

In figure 2 the different phases are plotted in the phase diagram of the extended Agassi. The first one, represented by a red sphere, corresponds to the spherical solution (φ = 0, β = 0) and it is limited to the area with χ < 1, Σ < 1, and Λ < 1. The value of the energy in this region is E/(jε) = −1. The limits of the region correspond to QPT surfaces, i.e. second order QPT for the vertical planes with χ = 1 and Σ = 1 (in red) and first order QPT for the plane Λ = 1 (in green). The second area, represented by a blue oval, corresponds to the HF deformed solution $(\cos \varphi =\tfrac{1}{\chi },\beta =0)$ and it is limited by the surfaces χ = 1 (in red), which implies a second order QPT, χ = Σ (in blue), which supposes a first order QPT, and ${\rm{\Lambda }}=\tfrac{1+{\chi }^{2}}{2\chi }$ (in green) that also implies a first order QPT. The value of the energy in this case is $E/(j\varepsilon )=-\tfrac{1+{\chi }^{2}}{2\chi }$. The third area, represented by a black oval, corresponds to the BCS deformed solution ($\varphi =0,\cos \beta =\tfrac{1}{{\rm{\Sigma }}}$) and the value of the energy is $E/(j\varepsilon )=-\tfrac{1+{{\rm{\Sigma }}}^{2}}{2{\rm{\Sigma }}}$. The region is limited by the plane Σ = 1 (in red), corresponding to a second order QPT, and the green surface ${\rm{\Lambda }}=\tfrac{1+{{\rm{\Sigma }}}^{2}}{2{\rm{\Sigma }}}$ and the plane χ = Σ (in blue) that correspond to first order QPTs. Finally, the fourth area, represented by the green crossed ovals, corresponds to the HF-BCS deformed solution $(\varphi =\tfrac{\pi }{2},\beta =\tfrac{\pi }{2})$ with energy E/(jε) = −Λ. This region is limited by the green surface which implies a first order QPT. Also note that for the surface χ = Σ the closed valley solution ($\cos \varphi \cos \beta =1/\chi$), represented by a thick ellipse, is also valid.

Take into account that the phases represented in the phase diagram are the ones corresponding to the deepest absolute minimum of the mean-field energy (including both surfaces A and B). However, in each region several phases can coexist, up to three (see [19] for details). In addition, in the line χ = Σ with ${\rm{\Lambda }}=\tfrac{1+{\chi }^{2}}{2\chi }$, four phases, HF, BCS, HF-BCS and the closed valley solutions, are degenerated. Finally, in the single point, χ = Σ = Λ = 1, the five solutions for the system are degenerated. This provides with a richer phase diagram than for other studied complex systems, such as the two-fluid LGM model [21], the proton–neutron IBM [22], or for Hamiltonians with up two three-body interactions [23].

A very convenient test of the kind of mean-field calculations presented so far is its comparison with the exact results. However, any exact diagonalization will be done for a finite size system, therefore finite size discrepancies between the exact and the mean-field results will always appear. These, of course, will be reduced as the size of the system increases, although they will always reach a limit due to the computational limitations.

A possible comparison between the exact and the mean-field results is a 1/j analysis. This will provide the value of the considered quantity in the thermodynamic limit and, therefore, should coincide with the mean-field value. Here, the approach presented in [24] is closely followed, where the ground-state energy is written as

$$\begin{eqnarray}&&\displaystyle \frac{{ \mathcal E }}{j}=a+\displaystyle \frac{b}{j}+\displaystyle \frac{c}{{j}^{2}}+\displaystyle \frac{d}{{j}^{3}}+O\left(\displaystyle \frac{1}{{j}^{4}}\right),\end{eqnarray} \tag{ 52 }$$

where it is assumed that the two-body coefficients are rescaled with the size of the system, i.e. they are multiplied by 1/j. The value of the energy in the thermodynamic limit will correspond, therefore, to the intercept, a.

Note that expression (52) is, in principle, only valid well apart from the critical point where the observables behave as j^α in a second-order phase transition and the value of α depends on the system under study.

This expansion is justified on the light of a shift transformation ${a}^{\dagger }={b}^{\dagger }+\lambda \sqrt{j}$, because for $\tfrac{{{ \mathcal O }}_{1}}{j}$, where ${{ \mathcal O }}_{1}$ is a one-body term, the highest order in j will be j⁰, while the lowest j⁻¹. In the case of a two-body term, ${{ \mathcal O }}_{2}$, appropriately scaled, $\tfrac{{{ \mathcal O }}_{2}}{{j}^{2}}$, the highest order in j is, again, j⁰ while the lowest j⁻². In our calculations, up to quadratic terms in the 1/j expansion of the studied observable will be considered.

In order to obtain the thermodynamic limit of the ground-state energy, its representation as a function of 1/j is plotted from j = 10 to j = 300 in steps of j = 10. The value of the intercept, obtained after performing a quadratic least squares fit, will correspond to the thermodynamic limit of the quantity under study. In principle, this analysis can be extended to any other excited state, although in this case it cannot be compared with the mean-field value.

The values of the effective order parameters are also compared using the same procedure. These are defined in terms of the expectation values of the following operators for the ground state

$$\begin{eqnarray}&&{{OP}}_{{J}^{2}}=\displaystyle \frac{\langle {\left({J}^{+}\right)}^{2}\rangle +\langle {\left({J}^{-}\right)}^{2}\rangle }{2{j}^{2}},\end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray}&&{{OP}}_{{A}_{0}^{2}}=\displaystyle \frac{\langle {A}_{0}^{+}{A}_{0}\rangle }{{j}^{2}},\end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray}&&{{OP}}_{{A}_{1}^{2}}=\displaystyle \frac{\langle {A}_{1}^{+}{A}_{1}\rangle +\langle {A}_{-1}^{+}{A}_{-1}\rangle }{2{j}^{2}}.\end{eqnarray} \tag{ 55 }$$

We first study in figure 3 the case χ = 0.5, Σ = 0.5, Λ = 0 (note that in all numerical calculations it is assumed ε = 1), which corresponds to a point in the phase diagram with φ = 0, β = 0, and E_gs/j = −1. The value of the energy in the thermodynamic limit is E_gs/j = −1.000 0034 ± 5 × 10⁻⁷ to be compared with −1. In the case of the order parameters, their obtained values in the thermodynamic limit are (−8.9 ± 1.4) × 10⁻⁷, (−3.1 ± 0.4) × 10⁻⁷ and (5.0 ± 0.9) × 10⁻⁶ for ${{OP}}_{{A}_{0}^{2}},{{OP}}_{{A}_{1}^{2}}$, and ${{OP}}_{{J}^{2}}$, respectively. These have to be compared with the mean-field values that are zero for the three cases. This supposes a difference between the mean-field and the extracted thermodynamic limit values of the order of 10⁻⁶ or smaller for all the magnitudes under study.

Zoom In Zoom Out Reset image size

In figure 4 the case χ = 0.5, Σ = 0.5, Λ = 2 is considered. It corresponds to a point in the phase diagram with φ = π/2, β = π/2, and E_gs = −Λ. The value of the energy in the thermodynamic limit is E_gs/j = −2.000 016 ± 3 × 10⁻⁶ to be compared with the mean field value −2. In the case of the order parameters, their values in the thermodynamic limit result to be 0.999 9999 ± 2 × 10⁻⁷, (−2.2 ± 0.4) × 10⁻⁶, and (2.0 ± 0.5) × 10⁻⁷ for ${{OP}}_{{A}_{0}^{2}},{{OP}}_{{A}_{1}^{2}}$, and ${{OP}}_{{J}^{2}}$, respectively. These have to be compared with the mean-field values that are 1, 0, and 0, respectively. This supposes a difference between the mean-field and the values in the thermodynamic limit smaller than 10⁻⁷ for the energy and of the order of 10⁻⁶ for the order parameters.

Zoom In Zoom Out Reset image size

In figure 5 the case χ = 0.5, Σ = 1.5, Λ = 0 is selected. This corresponds to a point in the phase diagram with $\varphi =0,\beta =\arccos (1/{\rm{\Sigma }})$, and ${E}_{{gs}}=-\tfrac{{{\rm{\Sigma }}}^{2}+1}{2{\rm{\Sigma }}}$. The value of the energy obtained for the thermodynamic limit is E_gs/j = −1.083 354 ± 3 × 10⁻⁶ to be compared with the mean-field value −1.083 333 (−13/12). In the case of the order parameters their values in the thermodynamic limit read as (−2.1 ± 0.3) × 10⁻⁵, 0.277 739 ± 4 × 10⁻⁶, and (1.6 ± 0.3) × 10⁻⁵ for ${{OP}}_{{A}_{0}^{2}},{{OP}}_{{A}_{1}^{2}}$, and ${{OP}}_{{J}^{2}}$, respectively. These results are to be compared with the mean-field values that are 0, 5/18 ≈ 0.277 778, and 0 respectively. This supposes, once more, a difference between the mean-field and the extracted thermodynamic limit of the order of 10⁻⁵ for the energy and for the order parameters.

Zoom In Zoom Out Reset image size

Once the technique has been presented, the comparison between the mean-field and the exact results in the large-j limit for selected trajectories is presented now. For specific values of the control parameters, a calculation with j ranging from 10 to 200 is performed extracting from a least squares fit to a function $a+\tfrac{b}{j}+\tfrac{c}{{j}^{2}}$ the value of the intercept, a, and its corresponding error. This value corresponds to the thermodynamic limit of the quantity under study. Same quantities as in the preceding discussion are studied, but to better appreciate the agreement between both, the mean-field and the exact results in the thermodynamic limit, the difference between both values is plotted in the following figures, together with the value of the statistical error extracted from the least squares fit. The selected trajectories move between two different phases crossing a QPT surface either of first or of second order. These trajectories are those studied in [19], although in the later case the exact calculations were performed for j = 50.

As a first case, the trajectory Σ = 0.5, Λ = 0 as a function of χ is presented in figure 6. As shown in panel (e), this trajectory crosses a second order QPT at χ = 1. The values for the difference between the large-j thermodynamic limit and the mean-field values for the ground state energy (panel a), and the order parameters associated to the expection values of J² (panel b), ${A}_{0}^{2}$(panel c), and ${A}_{1}^{2}$ (panel d), including their errors bars, are plotted. A remarkable agreement between the values in the thermodynamic limit and the mean-field ones is observed, although the agreement is clearly better in the spherical phase than in the HF deformed one. Moreover, the presented approach, as already mentioned, is not well suited for the region around the QPT, that is around χ ≈ 1, where the agreement worsen. Note that some points around the QPT may lay out of the plotted scale.

Zoom In Zoom Out Reset image size

The second trajectory studied is χ = 1.5, Λ = 0.5, as a function of Σ and the results are displayed in figure 7. As shown in panel (e), this trajectory crosses a first order QPT at Σ = 1.5. As in the preceding case, a remarkable agreement between the values in the thermodynamic limit and the mean-field ones is observed, although here the agreement is better for the BCS-deformed phase than for the HF-deformed region. Again, larger discrepancies, as expected, are observed in a small region close to the critical point at Σ = 1.5. Note that some points around the QPT may lay out of the plotted scale.

Zoom In Zoom Out Reset image size

The last trajectory analyzed is χ = 1.5, Σ = 2, as a function of Λ and the results are plotted in figure 8. As shown in panel (e), this trajectory crosses a first order QPT at Λ = 5/4. In this case the agreement is also noticeable. The difference in all studied magnitudes between the mean-field results and those obtained in the large-j thermodynamic limit is basically zero except in a small region around the critical point Λ = 5/4. Note that some points around the QPT may lay out of the plotted scale.

Zoom In Zoom Out Reset image size

An extended version of the Agassi model that includes an extra ${A}_{0}^{\dagger }{A}_{0}$ contribution (a non-standard pairing contribution) in the Hamiltonian has been reviewed. Therefore, the model depends on four free control parameters, ε, g, V, and h. However, a non-vanishing value for ε has always been considered. Hence, the number of effective free parameters is three: $V=\tfrac{\varepsilon \chi }{2j-1},g=\tfrac{\varepsilon {\rm{\Sigma }}}{2j-1}$, and $h=\tfrac{\varepsilon {\rm{\Lambda }}}{2j-1}$. The underlying algebra of the model, o(5), has been presented, and both the fermionic and the Schwinger boson representation of the model have been worked out. A mean-field approach for the extended Agassi model through a HFB mean-field approach is presented. In addition, the alternative use of a condensate of bosons in the Schwinger representation is discussed. The corresponding energy surfaces are obtained and the phase diagram of the model is established. The phase diagram presents four different regions: spherical, HF deformed, BCS deformed, and combined HF-BCS deformed phases. Moreover, there is a surface in which a special solution, called closed valley, exists. These phases are separated by several surfaces corresponding either to first or second order QPTs. In addition, there is a line in which four phases coexist (HF deformed, BCS deformed, combined HF-BCS deformed, and close valley deformed minimum) and are degenerated. There exists a single point χ = Σ = Λ = 1 in which the five phases are degenerated (spherical, HF deformed, BCS deformed, combined HF-BCS deformed, and close valley deformed minimum).

Finally, the exact results in the large-j limit of the model are extracted by using a quadratic 1/j expanssion and compared with the mean-field results. First, the exact large-j limit is obtained for specific values of the control parameters performing a 1/j analysis with j ranging from j = 10 to j = 300. Second, we have applied the technique for selected trajectories that cross QPT regions. In all cases, the large-j analysis provides very accurate values for the energy and the order parameters, and the correspondence with the mean-field values is noticeable. Therefore, the proposed technique seems to be suitable to calculate the thermodynamic limit of any other observables.

This work has been supported by the Spanish Ministerio de Economía y Competitividad and the European regional development fund (FEDER) under Projects No. FIS2017-88410-88410-P, FIS2014-53448-C2-2-P and FIS2015-63770-P, and by Consejería de Economía, Innovación, Ciencia y Empleo de la Junta de Andalucía (Spain) under Group FQM-160 and FQM-370.