Reproduction of exact solutions of Lipkin model by nonlinear random-phase approximation

TL;DR

This paper demonstrates that a nonlinear extension of the random-phase approximation (RPA) can exactly reproduce solutions of the Lipkin model, suggesting its potential for accurate many-body problem calculations.

Contribution

The authors show that nonlinear RPA reproduces exact solutions of the Lipkin model, establishing its equivalence to the Schrödinger equation for this system.

Findings

Exact solutions for N=2 analytically derived

Numerical validation for N=20

Nonlinear RPA matches the Schrödinger equation solutions

Abstract

It is shown that the random-phase approximation (RPA) method with its nonlinear generalization, which was previously considered as approximation, reproduces the exact solutions of the Lipkin model. The nonlinear RPA is based on an equation nonlinear on eigenvectors and includes many-particle-many-hole components in the creation operator of the excited states. We demonstrate the exact character of solutions analytically for the particle number $N$ = 2 and, numerically, for $N$ = 20. This finding indicates that the nonlinear RPA is equivalent to the exact Schr\"{o}dinger equation, which opens up new possibilities for realistic calculations in many-body problems.

Tables

Table 1. Table 1: Eigenstates, wave functions, total energies, excitation energies, and phonon-creation operators obtained for N = 2 𝑁 2 N=2 by the nonlinear higher RPA. The wave functions and energies are identical to those obtained by the exact shell model. E 10 o superscript subscript 𝐸 10 𝑜 E_{10}^{o} is given by Eq. ( 91 ). Each of the odd- and even- order subspaces has one excited state.

Eigenstate	Wave function	Total energy
Ground	$|{\Psi }_0⟩=\frac{V}{\sqrt{2E_{10}^o\left(E_{10}^o-\epsilon \right)}}\left(1-\frac{E_{10}^o-\epsilon }{2V}{J}_{+}^2\right)|{\psi }_0⟩$	$-E_{10}^o$
Odd-order excited	$Q_1^{o†}|{\Psi }_0⟩=\frac{1}{\sqrt{2}}{J}_{+}|{\psi }_0⟩$	0
Even-order excited	$Q_1^{e†}|{\Psi }_0⟩=\frac{V}{\sqrt{2E_{10}^o\left(E_{10}^o+\epsilon \right)}}\left(1+\frac{E_{10}^o+\epsilon }{2V}{J}_{+}^2\right)|{\psi }_0⟩$	$E_{10}^o$

Table 2. Table 2: Squared norm of components of lower and higher order with respect to J + subscript 𝐽 J_{+} (the symmetry-conserving basis) or J ~ + subscript ~ 𝐽 \tilde{J}_{+} (the symmetry-breaking basis) of the exact first-excited state with V = − 0.2 𝑉 0.2 V=-0.2 . The lower order in the symmetry conserving basis consists of the first and third (the higher is the fifth and seventh), and that in the symmetry-breaking basis is from the zeroth to the fourth (the higher is from the fifth to the eighth).

Basis	Squared norm of components
	Lower order	Higher order
Symmetry conserving	0.943	0.057
Symmetry breaking	0.716	0.284

Equations

$$H=\epsilon J_z+\frac{V}{2}\left(J_{+}^2+J_{-}^2\right),$$

$$J_z$$

$$J_{+}$$

$$[J_z,J_{+}]=J_{+},[J_z,J_{-}]=-J_{-},[J_{+},J_{-}]=2J_z.$$

$$J_z{J}_{+}^i|{\psi }_0⟩=f(i)J_{+}^i|{\psi }_0⟩,f(i)=i-\frac{N}{2},$$

$$J_{-}{J}_{+}^i|{\psi }_0⟩=g(i)J_{+}^{i-1}|{\psi }_0⟩,g(i)=-i^2+(N+1)i,$$

$$f(i)=g(i)=0,(i<0,i>2n),$$

$$G(i)$$

$$[H,Q_k^{†}]|{\Psi }_0⟩=E_{k0}{Q}_k^{†}|{\Psi }_0⟩.$$

$$Q_k|{\Psi }_0⟩=0,$$

$$Q_k^{o†}$$

$$Q_k^{e†}$$

$${\mathcal{J}}_{\pm }^i=J_{\pm }^i/\sqrt{G(i)}.$$

$$|{\Psi }_0⟩=\sum _{i=0}^n{\beta }_{2i}{\mathcal{J}}_{+}^{2i}|{\psi }_0⟩.$$

$$\displaystyle\left(\begin{array}[]{cc}A^{o}&B^{o}\\ B^{o}&A^{o}\end{array}\right)\left(\begin{array}[]{c}\bm{X}^{o}_{k}\\ \bm{Y}^{o}_{k}\end{array}\right)=E^{o}_{k0}\left(\begin{array}[]{cc}U^{o}&O\\ O&-U^{o}\end{array}\right)\left(\begin{array}[]{c}\bm{X}^{o}_{k}\\ \bm{Y}^{o}_{k}\end{array}\right),$$

$${\mathbf{X}}_k^o$$

$${\mathbf{Y}}_k^o$$

$$[A,B,C]=\frac{1}{2}[[A,B],C]+\frac{1}{2}[A,[B,C]],$$

$$A_{ij}^o$$

$$B_{ij}^o$$

$$U_{ij}^o$$

$$A_{ij}^o=A_{ji}^o,B_{ij}^o=B_{ji}^o.$$

$$U_{ij}^o=U_{ji}^o.$$

$$⟨{\Psi }_0|[{\mathcal{J}}_{-}^k,H,{\mathcal{J}}_{+}^l]|{\Psi }_0⟩$$

$$\displaystyle=\frac{\varepsilon}{2}\sum_{{\scriptsize\begin{array}[]{c}i\!=\!0\\ (0\!\leq\!2i\!+\!k\!-\!l\!\leq\!2n)\end{array}}}^{n}\frac{\beta_{2i}\beta_{2i+k-l}}{\sqrt{G(k)G(l)G(2i)G(2i+k-l)}}\Big{[}\{2f(2i+k)-f(2i)$$

$$\displaystyle-f(2i+k-l)\}G(2i+k)+\big{\{}-f(2i-l+k)+2f(2i-l)-f(2i)\big{\}}$$

$$\displaystyle\times g(2i)g(2i-1)\cdots g(2i-l+1)G(2i-l+k)\Big{]}$$

$$\displaystyle+\frac{V}{4}\sum_{{\scriptsize\begin{array}[]{c}i\!=\!0\\ (0\!\leq\!2i\!+\!k\!-\!2\!-\!l\!\leq\!2n)\end{array}}}^{n}\frac{\beta_{2i}\beta_{2i+k-2-l}}{\sqrt{G(k)G(l)G(2i)G(2i+k-2-l)}}\big{\{}G(2i+k)$$

$$-g(2i)g(2i-1)G(2i+k-2)-g(2i)g(2i-1)\cdots g(2i-l+1)G(2i-l+k)$$

$$\displaystyle+g(2i)g(2i-1)\cdots g(2i-l-1)G(2i-l-2+k)\big{\}}$$

$$\displaystyle+\frac{V}{4}\sum_{{\scriptsize\begin{array}[]{c}i\!=\!0\\ (0\!\leq\!2i\!+\!k\!+\!2\!-\!l\leq\!2n)\end{array}}}^{n}\frac{\beta_{2i}\beta_{2i+k+2-l}}{\sqrt{G(k)G(l)G(2i)G(2i+k+2-l)}}\big{\{}G(2i+k+2)$$

$$+g(2i)g(2i-1)\cdots g(2i-l+1)G(2i-l+2+k)$$

$$-g(2i+k)g(2i+k-1)\cdots g(2i+k-l+1)G(2i+k-l+2)$$

$$\displaystyle-g(2i+2)g(2i+1)\cdots g(2i+3-l)G(2i+2-l+k)\big{\}},$$

$$⟨{\Psi }_0|[{\mathcal{J}}_{-}^k,H,{\mathcal{J}}_{-}^l]|{\Psi }_0⟩$$

$$⟨{\Psi }_0|[{\mathcal{J}}_{-}^k,{\mathcal{J}}_{+}^l]|{\Psi }_0⟩=$$

$$\displaystyle-G(2i)g(2i-l+k)g(2i-l+k-1)\cdots g(2i-l+1)\big{\}}.$$

$$\displaystyle\left(\begin{array}[]{ccc}L^{k}_{11}&\cdots&L^{k}_{1n}\\ \vdots&&\vdots\\ L^{k}_{n1}&\cdots&L^{k}_{nn}\end{array}\right)\left(\begin{array}[]{c}\beta_{2}/\beta_{0}\\ \vdots\\ \beta_{2n}/\beta_{0}\end{array}\right)=\left(\begin{array}[]{c}-Y^{k}_{1}\\ \vdots\\ -Y^{k}_{2n-1}\end{array}\right),$$

$$L_{ij}^k=X_{2j-2i+1}^k\sqrt{\frac{G(i)G(2j)}{G(2j-2i+1)}}\frac{1}{G(2i-1)},$$

Full text

Random-phase approximation; Lipkin model; many-body problem.

I Introduction

The random-phase approximation (RPA) Bohm and Pines (1953); Sawada et al. (1957); Brout (1957); Wentzel (1957); Hubbard (1958); Gell-Mann and Brueckner (1957) and its quasiparticle generalization (QRPA) Marumori (1960); Arvieu and Vénéroni (1960) have been, for a long time, very important theoretical many-body methods in quantum chemistry, condensed matter physics and nuclear physics. Hence, it is natural to expect that an extension of the RPA will give a new more powerful method. Areas in need of more accurate methods of calculation include neutrino physics in connection with the search for the Majorana neutrino mass, constraints on which depend substantially on the nuclear matrix elements of neutrinoless double-$\beta$ decay Engel and Menéndez (2017).

The RPA approach in its current formulation including many refinements is an approximation which cuts off the excitations at the one-particle-one-hole (1p-1h) level. The extension of the RPA to include also the 2p-2h excitations, so called the second RPA, has been investigated and used by many authors Werthamer and Suhl (1962); Fano and Sawicki (1962); Sawicki (1962); Da Providencia (1965); Rowe (1968); Yannouleas et al. (1983); Takayanagi et al. (1988); Papakonstantinou and Roth (2010); Gambacurta et al. (2016); Takahara et al. (2004); Schuck and Tohyama (2016). Currently the frontier of this approach is the self-consistent second RPA Schuck and Tohyama (2016), which is still used only for schematic models (their equation is solved only approximately). The aim of this work is the inclusion in the RPA framework of states with arbitrary order of particle-hole excitations from the ground state and the nonlinearity of the eigenequation Rowe (1968); Catara et al. (1994, 1996); Takahara et al. (2004); Schuck and Tohyama (2016), and our equations are solved exactly. It turns out that such extension reproduces the exact solutions of the Lipkin model Lipkin et al. (1965). In what follows, the extended RPA is referred to as the nonlinear higher RPA. We call the novel creation operator of the excited state the phonon operator for simplicity. Note, however, that the boson commutation relation is not assumed.

In Sec. II we show the equations for applying our method to the Lipkin model, and the analytical (particle number of 2) and numerical (larger particle numbers) solutions are presented; these are the exact solutions. The result of the truncation approximation is also shown and compared with the shell model. In Sec. III the formulation using the symmetry-breaking basis for the large interaction strength is discussed and numerically investigated. Section IV is devoted to summary.

II Application to Lipkin model

II.1 Formulation

The single-particle space of the Lipkin model consists of two fermion levels, each of which has an $N$-fold degeneracy. The upper (lower) level has the energy of $\varepsilon/2$ ($-\varepsilon/2$). We assume, without loss of generality, that $N$ is even and equal to the particle number of the system. A parameter of $n=N/2$ is often used in this paper. The system state in which all particles are in the lower level is denoted by $|\psi_{0}\rangle$. The Hamiltonian of the Lipkin model is given by

[TABLE]

[TABLE]

The creation and annihilation operators of the fermion are denoted by $a^{\dagger}_{im}$ and $a_{im}$, ($i=0$: lower level, $i=1$: upper level). Index $m$ distinguishes the degenerated states. $J_{+}$, $J_{-}$, and $J_{z}$ satisfy the following commutation relations:

[TABLE]

$V$ is the strength of the interaction. Our purpose is to test our new method, therefore we are interested in the space affected by the interaction Lipkin et al. (1965). For this reason, the space relevant to us is spanned by state vectors $J_{+}^{i}|\psi_{0}\rangle$, $(i=0,\cdots,2n)$. This space splits into two subspaces. One is the odd-order subspace with respect to $J_{+}$, and another is the even-order subspace ($|\psi_{0}\rangle$ included).

The following formulae can be derived from the commutation relations (4):

[TABLE]

for $0\leq i\leq 2n$. We extend the definition of $f(i)$ and $g(i)$ to

[TABLE]

and introduce a function

[TABLE]

The phonon-creation operator $Q_{k}^{\dagger}$ ($k$ denotes an excited state) and the excitation energy $E_{k0}$ are determined by the equation of motion

[TABLE]

The ground state $|\Psi_{0}\rangle$ is determined by

[TABLE]

which we call the vacuum condition. $Q_{k}^{\dagger}|\Psi_{0}\rangle$ is the excited state, and its orthogonality to the ground state is guaranteed by the vacuum condition. The general framework is defined by Eqs. (17) and (18).

$Q_{k}^{\dagger}$ for the Lipkin model is set to

[TABLE]

with

[TABLE]

Here, $c_{k}$ is a c-number determined together with $X_{i}^{k}$’s and $Y_{i}^{k}$’s by solving the equations. The denominator in ${\cal J}_{\pm}^{i}$ is introduced to avoid the overflow in the numerical calculation of the Hamiltonian matrix elements for large $n$. The ground state can be written

[TABLE]

Equation (17) for the odd-order subspace can be cast into the matrix-vector form

[TABLE]

[TABLE]

$A^{o}$, $B^{o}$, and $U^{o}$ are $n\times n$ matrices, and the suffix $T$ stands for transpose. The equation for the even-order subspace can be written analogously. With the abbreviation of the symmetric double commutator

[TABLE]

the matrix elements in Eq. (31) are defined by

[TABLE]

for $i,j=1,\cdots,n$. The use of $|\Psi_{0}\rangle$ is the ultimate extension of the renormalized RPA Rowe (1968). The symmetric double commutator including $H$ is used for guaranteeing the symmetry of the Hamiltonian matrix:

[TABLE]

We also have

[TABLE]

The equations of the matrix elements and symmetry relations for the even-order subspace can be obtained analogously. Matrix elements of the even-order subspace are labeled with suffix $e$ (see below). The explicit equations of $A_{ij}^{o}$ and $A_{ij}^{e}$ can be obtained from

[TABLE]

for $k,l=1,\cdots,2n$. This equation is derived by making use of Eqs. (5)$-$(16). For $B_{ij}^{o}$ and $B_{ij}^{e}$, we can use

[TABLE]

and $U_{ij}^{o}$ and $U_{ij}^{e}$ are calculated by

[TABLE]

The vacuum condition (18) for $k$ belonging to the odd-order subspace yields

[TABLE]

i) $i=1,\cdots,n;\ j=i,\cdots,n$, (upper triangle including the diagonal line)

[TABLE]

ii) $i=2,\cdots,n;\ j=1,\cdots,i-1$, (lower triangle)

[TABLE]

If $X^{k}_{2i-1}$’s and $Y^{k}_{2i-1}$’s are given, Eq. (67) seems at first glance to indicate that $\beta_{2i}/\beta_{0}$’s depend on $k$. Actually, the solution is independent of $k$. An analysis related to this property is shown below using the numerical result.

For $k$ of the even-order subspace, the vacuum condition (18) yields the following equations (note that there are $n+1$ equations):

[TABLE]

Apparently, any of these equations determines $c_{k}$, if $\beta_{2i}$’s, $X^{k}_{2i}$’s, and $Y^{k}_{2i}$’s are given. It is confirmed numerically below that $c_{k}$ is independent of the choice of the equation. For the treatment of the vacuum condition by previous papers, see, e.g., Refs. Balian and Brezin (1969); Dukelsky and Schuck (1990).

The number of excited states is the same as that in the shell model (diagonalization of $H$ matrix represented by an orthonormal basis) as the phonon operators are constructed from the operators creating the orthonormal basis ($J_{+}^{i}$) and their hermite conjugates.

II.2 Analytical result

The equations for the odd-order subspace with $N=2$ are considered analytically. Since there is only one excited state, we write

[TABLE]

The vacuum condition (18) gives

[TABLE]

and the eigenequation (31) reads

[TABLE]

The matrix elements are given by

[TABLE]

where Eq. (74) is used. From the above equations the following equation for $x\equiv\alpha_{-}/\alpha_{+}$ is obtained:

[TABLE]

This algebraic equation has six solutions:

[TABLE]

We choose the physical solution

[TABLE]

satisfying

[TABLE]

Equations (74) and (89) give the ratio of the components of the exact ground state. The excitation energy is obtained

[TABLE]

which is identical to the exact one. The wavefunction $Q_{o}^{\dagger}|\Psi_{0}\rangle$ is equal to $\frac{1}{\sqrt{2}}J_{+}|\psi_{0}\rangle$. It is possible to reproduce the exact solutions of the even-order subspace without high-order algebraic equation by using $\beta_{0}$ and $\beta_{2}$ obtained by the odd-order subspace calculation. The analytical solutions for $N=2$ are summarized in Table 1. It is possible to confirm the explicit equations of $Q_{1}^{o\dagger}|\Psi_{0}\rangle$ and $Q_{1}^{e\dagger}|\Psi_{0}\rangle$ as well as $Q_{1}^{o}|\Psi_{0}\rangle=0$ and $Q_{1}^{e}|\Psi_{0}\rangle=0$ using the equations of this table.

If the equations of the self-consistent second RPA Takahara et al. (2004); Jemaï et al. (2013); Schuck and Tohyama (2016) with $c_{k}$ are solved exactly for $N=2$, our result should be obtained. The vacuum condition for the Lipkin model is solved in Ref. Dukelsky and Schuck (1990). The self-consistent RPA Jemaï et al. (2013); Delion et al. (2005) reproduces the ground and first excited (the odd-order subspace) states for $N=2$. The studies of Refs. Dukelsky and Schuck (1990); Delion et al. (2005) do not obtain the excited states in the even-order subspace because the phonon operators in these studies have only the 1p-1h components. That of Ref. Jemaï et al. (2013) includes the 2p-2h components in the phonon operator but does not obtain the exact even-order excited state for $N=2$ because $c_{k}$ is not used. The excitation energy of the RPA (with the 1p-1h phonon operator and the $A^{o}$, $B^{o}$, and $U^{o}$ calculated with $|\psi_{0}\rangle$) is $\sqrt{\varepsilon^{2}-V^{2}}$ (odd-order subspace), therefore, $V=\varepsilon$ is the breaking point of the RPA. It is seen analytically from the table that this problem does not occur in the nonlinear higher RPA.

II.3 Numerical result

For the initial calculation of the matrix elements in the nonlinear higher RPA equations we used an ansatz for the nonlinear higher RPA ground state as follows Šimkovic et al. (2000):

[TABLE]

where ${\cal N}$ is the normalization factor, with a small arbitrary value of $Y/X$ for generating the initial $\beta_{2i}$’s by expanding the exponential operator function [($J_{+})^{i}|\psi_{0}\rangle=0$ for $i>2n$] . As $Y/X$ is small, $\beta_{2}$, $\beta_{4},\cdots$ are very small. This initial guess was applied to the calculation for a small $V$, and for a slightly larger $V$ the solution for the slightly smaller $V$ was used as the initial guess. The solutions with increasing $V$ were obtained by repeating this manner.

The initial $\beta_{2i}$’s are used for calculating the matrix elements entering Eq. (31) through Eqs. (49)$-$(57), and $X^{k}_{2i-1}$’s and $Y^{k}_{2i-1}$’s are obtained for all $k$ by solving Eq. (31). In this procedure we use a technique to diagonalize $U^{o-1}(A^{o}+B^{o})U^{o-1}(A^{o}-B^{o})$ and obtain eigenvalue of $E_{k0}^{2}$ Ring and Schuck (1980). The orthonormalization condition is

[TABLE]

and the one for the even-order subspace can be written in the same way. Then, these $X^{k}_{2i-1}$’s and $Y^{k}_{2i-1}$’s are input to Eq. (67), and $\beta_{2i}/\beta_{0}$’s are obtained. Equation (67) with $k=1$ is used (here is the arbitrarity of the choice of $k$, as mentioned above, see also below). The component $\beta_{0}$ is determined by the normalization of $|\Psi_{0}\rangle$;

[TABLE]

The $\beta_{2i}$’s obtained from Eqs. (67) and (98) are input to Eq. (31) through Eqs. (49)$-$(57), and this procedure is iterated until the convergence is obtained. The converged $\beta_{2i}$’s are input to the eigenequation of the even-order subspace corresponding to Eq. (31), and $X^{k}_{2i}$’s and $Y^{k}_{2i}$’s are obtained for all $k$; the iteration is not necessary at this stage. The $\beta_{2i}$’s, $X^{k}_{2i}$’s, and $Y^{k}_{2i}$’s are input to Eq. (70) (again there is an arbitrarity of the choice of the equation as mentioned above), and $c_{k}$ is determined for all $k$.

In the numerical calculation, $\varepsilon$ is set equal to 1, and $N=2n=8$ is used. The spectrum is shown in Fig. 1 as a function of $-V$. We assign the excited-state numbers as $E_{10}^{e}<E_{20}^{e}<E_{30}^{e}<E_{40}^{e}$ and analogously for the odd-order subspace. The RPA breaking occurs at $V=-1/(N-1)=-0.143$, and we calculated up to about twice this strength. The accuracy check with $V=-0.15$ is shown by Fig 2, which shows that our calculation reproduces the exact result of the shell model with no truncation of the wavefunction space. Other excitation energies have the accuracies in the same range. We also confirmed that the components of the ground and excited states are equal to those of the exact calculation. $V=-0.15$ is always used in the analysis of this section.

The determination of $|\Psi_{0}\rangle$ is most sensitive to the lowest excited states through the vacuum condition. This is shown by Fig. 3, which illustrates the convergence of the self-consistent calculation to the exact result obtained using the vacuum conditions with $k=1$ and 2. The lowest excited state $k=1$ was used for the vacuum condition of the calculation of Fig. 2 because of this sensitivity. Figure 4 shows that $c_{k}$’s are obtained independently of the choice of the equation. This check is satisfactory for all $l$ specifying the equation.

II.4 Truncation approximation

For realistic calculations one cannot avoid approximation in any many-body approach. We compare the quality of the approximation under the truncation of the matrices between the nonlinear higher RPA and the shell model. As mentioned above, we treat $A^{o}$, $B^{o}$, and $U^{o}$ and those of the even-order subspace for solving the equations. Thus, the dimension of these matrices and that of the $H$ matrix of the shell model are referred to in the comparison. We performed calculations with $N=8$ and $V=-0.15$ (Fig. 5a). This figure illustrates $E_{10}^{o}$ of four methods. Method A is the exact shell model (the dimension of the matrix is $d=4$). Method B is the nonlinear higher RPA with $d=2$ ($n$ is still 4). Methods C and D are the shell model with different truncations: $d=2$ (C) and 3 (D). The results of methods B and C show that the nonlinear higher RPA is much better than the shell model with the same $d$. The shell model needs $d=3$ for obtaining the quality of our method with $d=2$. This is understood from the feature of the nonlinear higher RPA. When the highest order of $J_{+}$ in $Q_{k}^{o\dagger}$ is $2i-1$, that in the wavefunction is $4i-1$. The corresponding order of the shell model is $2i-1$. Thus, this advantage of the nonlinear higher RPA is expected well in the realistic calculations. Panel b of Fig. 5 shows the $V$ dependence of $E_{10}^{o}$. The truncation approximation B and D are very good up to $V\simeq-0.2$.

We made the same comparison for $N=20$ ($n=10$), $V$ = $-$0.072, and the result is shown in panel a of Fig. 6. Method A is the exact calculation ($d$ = 10), and method B is the nonlinear higher RPA with $d$ = 2 ($n$ is still 10). Methods C, D, and E are the shell model with $d$ = 2 (C), 3 (D), and 4 (E). The same tendency as Fig. 5 is seen. But it is clear in this example with larger $N$ that the nonlinear higher RPA calculation with half the matrix size of the shell model calculation has the comparable approximation quality with the shell model calculation as is expected from the above argument. The method B ($d=2$) is equivalent to the self-consistent (extended) second RPA Schuck and Tohyama (2016) (in our calculation the truncation is the only approximation). The $V$ dependence of $E_{10}^{o}$ is shown by panel b of Fig. 6. The truncation approximation B and E are good up to $V\simeq-0.06$.

III Representation with symmetry-breaking basis

In the previous section, the vacuum condition for the odd-order subspace determines the $\beta$’s, and that for the even-order subspace is used for determining $c_{k}$. In this section, we show that the ground state can also be determined without symmetry. In fact, the formulation with no symmetry has an opportunity to use for the Lipkin model because, if $|V|$ is large, the Hartree-Fock (HF) ground state breaks the symmetry.

III.1 Hartree-Fock basis

We show the result of the HF approximation for the Lipkin model without derivation. The HF equation is given by

[TABLE]

[TABLE]

where $|\tilde{\psi}_{0}\rangle$ denotes the HF ground state, and it is identical to $|\psi_{0}\rangle$ if $|V|$ is small (see below). The eigenvectors and eigenvalues are as follows:

[TABLE]

[TABLE]

The annihilation operators $\alpha_{1m}$ and $\alpha_{0m}$ of the particles of the HF eigenstate are derived by

[TABLE]

and these operators satisfy

[TABLE]

The equation for $\langle a_{1m}^{\dagger}a_{0m}\rangle$ $\equiv$ $\langle\tilde{\psi}_{0}|a_{1m}^{\dagger}a_{0m}|\tilde{\psi}_{0}\rangle$ can be derived:

[TABLE]

and we eventually obtain

[TABLE]

The HF solutions are determined by inserting Eq. (119) or (120) to Eqs. (106) and (109). The HF ground state energy $E_{gs}^{HF}$ $=$ $\langle\tilde{\psi}_{0}|H|\tilde{\psi}_{0}\rangle$ is found to be

[TABLE]

and its behavior is drawn in Fig. 7. For $|V|$ $>$ $\varepsilon/(N-1)$, the HF solution breaks the parity symmetry of the order with respect to $J_{+}$ ($\langle J_{+}\rangle\neq 0$). As seen from Eq. (118), the two HF solutions belong to the different branches. That is, one solution cannot be obtained from another one by changing the parameters.

III.2 Nonlinear higher RPA with symmetry-breaking basis

In this section the notation of $|\tilde{\psi}_{0}\rangle$ is used for the symmetry-breaking HF ground state [$|V|$ $>$ $\varepsilon/(N-1)$]. We introduce the operators using the symmetry-breaking basis;

[TABLE]

$\tilde{J}_{z}$, $\tilde{J}_{+}$, and $\tilde{J}_{-}$ satisfy the same commutation relations as those for $J_{z}$, $J_{+}$, and $J_{-}$. In this representation, the Hamiltonian is expressed

[TABLE]

[TABLE]

As seen from the linear terms of $\tilde{J}_{+}$ and $\tilde{J}_{-}$ in Eq. (125), the odd-order subspace $\{\tilde{J}_{+}^{2i+1}|\tilde{\psi}_{0}\rangle\}$ and even-order one $\{\tilde{J}_{+}^{2i}|\tilde{\psi}_{0}\rangle\}$ are not decoupled. The eigenequation of the nonlinear higher RPA in this representation can be derived analogously to Eq. (31). Its extension to the symmetry-breaking formulation is straightforward, and we omit the explicit equations.

The expression of the vacuum condition (18) without using the symmetry is obtained

[TABLE]

[TABLE]

[TABLE]

where $\tilde{\beta}$’s are the components of the ground state

[TABLE]

and $\tilde{X}_{i}^{k}$’s and $\tilde{Y}_{i}^{k}$’s are the amplitudes obtained from the eigenequation in the symmetry-breaking representation (omitted). We applied this symmetry-breaking formulation for $N=8$, $V=-0.5$ and confirmed that the exact solutions are obtained. For calculating $c_{k}$, $\tilde{\beta}_{i}$’s of the previous cycle in the iteration can be used. In fact, Eqs. (71) and (72) with the symmetry-conserving basis can be rewritten to an expression similar to Eq. (147), that is,

[TABLE]

[TABLE]

We investigated the truncation approximation with respect to $\tilde{J}_{+}$. The correct solutions of the truncated nonlinear higher RPA were not obtained in the region of $V$ of the symmetry breaking in the HF approximation. The iteration process converged to an unphysical solution, e.g., the first excitation energy $E_{10}$ is 1.165 for $V=-0.2$ and the truncation order of 6 ($N=8$) when the shell-model truncated at the same order gives 0.703 (the exact value is 0.527). Unphysical solutions are possible because of the nonlinearity of the eigenequation. The reason for this difficulty can be discussed by analyzing the exact wavefunctions. Table 2 shows the squared norm of the lower-order components and that of the higher-order components for $V=-0.2$. The components in the symmetry-breaking basis are more distributed to the higher order than those in the symmetry-conserving basis. Therefore, the neglected components in the symmetry-breaking basis are more important than those in the symmetry-conserving basis. In addition, when the amplitudes have a broad distribution, it is difficult to have an input wavefunction close to the solution. The reason for that unexpected amplitude distribution can be inferred from the fact that the exact solution does not have the phase transition; in this model, the symmetry breaking is an artifact of approximation. Thus, one of the reasons for this problem is the property of the Lipkin model.

IV Summary

We have shown that the nonlinear higher RPA reproduces the exact solutions of the Lipkin model. The reconstruction has been demonstrated analytically for $N=2$ and numerically for $N=8$. Our study is the first one which shows the reproduction of the exact solutions of the Lipkin model for arbitrary $N$ within an extension of the RPA method. The proper construction of the phonon operator is crucial. We examined the results carefully and conclude that every property of the solutions is consistent. The only approximation in the realistic applications is the truncation of the wavefunction space, thus, there is no possibility of the lack of physical effects due to the mathematical properties of the nonlinear higher RPA. Considering that the neck point of the shell model is its huge matrix dimension, the advantage of the nonlinear higher RPA under the truncation is encouraging.

We have also shown the formulation with the HF basis breaking the symmetry and reproduced the exact solution. However, the truncated solutions with this basis could not be obtained. It is an open question whether this problem occurs in the realistic systems having phase trnasitions.

For the feasibility of the realistic calculation, there are calculations of nuclei by the second RPA Catara et al. (1994); Papakonstantinou and Roth (2010). The feasibility of the iteration for solving the nonlinear second RPA is a matter of computational resource. Considering that the significant progress of the computers continues, the realistic application is a near-future task. It is also possible to consider a simplification Smetana et al. (2015) by introducing an approximate ansatz for the ground state of the nonlinear higher RPA.

Acknowledgments

This work is supported by the VEGA Grant Agency of the Slovak Republic under Contract No. 1/0922/16, the Ministry of Education, Youth and Sports of the Czech Republic under Contract No. LM2011027, RFBR Grant No. 16-02-01104, Underground laboratory LSM - Czech participation to European-level research infrastructure cz.02.1.01$/$0.0$/$0.0$/$16_013$/$0001733, and Grant No. HLP-2015-18 of the Heisenberg-Landau Program. M.I.K. acknowledges the support of Alexander von Humboldt Foundation. This work is also supported by Bogoliubov Laboratory of Theoretical Physics, Joint Institute for Nuclear Research through the stay of one of the authors there. The computer COMA of Center for Computational Sciences, University of Tsukuba was used through the Interdisciplinary Public Program of fiscal year 2016 (TKBNDFT) for the numerical calculation.