Excited-state quantum phase transitions in systems with two degrees of freedom: III. Interacting boson systems

TL;DR

This paper investigates excited-state quantum phase transitions in interacting boson systems, revealing complex classical stationary points and finite-size effects, thus advancing understanding of many-body quantum dynamics.

Contribution

It extends previous work on ESQPTs by analyzing a realistic interacting boson model with fixed angular momentum, highlighting the role of classical stationary points and phase space boundedness.

Findings

ESQPTs linked to classical stationary points

Finite-size effects due to partial separability

Complexity from kinetic and potential energy decomposition

Abstract

The series of articles [Ann. Phys. 345, 73 (2014) and 356, 57 (2015)] devoted to excited-state quantum phase transitions (ESQPTs) in systems with $f=2$ degrees of freedom is continued by studying the interacting boson model of nuclear collective dynamics as an example of a truly many-body system. The intrinsic Hamiltonian formalism with angular momentum fixed to $L=0$ is used to produce a generic first-order ground-state quantum phase transition with an adjustable energy barrier between the competing equilibrium configurations. The associated ESQPTs are shown to result from various classical stationary points of the model Hamiltonian, whose analysis is more complex than in previous cases because of (i) a non-trivial decomposition to kinetic and potential energy terms and (ii) the boundedness of the associated classical phase space. Finite-size effects resulting from a partial separability of both degrees of freedom are analyzed. The features studied here are inherent in a great majority of interacting boson systems.

Figures

Tables

Table 1. Table 1: Classical limit of the U( n + 1 𝑛 1 n\!+\!1 ) generators.

		symmetric generators	antisymmetric generators
$k$	$l$	${\hat{b}}_k^{†}{\hat{b}}_l+{\hat{b}}_l^{†}{\hat{b}}_k$	$i\left({\hat{b}}_k^{†}{\hat{b}}_l-{\hat{b}}_l^{†}{\hat{b}}_k\right)$
$>0$	$>0$	$N\left(q_k{q}_l+p_k{p}_l\right)$	$N\left(p_k{q}_l-q_k{p}_l\right)$
$>0$	$=0$	$N{q}_k\sqrt{2-\sum _{m=1}^n\left(q_m^2+p_m^2\right)}$	$N{p}_k\sqrt{2-\sum _{m=1}^n\left(q_m^2+p_m^2\right)}$
$=0$	$>0$	$N{q}_l\sqrt{2-\sum _{m=1}^n\left(q_m^2+p_m^2\right)}$	$-N{p}_l\sqrt{2-\sum _{m=1}^n\left(q_m^2+p_m^2\right)}$
$=0$	$=0$	$N\left(2-\sum _{m=1}^n\left(q_m^2+p_m^2\right)\right)$	$0$

Equations

$${\hat{H}}^{\lambda }=\sum _{k,k^{\prime }=0}^n{\epsilon }_{k{k}^{\prime }}^{\lambda }{\hat{b}}_k^{†}{\hat{b}}_{k^{\prime }}+\frac{1}{N-1}\sum _{\begin{array}{c}k,k^{\prime }\\ l^{\prime },l\end{array}\}=0}^n{v}_{k{k}^{\prime }l{l}^{\prime }}^{\lambda }{\hat{b}}_k^{†}{\hat{b}}_{k^{\prime }}^{†}{\hat{b}}_{l^{\prime }}{\hat{b}}_l,$$

$${\hat{b}}_k^{†}=\sqrt{\frac{N}{2}}({\hat{q}}_k-i{\hat{p}}_k),{\hat{b}}_k=\sqrt{\frac{N}{2}}({\hat{q}}_k+i{\hat{p}}_k),$$

$$\bigl{[}\hat{q}_{k},\hat{p}_{l}\bigr{]}=\frac{i}{N}\delta_{kl}$$

$$\lim_{N\to\infty}\frac{\hat{H}^{\lambda}}{N}\ \mapsto\ {\cal H}^{\lambda}\bigl{(}\{\hat{q}_{k},\hat{p}_{k}\}\bigr{)}\,,$$

$${\hat{a}}_0^{†}=\sqrt{\frac{\hat{N}}{{\hat{b}}_0^{†}{\hat{b}}_0}}{\hat{b}}_0^{†}{\hat{a}}_k^{†}={\hat{b}}_0\sqrt{\frac{1}{{\hat{b}}_0^{†}{\hat{b}}_0}}{\hat{b}}_k^{†}(\mathrm{for}k>0),$$

$${\hat{b}}_0^{†}=\sqrt{\frac{{\hat{a}}_0^{†}{\hat{a}}_0-\sum _{l=1}^n{\hat{a}}_l^{†}{\hat{a}}_l}{{\hat{a}}_0^{†}{\hat{a}}_0}}{\hat{a}}_0^{†},{\hat{b}}_k^{†}=\sqrt{\frac{1}{{\hat{a}}_0^{†}{\hat{a}}_0}}{\hat{a}}_0^{†}{\hat{a}}_k^{†}(k>0)$$

$$R^{2}\equiv\sum_{k=1}^{n}\bigl{(}q_{k}^{2}+p_{k}^{2}\bigr{)}\leq 2\equiv R_{0}^{2}\,,$$

$${\mathcal{H}}^{\lambda }(\alpha )=\frac{⟨\alpha |{\hat{H}}^{\lambda }|\alpha ⟩}{⟨\alpha |\hat{N}|\alpha ⟩}.$$

$$|\alpha ⟩=e^{-\frac{\Vert \alpha {\Vert }^2}{2}}{e}^{\Vert \alpha \Vert {\hat{B}}^{†}(\tilde{\alpha })}|0⟩,{\hat{B}}^{†}(\tilde{\alpha })=\sum _{k=0}^n{\tilde{\alpha }}_k{\hat{b}}_k^{†}.$$

$$|\tilde{\alpha },N⟩=\frac{1}{\sqrt{N!}}{\hat{B}}^{†}(\tilde{\alpha }{)}^N|0⟩,$$

$${\mathcal{H}}^{\lambda }(\tilde{\alpha },N)=\frac{1}{N}⟨\tilde{\alpha },N|{\hat{H}}^{\lambda }|\tilde{\alpha },N⟩.$$

$${\overline{\rho }}^{\lambda }(E)=\sum _i\overline{\delta }(E-E_i^{\lambda }),$$

$${\overline{j}}^{\lambda }(E)={\overline{\rho }}^{\lambda }(E){\overline{\varphi }}^{\lambda }(E)=\sum _i\frac{\partial E_i^{\lambda }}{\partial \lambda }\overline{\delta }(E-E_i^{\lambda }),$$

$$\frac{\partial E_i^{\lambda }}{\partial \lambda }=⟨{\psi }_i^{\lambda }|\frac{\partial {\hat{H}}^{\lambda }}{\partial \lambda }|{\psi }_i^{\lambda }⟩,$$

$${\overline{\rho }}^{\lambda }(E)\propto \frac{\partial }{\partial E}\underset{\mathbf{\Omega }}{\iint }{d}^f\mathbf{q}{d}^f\mathbf{p}\Theta \left[E-{\mathcal{H}}^{\lambda }(\mathbf{q},\mathbf{p})\right],$$

$$\frac{\partial\overline{\rho}^{\lambda}}{\partial E}\propto\left\{\begin{array}[]{ll}(-1)^{\frac{r}{2}}\ \Theta(E\!-\!E_{\rm c}^{\lambda})&{\rm for}\ r\!=\!0,2,4,\\ (-1)^{\frac{r+1}{2}}\ \ln|E\!-\!E_{\rm c}^{\lambda}|&{\rm for}\ r\!=\!1,3,\end{array}\right.$$

$$\frac{\partial\overline{j}^{\lambda}}{\partial E}\!\equiv\!\frac{\partial\,\overline{\partial_{\lambda}H}^{\lambda}}{\partial E}\propto\left\{\begin{array}[]{ll}\Theta(E\!-\!E_{\rm c}^{\lambda})&{\rm for}\ r\!=\!0,2,4,\\ \ln|E\!-\!E_{\rm c}^{\lambda}|&{\rm for}\ r\!=\!1,3.\end{array}\right.$$

$${\mathcal{H}}^{\lambda }(\mathbf{q},\mathbf{p})=\underset{{\mathcal{K}}_0^{\lambda }(\mathbf{q},\mathbf{p})}{\underbrace{{\mathcal{H}}^{\lambda }(\mathbf{q},\mathbf{p})-{\mathcal{H}}^{\lambda }(\mathbf{q},0)}}+\underset{{\mathcal{V}}_0^{\lambda }(\mathbf{q})}{\underbrace{{\mathcal{H}}^{\lambda }(\mathbf{q},0)}}.$$

$${\frac{\partial }{\partial q_k}{\mathcal{V}}_0^{\lambda }(\mathbf{q})|}_{\mathbf{q}={\mathbf{q}}_0^{\lambda }}=0.$$

$${\frac{\partial {\mathcal{K}}_0^{\lambda }(\mathbf{q},\mathbf{p})}{\partial q_k}|}_{\begin{array}{c}\mathbf{q}={\mathbf{q}}_{0i}^{\lambda }\\ \mathbf{p}=\pm {\mathbf{p}}_{0i}^{\lambda }({\mathbf{q}}_{0i}^{\lambda })\end{array}}+{\frac{\partial {\mathcal{V}}_0^{\lambda }(\mathbf{q})}{\partial q_k}|}_{\mathbf{q}={\mathbf{q}}_{0i}^{\lambda }}=0$$

$$E^{\lambda }(R,\mathbf{\varphi })=a^{\lambda }(R_0^2-R^2)+b^{\lambda }(\mathbf{\varphi })R(R_0^2-R^2{)}^{\frac{1}{2}}+c^{\lambda }(\mathbf{\varphi })R^2$$

$$+d^{\lambda }(\mathbf{\varphi })R^4+e^{\lambda }(\mathbf{\varphi })R^3(R_0^2-R^2{)}^{\frac{1}{2}}+f^{\lambda }(\mathbf{\varphi })R^2(R_0^2-R^2)$$

$$+g^{\lambda }(\mathbf{\varphi })R(R_0^2-R^2{)}^{\frac{3}{2}}+h^{\lambda }(R_0^2-R^2{)}^2.$$

$$E_{\partial \Omega }^{\lambda }(\mathbf{\varphi })=c^{\lambda }(\mathbf{\varphi })R_0^2+d^{\lambda }(\mathbf{\varphi })R_0^4,$$

$${\overline{\rho }}^{\lambda }(E)={\varrho }^{\lambda }(E)\left[1-\sigma \Theta \left[\pm (E-E_{\mathrm{c}}^{\lambda })\right]{\left|E-E_{\mathrm{c}}^{\lambda }\right|}^I\right],$$

$$I=\frac{2f-1}{K}+\frac{1}{M}-1,$$

$${\mathcal{H}}_{\mathrm{ex}}^{\lambda }(\tilde{\alpha },N-N_{\perp },N_{\perp })=$$

$$\frac{⟨{\tilde{\alpha }}_{\perp },N_{\perp }|\hat{B}(\tilde{\alpha }{)}^{N-N_{\perp }}{\hat{H}}^{\lambda }{\hat{B}}^{†}(\tilde{\alpha }{)}^{N-N_{\perp }}|{\tilde{\alpha }}_{\perp },N_{\perp }⟩}{(N-N_{\perp })!N},$$

$${\hat{H}}^{(\eta ,\chi )}=(1-\eta ){\hat{n}}_d-\frac{\eta }{4N}({\hat{Q}}^{\chi }\cdot {\hat{Q}}^{\chi })$$

$${\hat{H}}_1({\beta }_0^{\prime },\zeta )$$

$${\hat{H}}_2({\beta }_0^{\prime },\xi )$$

$${\hat{S}}^{†}({\beta }_0^{\prime })$$

$${\hat{D}}_{\mu }^{†}({\beta }_0^{\prime },\zeta )$$

Full text

I Introduction

Excited-state quantum phase transitions (ESQPTs) ref:Cejn06 ; ref:Capr08 ; ref:Cejn08 ; ref:Cejn09 ; ref:E1 ; ref:E2 ; ref:Stra16 represent an extension of the ground-state quantum phase transitions (QPTs) ref:Carr11 to the eigenstates with higher energy. ESQPTs were introduced as robust singularities in the spectra of energy levels in the nuclear interacting boson model (IBM) ref:Cejn06 ; ref:Hein06 ; ref:Mace06 , and have since been studied theoretically in numerous other many-body models like the Lipkin model, see, e.g., Refs. ref:Ribe07 ; ref:Opat18 , the molecular vibron model ref:Lare11 ; ref:Lare13 , the fermionic and bosonic pairing models ref:Skra11 ; ref:Rela16 , the extended Dicke model of superradiance ref:Bran13 ; ref:Basta14 ; ref:Kloc17a ; ref:Kloc17b , the Bose-Hubbard model of atom-molecule condensates ref:Rela14 ; ref:Bych18 , and in algebraic models of two-dimensional crystals ref:Diet13 ; ref:Iach15 ; ref:Diet17 . The ESQPTs were experimentally observed in molecules ref:Lare11 ; ref:Lare13 and in some artificial quantum systems like photonic crystals ref:Diet13 . Dynamical consequences of ESQPTs in the response to various kinds of driving ref:Ferna11 ; ref:Basti14 ; ref:Pueb15 ; ref:Kopy15 ; ref:Sant16 ; ref:Kopy17 ; ref:Kloc18 or in the entanglement and decoherence properties ref:Rela08 ; ref:Kloc17a are of potential relevance for quantum simulation and computation applications. Also the thermodynamic anomalies related to ESQPTs ref:Basta16 ; ref:Fern16 ; ref:Cejn17 are of fundamental interest.

For a general bound quantum system with quantum Hamiltonian $\hat{H}^{\lambda}$ depending on a control parameter $\lambda$, the ESQPT represents a singularity in the average density and slope of discrete energy levels $E_{i}^{\lambda}$ and in the overall structure of the corresponding eigenstates $|\psi_{i}^{\lambda}\rangle$ (in this paper, the dependence on the control parameter is always expressed by the superscript $\lambda$, which should not be confused with the power symbol). The singularity is caused by a quasi-stationary point of the underlying classical dynamics. A general classification of such singularities for non-degenerate stationary points in systems with an arbitrary number of degrees of freedom $f$ was given in Ref. ref:Stra16 . A systematic analysis of ESQPTs in systems with $f=2$ was initiated in Refs. ref:E1 and ref:E2 . In these articles, further referred as part I and part II, respectively, the focus was set to systems whose classical limit is expressed via coordinates $\boldsymbol{q}$ and momenta $\boldsymbol{p}$ in an unbounded phase space in which the Hamiltonian function ${\cal H}^{\lambda}(\boldsymbol{q},\boldsymbol{p})$ can be cast as a sum of a simple kinetic energy term ${\cal K}(\boldsymbol{p})\!\propto\!\boldsymbol{p}^{2}$ and an arbitrary potential ${\cal V}^{\lambda}(\boldsymbol{q})$. However, in the classical limit of quantum many-body systems the phase space is usually bounded and the Hamiltonian decomposition is more complicated. In such systems, the ESQPT analysis shows some peculiarities.

The present article (part III) continues the series of Refs. ref:E1 ; ref:E2 by considering a suitable example of a many-body system with $f=2$, namely the IBM ref:IBMbook with angular momentum set to $L=0$. The transition between the U(5) and O(6) dynamical symmetries of this model, which is connected with the second-order ground-state QPT, was at the beginning of the ESQPT development ref:Cejn06 ; ref:Hein06 ; ref:Mace06 . More recently, ESQPTs along the transition between the U(5) and SU(3) dynamical symmetries with the first-order ground-state QPT were studied for the $L=0$ subset of the IBM spectrum ref:Zhan16 , as well as with respect to states with $L>0$ ref:Zhan17 . Our present article also deals with the ESQPT structures accompanying the first-order QPT between spherical and axially deformed ground-state configurations. The focus on a first-order QPT Hamiltonian is in agreement with the general strategy of parts I and II. However, the setup of the model is different than in Refs. ref:Zhan16 ; ref:Zhan17 , that employed a simple “consistent-$Q$” IBM Hamiltonian. Instead we apply the formalism introduced in Refs. ref:Kirs85 ; ref:Levi87 in the context of resolution of the IBM dynamics into intrinsic and collective components. This approach turned out to be very useful in the study of first-order QPTs ref:Levi06 ; ref:Mace11 ; ref:Levi12 ; ref:Mace14 and we use here its advantage of providing an adjustable (preferably high) energy barrier between the competing ground-state configurations. This facilitates our study of typical ESQPT effects associated with generic first-order QPTs in quantum many-body systems.

The article is organized as follows: In Sec. II, we consider a general interacting boson system, construct its classical limit, and describe new general aspects of the ESQPT analysis that follow from the compactness of the relevant phase space and from non-separability of coordinate and momentum Hamiltonian terms. In Sec. III, as a case study demonstrating general effects introduced in Sec. II, we present a numerical analysis of ESQPTs in a family of IBM Hamiltonians describing the spherical–deformed nuclear shape-phase transition. In Sec. IV we give a brief summary and outlook.

II General interacting boson systems

II.1 Hamiltonian and classical limit

This article deals with systems of interacting bosons. It may be extended to systems of interacting fermions by identifying a particular fermionic collective algebra and realizing it in terms of bosonic operators. This is usually possible, but not discussed here. We further assume that the total number of particles (bosons) $N$ is conserved, so that it can naturally play the role of the size parameter $\aleph$ used in parts I and II ref:E1 ; ref:E2 , hence we set $\aleph\equiv N$.

In the following, we consider bosonic Hamiltonians with one- and two-body terms, namely

[TABLE]

where $\hat{b}_{k}^{{\dagger}}$ and $\hat{b}_{k}$ are creation and annihilation operators of bosons in $n\!+\!1$ single-particle states enumerated by $k=0,1,2,\dots,n$ (with $n\geq 1$). The $(n\!+\!1)^{2}$ bilinear products of boson creation and annihilation operators, or more precisely symmetric and antisymmetric Hermitian operators of the form $\hat{b}_{k}^{{\dagger}}\hat{b}_{l}+\hat{b}_{l}^{{\dagger}}\hat{b}_{k}$ and i\bigl{(}\hat{b}_{k}^{{\dagger}}\hat{b}_{l}-\hat{b}_{l}^{{\dagger}}\hat{b}_{k}\bigr{)}, represent generators of the system’s dynamical algebra U($n\!+\!1$). Coefficients $\varepsilon^{\lambda}_{kk^{\prime}}$ and $v^{\lambda}_{kk^{\prime}ll^{\prime}}$ in Eq. (1) are assumed to depend on an adjustable control parameter $\lambda$ and must always satisfy the constraints following from the Hermiticity requirement. The $N\!-\!1$ denominator ensures that the proportion between the one- and two-body terms is preserved for varying total number of bosons.

The Hamiltonian (1) can be cast in the coordinate–momentum form with the aid of transformation

[TABLE]

where $\{\hat{q}_{k},\hat{p}_{k}\}_{k=0}^{n}$ are operators of canonically conjugate coordinates and momenta. The commutator

[TABLE]

indicates that $N^{-1}$ represents an effective Planck constant of the system and $N\to\infty$ corresponds to its classical limit. The classical Hamiltonian ${\cal H}^{\lambda}$ associated with $\hat{H}^{\lambda}$ in Eq. (1) is obtained via the mapping

[TABLE]

where the substitution (2) is made in $\hat{H}^{\lambda}$ with the confidence that for ${N\to\infty}$ all the $\hat{q}_{k}$ and $\hat{p}_{l}$ terms mutually commute and become just ordinary numbers $q_{k}$ and $p_{l}$. These variables constitute a classical phase space, which is finite due to the conservation of ${\hat{N}=\sum_{k=0}^{n}\hat{b}_{k}^{{\dagger}}\hat{b}_{k}}$. Indeed, the relation expressing this constraint, {\sum_{k=0}^{n}\bigl{(}q_{k}^{2}+p_{k}^{2}\bigr{)}=2}, implies that this kind of phase space coincides with the surface of a $2(n\!+\!1)$-dimensional sphere with radius $\sqrt{2}$. But this is not yet the space applied in most of calculations.

One path to the proper classical limit of Hamiltonian (1) leads via a generalization of the Holstein-Primakoff transformation ref:Blai78 . The constraint on conserved observable $\hat{N}$ is used for a complete elimination of one degree of freedom, e.g. that associated with boson $b_{0}$. In particular, bosons $b_{k}$ are transformed to new bosons $a_{k}$,

[TABLE]

with the inverse transformation given by:

[TABLE]

(note that the operator denominators in these expressions do not produce singularities). Since the $a_{k}$ operators satisfy bosonic commutation relations, the transformation is canonical. However, the total number of the new bosons does not match with $\hat{N}$ and is not conserved by Hamiltonian (1). Instead we have ${\hat{a}_{0}^{{\dagger}}\hat{a}_{0}=\hat{N}}$ and ${\sum_{l=1}^{n}\hat{a}_{l}^{{\dagger}}\hat{a}_{l}=\hat{N}-\hat{b}_{0}^{{\dagger}}\hat{b}_{0}}$. This means that the conserved number $\hat{N}$ is “stored” solely in boson $a_{0}$ while the active degrees of freedom are represented by the remaining bosons $a_{k}$ with ${k=1,...,n}$. Hence we end up with a system with ${f=n}$ degrees of freedom.

Applying the coordinate–momentum decomposition (2) to all $a_{k}$ bosons from Eq. (5), we obtain the classical-limit expressions of generators of the dynamical group ${{\rm U}(n+1)}$ in terms of coordinates and momenta ${\{q_{k},p_{k}\}_{k=1}^{n}\equiv\{\boldsymbol{q},\boldsymbol{p}\}}$ associated with bosons $\{a_{k}\}_{k=1}^{k}$. These expressions are given in Table 1. We see that the classical representation is the same as if we applied the decomposition (2) to the original bosons $b_{k}$ with a substitution ${\hat{b}_{0}^{{\dagger}},\hat{b}_{0}\mapsto\sqrt{2-\sum_{l=1}^{f}(q_{l}^{2}+p_{l}^{2})}}$. Note that the new coordinates and momenta must satisfy the condition:

[TABLE]

which means that the phase space with ${f=n}$ is the interior of a $2f$-dimensional sphere with radius ${R_{0}=\sqrt{2}}$. Hereafter this domain is denoted by symbol ${\bf\Omega}$.

Decomposing an arbitrary quantum Hamiltonian (1) into the Hermitian generators of ${{\rm U}(n+1)}$, we can use the expressions in Table 1 to write down its classical counterpart. A common alternative to this construction is based on the use of various kinds of coherent states ref:Zhan90 ; ref:Iach14 , which finally does the same job. In this case, the classical Hamiltonian ${\cal H}^{\lambda}$ is obtained from the energy-per-boson expectation value

[TABLE]

in a coherent state $|\alpha\rangle$, with ${\alpha\equiv\{\alpha_{k}\}_{k=0}^{n}}$ denoting a set of complex parameters composed of real coordinate and momentum variables $q_{k}$ and $p_{k}$. In particular, one can use the Glauber coherent states given by:

[TABLE]

Here $|0\rangle$ represents the boson vacuum and $\hat{B}^{{\dagger}}(\tilde{\alpha})|0\rangle$ is a normalized single-boson state, which is in the basis $\{\hat{b}_{k}|0\rangle\}_{k=0}^{n}$ expressed by expansion coefficients ${\{\tilde{\alpha}_{k}\equiv\alpha_{k}/\lVert\alpha\rVert\}}$ with the normalization denominator ${\lVert\alpha\rVert=\sqrt{\sum_{k=0}^{n}|\alpha_{k}|^{2}}}$. For the state (9), the average number of bosons is ${\langle N\rangle_{\alpha}\!\equiv\!\langle\alpha|\hat{N}|\alpha\rangle=\lVert\alpha\rVert^{2}}$. Substituting ${\alpha_{k}=\sqrt{\frac{\langle N\rangle_{\alpha}}{2}}(q_{k}+ip_{k})}$ for parameters with ${k=1,...,n}$ and setting ${\alpha_{0}=\sqrt{\langle N\rangle_{\alpha}-\sum_{l=1}^{n}|\alpha_{k}|^{2}}}$, one arrives at the same expression as obtained by the above procedure with $N$ replaced by $\langle N\rangle_{\alpha}$.

The projective coherent states,

[TABLE]

i.e., condensates of a fixed number $N$ of bosons in the state $\hat{B}^{{\dagger}}(\tilde{\alpha})|0\rangle$, yield equivalent results with

[TABLE]

Note that the same can be achieved also with the other standard forms of coherent states (so-called algebraic and group coherent states) associated with the factorization ${{\rm U}(n+1)/[{\rm U}(n)\otimes{\rm U}(1)]}$ of the dynamical algebra ref:Iach14 .

The classical limit derived above holds for all boson systems that conserve $\hat{N}$. Additional conservation laws may further reduce the number of effective degrees of freedom. In this article series we investigate systems with ${f=2}$. The bosonic model discussed in Sec. III has ${n=5}$, which corresponds to the full number of the model’s degrees of freedom. However, the invariance of the Hamiltonian under rotations in the 3-dimensional space makes it possible to fix the angular momentum to ${L=0}$, which reduces the number of active degrees of freedom by three.

II.2 ESQPT analysis

II.2.1 ESQPTs due to non-degenerate stationary points

Consider a general bound quantum system with Hamiltonian $\hat{H}^{\lambda}$ and discrete energy spectrum $E_{i}^{\lambda}$ enumerated by integer index $i$ and depending on a real control parameter $\lambda$. As discussed in parts I and II ref:E1 ; ref:E2 , the ESQPT analysis relies on a decomposition $\bullet=\overline{\bullet}+\widetilde{\bullet}$ of the relevant quantities into smooth and oscillatory components. In particular, we define the smoothed level density

[TABLE]

where $\overline{\delta}$ is a smoothed $\delta$ function (e.g. a Gaussian of width exceeding the typical spacing between levels), and the smoothed level flow

[TABLE]

where $\overline{\phi}^{\lambda}(E)$ represents a “velocity field” of the spectrum. It was shown in part I ref:E1 that the smoothed level density and flow satisfy the ordinary continuity equation ensuring that levels do not disappear when $\lambda$ is varied. Substituting into Eq. (13) the Hellmann-Feynman formula,

[TABLE]

with $|\psi_{i}^{\lambda}\rangle$ standing for the eigenvector associated with level $E_{i}^{\lambda}$, we find that $\overline{j}^{\lambda}(E)\!=\!\overline{\partial_{\lambda}H}^{\lambda}(E)$, where the symbol on the right-hand side denotes the smoothed expectation value (strength function) of observable $\frac{\partial}{\partial\lambda}\hat{H}^{\lambda}$.

The above-introduced smoothed quantities play a crucial role in the definition and classification of ESQPTs. The smoothed level density coincides with the semiclassical density of states obtained by integration of the classical Hamiltonian ${\cal H}^{\lambda}(\boldsymbol{q},\boldsymbol{p})$ over the phase space, namely

[TABLE]

where $\Theta(x)$ represents the step function ($\Theta\!=\!0$ for $x\!<\!0$ and $\Theta\!=\!1$ for $x\!\geq\!0$) and ${\bf\Omega}$ is the domain of allowed $\{\boldsymbol{q},\boldsymbol{p}\}$ values. A stationary point of the classical Hamiltonian yields a singularity of $\overline{\rho}^{\lambda}(E)$ as a function of energy, and this singularity propagates via Eqs. (13) and (14) also to energy dependencies of $\overline{j}^{\lambda}(E)$ and $\overline{\partial_{\lambda}H}^{\lambda}(E)$.

The classification of ESQPTs caused by non-degenerate stationary points (i.e., those with locally quadratic dependencies of ${\cal H}^{\lambda}$ on all $\boldsymbol{q}$ and $\boldsymbol{p}$ components) was given in Ref. ref:Stra16 . The ESQPT non-analyticity associated with such points appears in the $(f\!-\!1)$th derivative of the respective quantities and depends on the index $r$ of the stationary point (the number of negative eigenvalues of the Hessian matrix). The $\lambda$-dependent energy ${\cal H}^{\lambda}=E^{\lambda}_{\rm c}$ of the stationary point forms a critical borderline (phase separatrix) in the plane $\lambda\times E$. For systems with ${f=2}$, the non-degenerate stationary points cause singularities in the first derivative of the smoothed level density,

[TABLE]

and in the first-derivative of the smoothed level slope:

[TABLE]

Note that the sign prefactors in Eq. (17) are not determined and may differ from those in Eq. (16). Degenerate stationary points (flatter than quadratic) may cause even sharper singularities, which however are not generally classified ref:Stra16 .

Our intention is to perform an ESQPT analysis of general interacting boson systems introduced in Sec. II.1. We saw (cf. Table 1) that the classical-limit Hamiltonian ${\cal H}^{\lambda}$ of such a system with one- plus two-body interactions is a quartic function of coordinates and momenta $\{\boldsymbol{q},\boldsymbol{p}\}$ defined within a compact phase-space domain ${\bf\Omega}$ from Eq. (7). It turns out that to find stationary points of this setup is more complicated than in the cases studied in parts I and II ref:E1 ; ref:E2 for two reasons: First, the analysis cannot be reduced just to the coordinate space due to complicated kinetic terms of ${\cal H}^{\lambda}$. Second, the compactness of the domain ${\bf\Omega}$ may create additional ESQPT singularities not connected to stationary points. These issues are outlined in the following two subsections.

II.2.2 Role of coordinate-dependent kinetic terms

Stationary points of a common Hamiltonian ${{\cal H}^{\lambda}(\boldsymbol{q},\boldsymbol{p})={\cal K}^{\lambda}(\boldsymbol{p})+{\cal V}^{\lambda}(\boldsymbol{q})}$, with the kinetic term ${\cal K}^{\lambda}$ being a quadratic function of momenta and the potential ${\cal V}^{\lambda}$ an arbitrary function of coordinates, are always found at ${(\boldsymbol{q},\boldsymbol{p})=(\boldsymbol{q}^{\lambda}_{0},0)}$, where $\boldsymbol{q}^{\lambda}_{0}$ are stationary points of the potential. This reduces the stationary-point analysis from the full $2f$-dimensional phase space just to the corresponding $f$-dimensional configuration space. The study of ESQPTs benefits from this reduction, see parts I and II ref:E1 ; ref:E2 . However, Hamiltonians corresponding to generic interacting boson systems contain intermixed combinations of coordinates and momenta that typically yield coordinate-dependent kinetic terms. This makes the analysis more complicated.

In general, the potential energy term ${\cal V}^{\lambda}(\boldsymbol{q})$ is extracted from the total Hamiltonian ${\cal H}^{\lambda}(\boldsymbol{q},\boldsymbol{p})$ of any form by setting ${\dot{\boldsymbol{q}}=0}$, so ${\frac{\partial}{\partial p_{k}}{\cal H}^{\lambda}=0}$ for all ${k=1,...,f}$. This condition is satisfied for some momenta ${\boldsymbol{p}=\boldsymbol{p}_{0}^{\lambda}(\boldsymbol{q})}$ and we can set ${{\cal V}^{\lambda}(\boldsymbol{q})\equiv{\cal H}^{\lambda}(\boldsymbol{q},\boldsymbol{p}^{\lambda}_{0}(\boldsymbol{q}))}$ and ${{\cal K}^{\lambda}(\boldsymbol{q},\boldsymbol{p})\equiv{\cal H}^{\lambda}(\boldsymbol{q},\boldsymbol{p})-{\cal V}^{\lambda}(\boldsymbol{q})}$. However, there may exist more than one solutions $\boldsymbol{p}^{\lambda}_{0}(\boldsymbol{q})$, which makes the above decomposition ambiguous.

The construction described in Sec. II.1 ensures that the classical limit of any bosonic Hamiltonian (1) contains only even powers of momentum $\boldsymbol{p}$. Therefore, the set of nonlinear equations following from the ${\dot{\boldsymbol{q}}=0}$ condition has solutions of two types: (i) the trivial solution ${\boldsymbol{p}^{\lambda}_{0}=0}$, which is always present and independent of $\boldsymbol{q}$, and (ii) non-trivial solutions, which come in sign conjugated pairs ${\boldsymbol{p}_{0}^{\lambda}=\pm\boldsymbol{p}_{0i}^{\lambda}(\boldsymbol{q})}$ enumerated by index $i$. The latter solutions may exist only in a limited domain of $\boldsymbol{q}$ or may be absent altogether, so their number cannot be set from general arguments. The trivial solution is usually employed in the decomposition of the Hamiltonian to the kinetic and potential terms:

[TABLE]

Specifying the solution $\boldsymbol{p}^{\lambda}_{0}(\boldsymbol{q})$ of $\dot{\boldsymbol{q}}\!=\!0$, one needs to complete the determination of stationary points by solving the ${\dot{\boldsymbol{p}}=0}$ condition, so ${\frac{\partial}{\partial q_{k}}{\cal H}^{\lambda}=0}$ for all $k=1,...,f$. This is certainly achieved by the usual procedure of finding $\boldsymbol{q}^{\lambda}_{0}$ such that

[TABLE]

Indeed, the points ${(\boldsymbol{q},\boldsymbol{p})=(\boldsymbol{q}^{\lambda}_{0},0)}$ are stationary points of the whole Hamiltonian as we obviously have ${\frac{\partial}{\partial q_{k}}{\cal K}^{\lambda}_{0}=0}$ at $\boldsymbol{p}\!=\!0$ for any $\boldsymbol{q}$. However, within the domains of existence of the non-trivial solutions $\pm\boldsymbol{p}_{0i}^{\lambda}(\boldsymbol{q})$ of ${\dot{\boldsymbol{q}}=0}$, there may be some additional solutions of the ${\dot{\boldsymbol{p}}=0}$ problem. These read as ${(\boldsymbol{q},\boldsymbol{p})=(\boldsymbol{q}_{0i}^{\lambda},\pm\boldsymbol{p}_{0i}^{\lambda}(\boldsymbol{q}_{0i}^{\lambda}))}$ and have to satisfy

[TABLE]

for all ${k=1,...,f}$. This condition follows directly from the full Hamiltonian decomposition in Eq. (18).

We can therefore conclude that Hamiltonians with coordinate-dependent kinetic terms may generate stationary points beyond the analysis of the potential energy term alone. In the study of ESQPTs in such systems, both Eqs. (19) and (20) must be taken into account.

II.2.3 Effects of phase-space boundary

The boundedness of the phase space implies the existence of an upper limit $E^{\lambda}_{\rm max}$ of the energy spectrum, which for $N\to\infty$ coincides with the maximal value ${\cal H}^{\lambda}_{\rm max}$ of the classical Hamiltonian within the domain ${\bf\Omega}$. For ${E>{\cal H}^{\lambda}_{\rm max}}$, formula (15) ensures ${\overline{\rho}^{\lambda}(E)=0}$. It is therefore legitimate to expect that the level density $\overline{\rho}^{\lambda}(E)$ exhibits the following broad behavior: it grows from zero at the minimal (ground-state) energy $E^{\lambda}_{\rm min}$, reaches a maximum somewhere in the middle of the spectrum, and then decreases back to zero above $E^{\lambda}_{\rm max}$. Note that neither the increasing, nor the decreasing part of $\overline{\rho}^{\lambda}(E)$ must be monotonic. We know (see Secs. II.2.1 and II.2.2) that if the domain ${\bf\Omega}$ contains stationary points of ${\cal H}^{\lambda}$, they generate the respective types of ESQPTs. Here we are going to show that the phase space boundary itself, i.e., the hypersurface ${\bf\partial\Omega}$, may create some extra singularities.

To analyze this effect, we introduce in the $2f$-dimensional phase space the set of hyperspherical coordinates ${R\in[0,\infty)}$ (radius) and ${\{\phi_{l}\}_{l=1}^{2f-1}\equiv\boldsymbol{\phi}}$ (angles). The energy landscape derived from the classical Hamiltonian ${\cal H}^{\lambda}(\boldsymbol{q},\boldsymbol{p})$ for ${R\in[0,R_{0}]}$ (i.e., inside ${\bf\Omega}$) then reads:

[TABLE]

Here, $a^{\lambda}$, $b^{\lambda}$, …, $h^{\lambda}$ are some coefficients (constants or functions of angles $\boldsymbol{\phi}$), which are uniquely determined from the set of constants $\varepsilon^{\lambda}_{kk^{\prime}}$ and $v^{\lambda}_{kk^{\prime}ll^{\prime}}$ in the quantum Hamiltonian (1). Note that one- and two-body parts of the Hamiltonian yield terms $\sim R^{2}$ and $\sim R^{4}$, respectively, and that $(R_{0}^{2}-R^{2})^{\frac{1}{2}}$ replaces both creation and annihilation operators of the $b_{0}$ boson. Values $R\in(R_{0},\infty)$ (outside ${\bf\Omega}$) are forbidden and one can visualize this by embedding the system into a sharp energy well yielding the energy $E^{\lambda}(R,\boldsymbol{\phi})$ equal to (21) for ${R\in[0,R_{0}]}$ and to infinity for ${R\in(R_{0},\infty)}$.

On the hypersurface ${\bf\partial\Omega}$, the energy (21) takes values

[TABLE]

which depend only on angles $\boldsymbol{\phi}$. It turns out that the boundary generates a spectral singularity if there exists a point $\boldsymbol{\phi}^{\lambda}_{0}$ such that ${\frac{\partial}{\partial\phi_{l}}E^{\lambda}_{\bf\partial\Omega}(\boldsymbol{\phi})|_{\boldsymbol{\phi}=\boldsymbol{\phi}^{\lambda}_{0}}=0}$ for all ${l=1,...,2f-1}$. At this point, whose energy is ${E^{\lambda}_{\bf\partial\Omega}(\boldsymbol{\phi}^{\lambda}_{0})\equiv E^{\lambda}_{\rm c}}$, the phase-space integration in Eq. (15) yields a non-analytic contribution to the level density $\overline{\rho}^{\lambda}(E)$. We stress that the point $(R,\boldsymbol{\phi})\!=\!(R_{0},\boldsymbol{\phi}^{\lambda}_{0})$ does not have to be a stationary point of the whole Hamiltonian as we set no constraint on $\frac{\partial}{\partial R}E^{\lambda}(R,\boldsymbol{\phi})$.

For the sake of concreteness, let us consider a local extreme (minimum or maximum) of the surface energy function $E^{\lambda}_{\bf\partial\Omega}(\boldsymbol{\phi})$. The extreme is assumed to be separable, i.e., expressed through ${E^{\lambda}_{\bf\partial\Omega}\propto\sum_{l=1}^{2f-1}(\delta r_{l})^{K_{l}}}$, where $\delta r_{l}$ denote some local orthogonal coordinates on ${\bf\partial\Omega}$ in a vicinity of $\boldsymbol{\phi}^{\lambda}_{0}$. The powers $K_{l}$ specify energy dependencies in the respective directions: ${K_{l}=2}$ indicates a quadratic extreme, ${K_{l}=4}$ quartic etc. ${K_{l}\to\infty}$ leads to a totally flat (locally constant) dependence. We introduce an average inverse power $K$ given by: ${K^{-1}=\sum_{l=1}^{2f-1}K_{l}^{-1}/(2f-1)}$. It can be shown (the analysis is not presented here) that as the energy $E$ crosses the energy $E^{\lambda}_{\rm c}$ of the extreme, the smoothed level density close to $E^{\lambda}_{\rm c}$ varies according to

[TABLE]

where $\varrho^{\lambda}(E)$ is an unspecified analytic function and ${\sigma={\rm sgn}\frac{\partial}{\partial R}E^{\lambda}|_{(R,\boldsymbol{\phi})=(R_{0},\boldsymbol{\phi}_{0}^{\lambda})}}$. The sign $+$ or $-$, respectively, corresponds to a local minimum or maximum of the surface energy function and $M\!\in\!\{\frac{1}{2},1,\frac{3}{2},2\}$ stands for the lowest power in the dependence of $E^{\lambda}$ on ${(R-R_{0})}$ at ${\boldsymbol{\phi}=\boldsymbol{\phi}^{\lambda}_{0}}$. This indicates that $\overline{\rho}^{\lambda}(E)$ exhibits at $E^{\lambda}_{\rm c}$ a non-analyticity connected with its $\lceil I\rceil$th derivative, where $\lceil x\rceil$ denotes the ceiling function. The derivative is discontinuous if $I$ is integer, or divergent otherwise. Note that for ${K=2}$ and ${M=\frac{1}{2}}$ (assumingly the most generic case as it does not require an accidental vanishing of the lowest order expansion coefficients in any of the energy dependences) we expect a divergent derivative $\frac{\partial^{f+1}}{\partial E^{f+1}}\overline{\rho}^{\lambda}(E)|_{E=E^{\lambda}_{\rm c}}$. On the other hand, the most “conspirative” case ${K=\infty}$ and ${M\geq 1}$ leads to discontinuous or divergent $\overline{\rho}^{\lambda}(E)|_{E=E^{\lambda}_{\rm c}}$, independently of $f$.

The analysis of other types of stationary points of the boundary function $E^{\lambda}_{\bf\partial\Omega}(\boldsymbol{\phi})$ needs to be performed case by case. Note that the energy $E^{\lambda}_{\bf\partial\Omega}(\boldsymbol{\phi}_{0})$ may coincide with that of a regular stationary point inside ${\bf\Omega}$, which enhances the corresponding singularity in the spectrum.

II.2.4 Finite-size effects due to partial separability

In part II ref:E2 we have shown that finite-size effects in the oscillatory component $\widetilde{\rho}^{\lambda}(E)$ of the level density can be particularly strong in ${f=2}$ systems that exhibit partial (effective) decoupling of both degrees of freedom. If one degree of freedom (mode of motion) in such a decoupled system is much faster than the other, the quantum energy spectrum consists of separated bands of states corresponding to distinguished excitations of both modes. Excitations of the fast mode form a sequence of levels with large spacings, each of them being a base for excitations of the slow mode that give rise to repeated bands of levels with small spacings. If, in such a situation, one of the modes undergoes an ESQPT, we observe a peculiar pattern of finite-size precursors that carry strong features of the underlying ${f=1}$ dynamics. A qualitative analysis of these phenomena has been presented in terms of an effective Hamiltonian describing one-dimensional dynamics based on a boosted momentum of the other degree of freedom, see formula (28) in part II ref:E2 .

In a general interacting boson system, these considerations can be put on a more rigorous ground with the aid of generalized projective coherent states, which are outlined in Appendix A. Suppose that the fast mode of motions of the system, that is at least partially separated from the other(s), is associated with the boson creation operator $\hat{B}^{{\dagger}}(\tilde{\alpha}_{\perp})$. This boson plays the role of a single-mode excitation quantum and must be perpendicular to the ground-state condensate boson, represented by operator $\hat{B}^{{\dagger}}(\tilde{\alpha})$ from Eq. (10). Therefore, in absence of other types of excitations, the states of the system with the total number of bosons $N$ have the form $\propto\hat{B}^{{\dagger}}(\tilde{\alpha}_{\perp})^{N_{\perp}}\hat{B}^{{\dagger}}(\tilde{\alpha})^{N\!-\!N_{\perp}}|0\rangle$, where $N_{\perp}$ is the number of excitation quanta.

We know that parameters $\tilde{\alpha}$ from the expansion of $\hat{B}^{{\dagger}}(\tilde{\alpha})$ are associated with coordinates and momenta of the system. The classical counterpart associated with a quantum Hamiltonian $\hat{H}^{\lambda}$ is expressed solely in terms of the ground-state boson $\hat{B}^{{\dagger}}(\tilde{\alpha})$, namely as ${{\cal H}^{\lambda}(\tilde{\alpha},N)=\langle 0|\hat{B}(\tilde{\alpha})^{N}\hat{H}^{\lambda}\hat{B}^{{\dagger}}(\tilde{\alpha})^{N}|0\rangle/(N!N)}$, see Eq. (11). An excited energy surface, which underlies motions of any character on top of the state with $N_{\perp}$ excitation quanta of the above type, can be written as:

[TABLE]

where ${|\tilde{\alpha}_{\perp},N_{\perp}\rangle=\frac{1}{\sqrt{N_{\perp}!}}\hat{B}^{{\dagger}}(\tilde{\alpha}_{\perp})^{N_{\perp}}|0\rangle}$. Since coefficients $\tilde{\alpha}_{\perp}$ depend via the orthogonality constraint on $\tilde{\alpha}$, the above expression is written as a function of only the latter variables, i.e., standard coordinates and momenta. Formula (24) is a special case of Eq. (37). It will be applied below in the analysis of oscillatory structures in the spectrum of the IBM.

III The IBM as a case study

III.1 First-order QPT Hamiltonian

In what follows, we explore the general ESQPT signatures in the the interacting boson model of collective nuclear dynamics ref:IBMbook . The model, describing rotations and quadrupole vibrations of nuclei, is formulated in terms of bosons $s$ and $d$ with angular momenta 0 and 2, respectively, and has ${f=5}$ degrees of freedom. Considering bilinear combinations of boson creation and annihilation operators, one builds the U(6) dynamical group which branches into three dynamical-symmetry chains starting with subgroups U(5), O(6) and SU(3). All Hamiltonians are rotationally invariant under the group O(3), which ensures conservation of the total angular momentum quantum number $L$. The conserved total number of bosons $N$, resulting from selection of a single irrep of U(6), can undergo the limiting process ${N\to\infty}$, which leads to semiclassical description and sharpening of quantum critical effects.

The IBM was instrumental in setting the foundations for both QPT and ESQPT phenomena ref:Diep80 ; ref:Feng81 ; ref:Cejn09 . The original papers on ESQPTs ref:Cejn06 ; ref:Hein06 ; ref:Mace06 considered the transition between spherical and axially-unstable quadrupole shapes of the ground state [U(5) and O(6) dynamical symmetries]. In this transition, the model is fully integrable all the way and exhibits a second-order ground-state QPT. The seniority quantum number $\tau$ is conserved and the constraint ${\tau=L=0}$ reduces the number of effective degrees of freedom to ${f=1}$. This makes the ESQPT signatures very well visible—they affect in the same way all individual energy levels satisfying the above constraint and are apparent already in the zeroth derivative of the smoothed level density and slope.

Here we focus on the transition between spherical and axially-deformed ground-state shapes, which requires a more complex treatment. The irrotational constraint $L=0$ yields the effective number of degrees of freedom ${f=2}$ (no additional conserved quantity except $N$ exists), which implies that the ESQPT spectral singularities appear in the first derivative of the smoothed level density and slope. Studying the ESQPTs accompanying the first-order ground-state QPT, we parallel the ESQPT analysis done in parts I and II ref:E1 ; ref:E2 and extend it into the setting of a realistic interacting many-body system.

The structural evolution between spherical and axially-deformed equilibrium configurations of the IBM was studied from the ESQPT perspective in Refs. ref:Zhan16 ; ref:Zhan17 . The Hamiltonian used there was of the well known form

[TABLE]

depending on parameters ${\eta\in[0,1]}$ and ${\chi\in[-\frac{\sqrt{7}}{2},+\frac{\sqrt{7}}{2}]}$. Here, ${\hat{n}_{d}=(\hat{d}^{{\dagger}}\cdot\hat{\tilde{d}})}$ is the $d$-boson number operator and ${\hat{Q}^{\chi}_{\mu}=[\hat{d}^{{\dagger}}\hat{s}+\hat{s}^{{\dagger}}\hat{\tilde{d}}]^{(2)}_{\mu}+\chi[\hat{d}^{{\dagger}}\hat{\tilde{d}}]^{(2)}_{\mu}}$ the quadrupole operator. For ${\chi=-\frac{\sqrt{7}}{2}}$, ${\chi=+\frac{\sqrt{7}}{2}}$ or ${\chi=0}$, the Hamiltonian depending on $\eta$ describes the U(5)–SU(3), U(5)–$\overline{{\rm SU(3)}}$ or U(5)–O(6) transition, respectively, while for ${\eta=1}$, the Hamiltonian depending on $\chi$ describes the $\overline{{\rm SU(3)}}$–O(6)–SU(3) transition. Note that $\overline{{\rm SU(3)}}$ differs from SU(3) by the type of equilibrium deformation (oblate instead of prolate). We use the symbol ${[\hat{A}^{(\ell_{1})}\hat{B}^{(\ell_{2})}]^{(\ell)}_{m}}$ for the coupling of spherical tensors with angular momenta $\ell_{1}$ and $\ell_{2}$ to one with angular momentum $\ell$ and projection $m$, while the scalar product is written as ${(\hat{A}^{(\ell)}\cdot\hat{B}^{(\ell)})=(-1)^{\ell}\sqrt{2\ell+1}[\hat{A}^{(\ell)}\hat{B}^{(\ell)}]^{(0)}_{0}}$. The tilde above an annihilation operator means ${\hat{\tilde{x}}^{(\ell)}_{m}=(-1)^{\ell+m}\hat{x}^{(\ell)}_{-m}}$, ensuring proper transformation properties under rotations.

In this paper, we employ a more sophisticated Hamiltonian adopted from the formalism resolving the intrinsic and collective motions within the IBM, see Refs. ref:Kirs85 ; ref:Levi87 . This Hamiltonian enables us to generate a ground-state shape-phase transition with an adjustable barrier between the spherical and deformed configurations ref:Levi06 ; ref:Mace11 ; ref:Levi12 ; ref:Mace14 , which is in contrast to the “consistent-$Q$” Hamiltonian (25) yielding a very small barrier. In particular, we use two Hamiltonian forms,

[TABLE]

where $\zeta\in[0,1]$ and $\xi\in[0,1]$ are two control parameters and $\hat{S}^{{\dagger}}(\beta_{0}^{\prime})$ and $\hat{D}^{{\dagger}}(\beta_{0}^{\prime},\zeta)$ two specific boson pair operators with ${\ell=0}$ and 2:

[TABLE]

$\beta_{0}^{\prime}$ is an adjustable parameter that modifies the shape of the potential energy functions associated with Hamiltonians (LABEL:H1n) and (27). It affects the size of the equilibrium deformation, rigidity of spherical and deformed potential wells, and the height of the barrier between these configurations in the phase-transitional domain. Note that for the purposes of the present paper we modify the original notation of Refs. ref:Mace11 ; ref:Levi12 ; ref:Mace14 , where the control parameter of Eq. (LABEL:H1n) is denoted as ${\rho\equiv\zeta/\beta_{0}^{\prime}}$ while operators ${\hat{P}_{0}^{{\dagger}}(\beta_{0}^{\prime})\equiv\hat{S}^{{\dagger}}}(\beta_{0}^{\prime})$, ${\hat{R}_{2\mu}^{{\dagger}}(\rho)\equiv\hat{D}_{\mu}^{{\dagger}}(\beta_{0}^{\prime},\zeta=\beta_{0}^{\prime}\rho)/\beta_{0}^{\prime}}$ and ${\hat{P}^{{\dagger}}_{2\mu}(\beta_{0}^{\prime})\equiv\hat{D}_{\mu}^{{\dagger}}(\beta_{0}^{\prime},\zeta=1)}$ are employed instead of those in Eqs. (28) and (29).

Hamiltonian ${\hat{H}_{1}(\beta_{0}^{\prime},\zeta=0)}$ represents the U(5) limit of the model, while ${\hat{H}_{2}(\beta_{0}^{\prime}=\sqrt{2},\xi=1)}$ corresponds to the SU(3) limit. The first-order critical-point Hamiltonian between these two dynamical symmetries reads ${\hat{H}_{\rm c}(\beta_{0}^{\prime})\equiv\hat{H}_{1}(\beta_{0}^{\prime},\zeta=1)=\hat{H}_{2}(\beta_{0}^{\prime},\xi=0)}$. Hence the transition from the U(5) limit to the critical point is realized by Hamiltonian (LABEL:H1n) with $\zeta$ varying from 0 to 1, and the subsequent transition from the critical point to the SU(3) limit is driven by Hamiltonian (27) with $\xi$ changing from 0 to 1 and $\beta_{0}^{\prime}$ converging to $\sqrt{2}$. In agreement with the notation used in Sec. II, we introduce a single control parameter $\lambda\in[0,\infty)$ and define the Hamiltonian as follows:

[TABLE]

Since $\beta_{0}^{\prime}$ is not necessarily restricted to $\sqrt{2}$ and $\lambda$ may exceed the value 2 (corresponding to ${\xi=1}$), we do not insist on the SU(3) dynamical symmetry as an endpoint of the transition. We stress that although both parts of this formula coincide at the critical point ${\lambda=1\equiv\lambda_{\rm c}}$, implying that all eigenstates are the same there, the slopes $\frac{\partial}{\partial\lambda}E_{i}^{\lambda}$ of excited levels differ as one approaches to $\lambda_{\rm c}$ from the right and from the left. The energy of the ground state is set to ${E_{0}^{\lambda}=0}$ independently of $\lambda$.

The Hamiltonian described above was studied in detail from the viewpoint of classical and quantum chaos and with respect to partial and quasi dynamical symmetries ref:Mace11 ; ref:Levi12 ; ref:Mace14 . Its QPT and ESQPT analysis relies upon the construction of the classical limit. With the aid of the Glauber coherent-state formalism in the form presented in Ref. ref:Hatc82 (cf. Sec. II.1 above) the classical counterparts of Hamiltonians (LABEL:H1n) and (27) constrained to $L=0$ can be cast in the following way,

[TABLE]

where $T=p_{\beta}^{2}+\beta^{-2}p_{\gamma}^{2}$ denotes the standard kinetic term and ${\cal H}_{d,0}=\tfrac{1}{2}(T+\beta^{2})$ stands for the classical limit of $\hat{n}_{d}$ restricted to $L=0$.

Polar coordinates ${\beta\in[0,\sqrt{2}]}$ and ${\gamma\in[0,2\pi)}$ are related to the Bohr deformation variables, and can be represented by equivalent Cartesian coordinates ${x=\beta\cos\gamma}$ and ${y=\beta\sin\gamma}$. Together with the canonically conjugate momenta ${p_{\beta}\in[0,\sqrt{2}]}$ and ${p_{\gamma}\in[0,1]}$, or ${p_{x}=p_{\beta}\cos\gamma-\frac{p_{\gamma}}{\beta}\sin\gamma\in[0,\sqrt{2}]}$ and ${p_{y}=p_{\beta}\sin\gamma+\frac{p_{\gamma}}{\beta}\cos\gamma\in[0,\sqrt{2}]}$, they span a compact classical phase space domain ${\bf\Omega}$, in agreement with Eq. (7). The Hamiltonian exhibits the well known discrete symmetry under the shift ${\gamma\to\gamma+\frac{2\pi}{3}}$. Note that the present bounded form of the deformation parameter $\beta$ can be transformed to a boundless form ${\beta^{\prime}=\beta/\sqrt{2-\beta^{2}}\in[0,\infty)}$ ref:Klei81 , which is more common in some literature. The value ${\beta_{0}^{\prime}=\sqrt{2}}$ of the shape parameter is associated with the equilibrium deformation ${\beta_{0}=2/\sqrt{3}}$ in the SU(3) limit.

III.2 ESQPT analysis

III.2.1 Smooth level density

Results of our numerical analysis of the smooth component $\overline{\rho}^{\lambda}(E)$ of the level density associated with the IBM Hamiltonian (30) are shown in Figs. 1–4. We use two values of the Hamiltonian shape parameter $\beta_{0}^{\prime}$, namely ${\beta_{0}^{\prime}=\sqrt{2}}$ in Figs. 1 and 3, and ${\beta_{0}^{\prime}=1.7}$ in Figs. 2 and 4. Both cases represent the first-order QPT between spherical and deformed ground-state shapes with the same critical value ${\lambda_{\rm c}=1}$. The phase-coexistence interval $(\lambda^{*},\lambda^{**})$, demarcated by the spinodal and antispinodal points, coincides with ${(0.707,1.333)}$ for ${\beta_{0}^{\prime}=\sqrt{2}}$ and with ${(0.460,1.257)}$ for ${\beta_{0}^{\prime}=1.7}$. As explained above, the U(5) dynamical symmetry is located at ${\lambda=0}$ for all $\beta_{0}^{\prime}$ values, and the point ${\lambda=2}$ of the ${\beta_{0}^{\prime}=\sqrt{2}}$ case is associated with the SU(3) dynamical symmetry. Samples of the potential energy functions ${\cal V}^{\lambda}(x,y)$ along the ${y=0}$ line can be seen in the left columns of Figs. 3 and 4.

Figures 1–4 depict not the level density itself, but its first derivative with respect to energy, $\frac{\partial}{\partial E}\overline{\rho}^{\lambda}(E)$. This is in agreement with the general classification of ESQPTs in ${f=2}$ systems, where we generically expect singularities in the first derivative of the level density and slope—see Sec. II.2.1. The smooth level density is calculated from the integral formula (15) using the classical Hamiltonian forms (31) and (32). The integration is performed numerically in the whole above-defined finite phase space of the system.

The most complete information on the level density is contained in Figs. 1 and 2. The top panels of these figures show $\frac{\partial}{\partial E}\overline{\rho}^{\lambda}(E)$ in the plane of control parameter ${\lambda\in[0,3.2]}$ and excitation energy ${E\in[0,3]}$ encoded in the shades of blue and red colors. Blue color demarcates the domains of increasing level density (positive derivative), red indicates the domains of decreasing level density (negative derivative). The darker is the shade, the larger absolute value takes the derivative at the corresponding place. White areas represent the domains where the level density is flat (has zero derivative). In these domains, $\overline{\rho}^{\lambda}(E)$ is either equal to zero, which is the case beyond the finite energy domain $[E^{\lambda}_{\rm min},E^{\lambda}_{\rm max}]$ of the spectrum (i.e., below the blue and above red areas), or is fixed at a non-zero constant, which happens in some regions in the middle of the spectrum (between the blue and red areas).

The middle panels of both Figs. 1 and 2 display critical curves where the level density derivative from the respective upper panel exhibits any kind of singularity. These curves, which form the ESQPT critical borderlines in the $\lambda\times E$ plane, correspond to energies of various stationary points of classical dynamics. The stationary points were determined by an exhaustive numerical procedure using the classical Hamiltonian forms (31) and (32). Note that analytic determination was possible in some cases for ${\beta_{0}^{\prime}=\sqrt{2}}$, see Ref. ref:Mace14 .

The dependencies of $\frac{\partial}{\partial E}\overline{\rho}^{\lambda}(E)$ on energy $E$ for four selected values of $\lambda$ are shown in the right columns of Figs. 3 and 4. The selected values correspond to the cuts (a)–(d) indicated in Figs. 1 and 2. In the level density plots shown in the respective panels (a)–(d) of Figs. 3 and 4 we observe upward and downward jumps as well as positive and negative divergences of the level density derivative. These singularities appear precisely at the critical borderlines from the middle panels of Figs. 1 and 2.

Classification of the ESQPT singularities in Figs. 1–4 is encoded in the types of the corresponding lines. The same line type is used consistently in the middle and bottom panels of Figs. 1 and 2, as well as in Figs. 3 and 4. The key is as follows:

• (i)

The full blue line corresponds to stationary points with index ${r=0}$ (minima of the Hamiltonian function). These cause upward jumps of the level density derivative.

• (ii)

The dashed blue lines indicate stationary points with index ${r=1}$ (saddles of the Hamiltonian function). They give rise to a positive logarithmic divergence of the level density derivative.

• (iii)

The dashed-dotted black line represents stationary points with index ${r=2}$ (saddles of the Hamiltonian function). These yield downward jumps of the level density derivative.

• (iv)

The dashed red line demarcates stationary points with index ${r=3}$ (saddles of the Hamiltonian function). Such points result in a negative logarithmic divergence of the level density derivative.

• (v)

The full red line corresponds to stationary points with index ${r=4}$ (maxima of the Hamiltonian function). They cause upward jumps of the level density derivative, similar to type (i).

• (vi)

The thick green lines are connected with the boundary $\partial{\bf\Omega}$ of the phase space. The corresponding singularities are divergences or jumps of the level density derivative.

It turns out that most of the above level-density singularities follow from stationary points of the respective potential energy functions ${\cal V}^{\lambda}(x,y)$ (cf. the left panels of Figs. 3 and 4). These form essentially the same types of ESQPTs as studied in parts I and II ref:E1 ; ref:E2 . However, some of the singularities result from the new features discussed in the present paper, either from the nontrivial kinetic energy term of the Hamiltonian (see Sec. II.2.2), or from the boundedness of the compact phase space domain (see Sec. II.2.3).

The “kinetic” type of ESQPTs can be recognized with the aid of the bottom panels of Figs. 1 and 2, which show the coordinate $x$ and momentum $p_{x}$ corresponding to the stationary points responsible for the ESQPTs in the respective middle panels ($y$ and $p_{y}$ are not shown as they contain similar information). Nonzero momentum values indicate the “kinetic” stationary points satisfying Eq. (20) with ${\boldsymbol{p}_{0}\neq 0}$. This concerns a segment of the line (iv) in Fig. 1, and a segment of line (iii) and whole lines (ii) and (iv) in Fig. 2.

The “boundary” type of ESQPTs is connected with the critical borderlines in the middle panels of Figs. 1 and 2 which have no counterparts in the bottom panels. The lines of type (vi) correspond to the overall minimum $E_{{\bf\partial\Omega}{\rm min}}^{\lambda}$ and maximum $E_{{\bf\partial\Omega}{\rm max}}^{\lambda}$ of the boundary energy function (22), which are independently of $\beta_{0}^{\prime}$ given by the following simple analytic prescriptions:

[TABLE]

We note that the associated singularities of the level density (visible already in the first derivative, see Figs. 3 and 4) are beyond the classification following from Eq. (23) because of a high degree of degeneracy and non-separability of these points. Additional singularities resulting from less degenerate stationary points of $E_{{\bf\partial\Omega}}^{\lambda}(\boldsymbol{\phi})$ may be present in higher level density derivatives—their detection would be already beyond our numerical reach.

Comparing the results for ${\beta_{0}^{\prime}=\sqrt{2}}$ (Figs. 1 and 3) and ${\beta_{0}^{\prime}=1.7}$ (Figs. 2 and 4), we see an overall similarity of phase diagrams, but note that the latter case exhibits a somewhat more complex system of ESQPTs. The ${\beta_{0}^{\prime}=1.7}$ energy surface generates a larger number of critical borderlines and a more sophisticated bifurcation pattern of stationary points. In particular, we observe three ESQPTs associated with the “kinetic” type of stationary points, instead of just one such ESQPT present in the ${\beta_{0}^{\prime}=\sqrt{2}}$ phase diagram. Our numerical study of a wider range of parameter $\beta_{0}^{\prime}$ showed that complexity of the phase diagram keeps increased within the whole interval ${\beta_{0}^{\prime}\in(\sqrt{2},\sqrt{3})}$. If $\beta_{0}^{\prime}$ decreases below $\sqrt{2}$, the ESQPTs at $(\lambda,E)=(0,1)$ become degenerate (see Fig. 1), which simplifies the picture at small values of $\lambda$. On the other hand, if $\beta_{0}^{\prime}$ increases above $\sqrt{3}$, the spinodal point $\lambda^{*}$ disappears (implying that the spherical and deformed potential minima coexist at ${\lambda=0}$), which again leads to a less complicated structure of critical borderlines.

III.2.2 Oscillatory level density

While the smooth component $\overline{\rho}^{\lambda}(E)$ of the level density captures the constitutive features that define the ESQPT of a given type in the infinite-size limit, the oscillatory component $\widetilde{\rho}^{\lambda}(E)$ carries information on finite-size effects. The upper panel of Fig. 5 displays the spectrum of ${L=0}$ eigenstates of the IBM Hamiltonian (30) with ${N=50}$ for ${\beta_{0}^{\prime}=\sqrt{2}}$ and ${\lambda\in[0,2.2]}$. This is a finite-$N$ realization of the semiclassical level density discussed in the previous paragraph, cf. Fig. 1. The lower panel of Fig. 5 shows the oscillatory component of the level density derived from the spectrum in the upper panel. It was calculated from the formula

[TABLE]

where $\delta_{\sigma}(x)$ stands for a Gaussian of a width $\sigma$ smaller than a typical spacing between energy levels $E_{i}^{\lambda}$, so ${\sigma<1/\overline{\rho}^{\lambda}(E)}$, and the smooth component $\overline{\rho}^{\lambda}(E)$ is evaluated from the semiclassical integral in Eq. (15). We observe a regular alternation of denser and sparser parts of the spectrum, indicated by darker and lighter shades, respectively.

It turns out that the structures in Fig. 5 exhibit the separability-enhanced precursors of ESQPTs discussed in Sec. II.2.4. In part II of this series ref:E2 we showed, using a simpler model with partially separated $x$ and $y$ degrees of freedom, that stationary points of an effective potential for $x$ (obtained by fixing the averages of variables related to the $y$ degree of freedom) give rise to precursors of ESQPT structures in spectra of excitations in $x$ direction built on individual excitations in $y$ direction. The effective number of degrees of freedom associated with these structures was ${f_{\rm eff}=1}$ rather than ${f=2}$.

A similar effect can be demonstrated in the present model, in which the two partially separated degrees of freedom are the $\beta$ and $\gamma$ vibrational modes. For the ${\beta_{0}^{\prime}=\sqrt{2}}$ Hamiltonian (30) it turns out that the energy ratio between the $\beta$ and $\gamma$ phonon energies is ${\varepsilon_{\beta}/\varepsilon_{\gamma}=(2\lambda-1)/3}$ ref:Mace14 . So the $\gamma$ mode is markedly faster than the $\beta$ mode for $\lambda\lesssim 1.5$, and in this parameter domain we expect bands of $\beta$-vibrational states built on individual $\gamma$-vibrational states in the spectrum. Hence the formula (24) can be applied with the condensate boson $\hat{B}^{{\dagger}}(\tilde{\alpha})$ from Eq. (38) and the orthogonal excitation boson $\hat{B}^{{\dagger}}(\tilde{\alpha}_{\perp})$ associated with $\hat{B}^{{\dagger}}_{\gamma}(\gamma)$ from Eq. (39). As the coefficients in these equations have no imaginary components (no momenta), the excited energy surfaces resulting from this construction represent surfaces of potential energy depending on coordinates $\beta$ and $\gamma$, or equivalently $x$ and $y$.

Excited potentials (their ${y=0}$ cuts) for three choices of parameter $\lambda$ and several numbers $N_{\gamma}$ of the $\gamma$-vibrational quanta are shown in Fig. 6 (the reason why $N_{\gamma}$ takes only even values is explained in Appendix A). The values of $\lambda$ are chosen at the spinodal point ${\lambda^{*}=\frac{1}{\sqrt{2}}}$ (the point, where the ground-state potential develops a secondary minimum corresponding to a deformed shape), critical point ${\lambda_{\rm c}=1}$ (where both minima of the ground-state potential become degenerate), and the antispinodal point ${\lambda^{**}=\frac{4}{3}}$ (where the minimum associated with the spherical shape disappears). We observe that with increasing $N_{\gamma}$ the minimum corresponding to the spherical shape becomes deeper in comparison with the minimum associated with the deformed shape. As a consequence, the phase evolution of excited states is retarded with respect to the ground state (the critical value of parameter $\lambda$ at which both minima of the excited potential get degenerate increases above the ground-state critical point $\lambda_{\rm c}$). Although we do not evaluate the dependence of the excited potential energy surfaces on the $y$ variable, we assume that for lower excitation energies the ${y=0}$ cut represents the most substantial information due to expected averaging of excited dynamics at ${\langle\gamma\rangle=0}$.

Energies of the lowest stationary points of the excited potentials from Fig. 6 are displayed by the colored curves in the lower panel of Fig. 5 for the whole interval ${\lambda\in[0,2.2]}$. Each line type corresponds to the same number $N_{\gamma}$ of the $\gamma$-excitation quanta (the color code is the same as in Fig. 6). The three curves of the same line type with increasing energy represent, consecutively, the primary minimum, secondary minimum and intermediate maximum (probably a saddle point of the whole dependence on $x$ and $y$) of the respective excited potential.

We see that at lower excitation energies there is a rather strong correlation of the colored curves with the oscillatory structures in Fig. 5. This illustrates the above-discussed finite-size effects induced by partial separability. In general, minima of an effective ${f_{\rm eff}=1}$ potential lead to a step-like increases of the level density while a local maximum indicates a peak in the level density. Such structures can be indeed distinguished in Fig. 6. We clearly observe their gradual shift to higher parameter values, which agrees with the above-discussed evolution of excited potentials. For higher energies, the correspondence gradually deteriorates. This indicates that the large-$N$ approximation involved in the foundations of the present generalized coherent state method becomes worse as the ratio $N_{\gamma}/N$ increases. Moreover, for highly excited potentials the whole ${(x,y)}$ dependence may become increasingly relevant.

We point out that the spectra of our model contain, besides the structures in $\widetilde{\rho}^{\lambda}(E)$ related to partial separability, also some additional oscillatory structures with a different origin. These can be noticed in both panels of Fig. 5 in the domain ${\lambda>1}$ for energies ${E<\lambda-1}$, mostly close to the values of $\lambda$ for which the energy ratio ${\varepsilon_{\beta}/\varepsilon_{\gamma}}$ becomes a simple rational number (particularly ${1,\frac{2}{3},\frac{1}{2},\frac{2}{5},\frac{1}{3}}$). It was shown in Ref. ref:Mace14 (cf. Fig. 14 therein) that these clustering effects are related to splitting of classical phase-space tori near the above resonances via a mechanism following from the classical chaos theory.

III.2.3 Links to previous results

To conclude the present ESQPT case study, let us briefly address the relation of our results to those of Ref. ref:Zhan16 . The cited work applied the simplified IBM Hamiltonian (25) along several symmetry-connecting paths in parameters $\eta$ and $\chi$. It turned out that energy dependencies of expectation values of the $d$-boson number ${\hat{n}_{d}}$ in individual excited states with ${L=0}$ exhibit some special structures near the ESQPT borderlines, considered in Ref. ref:Zhan16 only for cases of the potential maximum related to spherical deformation and the saddle point related to oblate deformation.

These structures can be justified in terms of the theoretical background outlined in Sec. II.2.1. In particular, we saw that the ESQPT singularities of the level density are accompanied by the same types of singularities in the level flow, i.e., for ${f=2}$ systems generically by nonanalyticities of the smoothed slope of the spectrum. For the Hamiltonian (25) it follows from Eq. (14) that

[TABLE]

where $|\psi_{i}^{(\eta,\chi)}\rangle$ is the $i$th eigenstate with energy $E_{i}^{(\eta,\chi)}$ at the given point $(\eta,\chi)$. Since energy derivatives on the left-hand side of these formulas (if averaged out over the ${L=0}$ eigenstates within a narrow interval around the given energy) exhibit precursors of the ESQPT singularities, the same must be true for the expectation values on the right-hand side. The latter involve the ${\hat{n}_{d}}$ expectation value, either explicitly as in the first formula, or via expansions of operators $(\hat{Q}^{\chi}\!\cdot\!\hat{Q}^{\chi})$ and $\frac{\partial}{\partial\chi}(\hat{Q}^{\chi}\!\cdot\!\hat{Q}^{\chi})$ into the Casimir invariants of all subalgebras of U(6); cf. formula (5) and Table I in Ref. ref:Cejn98 . Hence the locally averaged $\hat{n}_{d}$ expectation values naturally show some anomalous energy dependencies near the ESQPT critical energies.

It shall be pointed out that anomalous behavior at ESQPT energies may be anticipated for locally smoothed expectation values of almost any generic observable $\hat{A}$. Indeed (returning to a general Hamiltonian $\hat{H}^{\lambda}$), since the slope ${\frac{\partial}{\partial\lambda}E_{i}^{\lambda}}$ of individual levels can be decomposed into a sum of terms ${\langle\psi_{i}^{\lambda}|\hat{A}|\psi_{i}^{\lambda}\rangle}$ and ${\langle\psi_{i}^{\lambda}|\frac{\partial}{\partial\lambda}\hat{H}^{\lambda}-\hat{A}|\psi_{i}^{\lambda}\rangle}$, it is natural to assume that the anomaly in the smoothed slope is shared by both these terms. Note that instead of spectral averages of expectation values one can consider the whole Peres lattices ${E_{i}^{\lambda}\times\langle\psi_{i}^{\lambda}|\hat{A}|\psi_{i}^{\lambda}\rangle}$ constructed for all Hamiltonian eigenstates. Examples (from another model) of ESQPT anomalies simultaneously reflected by Peres lattices of several observables can be found in Re. ref:Kloc17a .

The present analysis of ESQPT effects in the IBM differs from one of Ref. ref:Zhan16 in several aspects. First, we focus on the description of ESQPTs in terms of the level density instead of expectation values. Peres lattices of observable ${\propto\sqrt{\hat{n}_{d}}}$ for the present Hamiltonian were studied in Ref. ref:Mace14 . Second, our aim was to localize and classify spectral effects of all singularities in classical dynamics and to anchor them in the general theory developed in parts I and II ref:E1 ; ref:E2 and Sec. II. Using this exhaustive treatment, we demonstrate that for a genuine ${f=2}$ many-body system the singularities appear in large numbers and are of rather diverse types. Third, we used the more complex IBM Hamiltonian from Eq. (30), which describes the first-order transition between spherical and deformed nuclear shapes with flexible properties. This allowed us to continue our previous study of ESQPT structures accompanying a generic first-order ground-state QPT.

Note that Ref. ref:Zhan17 extends the study of Ref. ref:Zhan16 to states with ${L>0}$. However, relaxing the angular momentum constraint, one generally leaves the ${f=2}$ subset of IBM dynamics. So the ${L>0}$ results are beyond the reach of the present approach. In particular, for ${f>2}$ systems one generally expects ESQPT singularities in higher derivatives of the level density and flow.

IV Summary and outlook

This paper is part III of a series devoted to the study of ESQPTs in bound quantum systems with ${f=2}$ effective degrees of freedom. In parts I and II ref:E1 ; ref:E2 , we considered systems with standard Hamiltonians of the type ${\hat{H}^{\lambda}=\frac{1}{2}\hat{\boldsymbol{p}}^{2}+V^{\lambda}(\hat{\boldsymbol{q}})}$. Here we have extended our view to genuine quantum many-body systems with Hamiltonians written in terms of particle creation and annihilation operators. In particular, we have studied bound systems composed of several mutually interacting bosonic species, assuming that the total number of particles is conserved. The analysis has been first performed for a general system of the above type, and then illustrated by a concrete example of the interacting boson model of nuclear collective dynamics.

The present work can be seen as a realistic implementation of the results discussed in parts I and II ref:E1 ; ref:E2 , but also as a study of new ESQPT-related effects that emerge in the many-body context. In particular, we have discussed additional singularities of the level density that appear in connection with (i) coordinate-dependent kinetic terms of a generic many-body Hamiltonian and (ii) finiteness of the phase space resulting from the conserved number of particles. We have also shown that finite-size oscillatory structures in the level density, caused by partial separability of dynamics (as thoroughly studied in part II ref:E2 ), can be effectively described by excited energy surfaces resulting from the formalism of generalized coherent states. These structures, which originate in dynamics reduced to a single degree of freedom, can in some cases and for moderate values of $N$ overwhelm the precursors of the ${N\to\infty}$ ESQPT signatures in the smooth component of the level density.

It should be pointed out that one may hardly hope in a direct experimental evidence of the effects discussed here in collective dynamics of nuclei at low energies, where the IBM is applicable. The density of states is too low in these domains to provide an unambiguous signature of an ESQPT, while at higher energies with larger state densities, the number of effective degrees of freedom grows very rapidly in nuclei. On the other hand, the nuclear IBM offers a powerful tool to verify the conclusions of our general analysis, which may turn relevant in other contexts of nuclear physics.

The present series of works can be naturally continued by investigating other types of quantum many-body systems with two effective degrees of freedom. First, the analysis can be done for some bound fermionic systems. If the ${f=2}$ collective dynamics of such a system is described by means of a finite dynamical algebra, there might be a close link to the present work enabled by a suitable bosonic representation of the algebra ref:Klei91 . Second, a similar study as presented here can be performed in the context of extended lattice systems with collective excitations, like two-dimensional crystals ref:Diet13 ; ref:Iach15 ; ref:Diet17 .

Acknowledgments

M.M. thanks F. Iachello for insightful discussions and acknowledges support by the Czech MEYS project RVO:68081731. P.S. and P.C. acknowledge funding by the Charles University Project UNCE/SCI/013. A.L. appreciates support by the Israel Science Foundation.

Appendix A Generalized projective coherent states

As an extension of Eq. (10), we define a generalized coherent state given by the following prescription:

[TABLE]

Here, $\hat{B}^{{\dagger}}(\tilde{\alpha})$ is a boson creation operator from Eq. (9), with ${\{\tilde{\alpha}\equiv\{\alpha_{k}/\lVert\alpha\rVert\}_{k=1}^{n}}$ denoting normalized expansion coefficients, while operators ${\hat{B}_{\perp l}^{{\dagger}}(\tilde{\alpha})\equiv\hat{B}^{{\dagger}}(\tilde{\alpha}_{\perp l})}$ with ${l=1,...,m}$ create normalized bosons in states perpendicular to $\hat{B}^{{\dagger}}(\tilde{\alpha})$ and to each other (expansion coefficients $\tilde{\alpha}_{\perp l}$ satisfy the orthonormality constraints ${\tilde{\alpha}_{\perp l}\cdot\tilde{\alpha}=0}$ and ${\tilde{\alpha}_{\perp l}\cdot\tilde{\alpha}_{\perp l^{\prime}}=\delta_{ll^{\prime}}}$, and therefore depend on $\tilde{\alpha}$). Powers $N_{l}$ represent the numbers of bosons in these states. As the total number of bosons $N$ is fixed, we can write ${N_{0}=N-N_{1}-N_{2}-\cdots-N_{m}}$.

The states (36) may be used in an extended variational procedure ref:Capr05 to approximate the eigensolutions of an interacting boson Hamiltonian $\hat{H}$ conserving the total particle number. The ground state for large $N$ is captured by a condensate $\propto\hat{B}^{{\dagger}}(\tilde{\alpha})^{N}|0\rangle$ from Eq. (10), that is by the state (36) with ${\{N_{0},N_{1},...,N_{m}\}=\{N,0,...,0\}}$. The optimal set of parameters $\tilde{\alpha}$ results from the minimization of the energy surface (11). Excited states are then described with the aid of the supplementary bosons $\hat{B}^{{\dagger}}_{\perp l}(\tilde{\alpha})$, which depend via the orthogonality constraint on the ground-state parameters $\tilde{\alpha}$. These bosons must be selected so that they represent the proper excitation modes of the system. By increasing the ${l>0}$ occupation numbers $\{N_{l}\}$ in Eq. (36), we obtain the bands of excited states of various types.

In analogy to the use of condensate states (10) for the construction of classical Hamiltonian function ${\cal H}(\tilde{\alpha},N)$ in the phase space, see Eq. (11), the generalized coherent states (36) allow one to define classical energy surfaces corresponding to individual excited states. They are given by expressions

[TABLE]

which prescribe specific dependences of energy on coordinates and momenta (encoded in variables $\tilde{\alpha}$) for the states with given numbers $\{N_{l}\}$ of elementary excitations. Let us stress that the surfaces ${\cal H}_{\rm ex}$ are not designed for a variational determination of excited states as those would not in general be orthogonal to the variational ground state and to each other. However, it turns out (see Sec. III.2.2) that stationary points of these surfaces have a strong impact on the oscillatory structures of the level density.

Applying these ideas in the framework of the IBM, we first define the ground-state condensate boson:

[TABLE]

Note that variables ${\beta\in[0,\sqrt{2}]}$ and ${\gamma\in[0,2\pi)}$ correspond to the standard nuclear shape parameters written in the finite-range convention ref:Klei81 . The components with $\hat{d}_{\pm 1}^{{\dagger}}$ bosons vanish due to the transformation to the intrinsic frame ref:IBMbook . Making associations ${\hat{b}^{{\dagger}}_{0}\equiv\hat{s}^{{\dagger}}}$ and $\{\hat{b}^{{\dagger}}_{1},...,\hat{b}^{{\dagger}}_{5}\}\equiv\{\hat{d}_{-2}^{{\dagger}},...,\hat{d}_{+2}^{{\dagger}}\}$, the expansion coefficients in Eq. (10) read as follows: $\tilde{\alpha}=\left(s(\beta),\frac{1}{2}\beta\sin\gamma,0,\frac{1}{\sqrt{2}}\beta\cos\gamma,0,\frac{1}{2}\beta\sin\gamma\right)$. Reality of these coefficients indicates that the function ${\cal H}(\tilde{\alpha},N)$ in Eq. (11) is not the full Hamiltonian, but only the potential energy depending on shape parameters $\beta,\gamma$.

In addition to the condensate boson (38), we introduce bosons associated with $\beta$ and $\gamma$ vibrational modes ref:Levi87 :

[TABLE]

The single-particle states associated with $\hat{B}^{{\dagger}}(\beta,\gamma)$, $\hat{B}^{{\dagger}}_{\beta}(\beta,\gamma)$ and $\hat{B}^{{\dagger}}_{\gamma}(\gamma)$ are normalized and mutually orthogonal. Note that one can introduce yet three other collective bosons, namely ${\tfrac{1}{\sqrt{2}}(\hat{d}^{\dagger}_{+1}\pm\hat{d}^{\dagger}_{-1})}$ and ${\tfrac{1}{\sqrt{2}}(\hat{d}^{\dagger}_{+2}-\hat{d}^{\dagger}_{-2})\equiv\hat{B}^{\dagger}_{z}}$, which together with the above bosons form a complete set replacing the original six bosons $\hat{s}^{{\dagger}}$ and $\{\hat{d}_{m}^{{\dagger}}\}$. In the large $N$ limit, and in absence of an axial symmetry of the intrinsic shape (i.e., for $\gamma$ not equal to integer multiples of $\frac{\pi}{3}$), the additional bosons characterize collective rotations, while bosons (38) and (39) approximate intrinsic excitations. In particular, multiple $\gamma$ excitations can be written as ${\propto\hat{B}^{{\dagger}}_{\gamma}(\gamma)^{N_{\gamma}}\hat{B}^{{\dagger}}(\beta,\gamma)^{N-N_{\gamma}}|0\rangle}$ with ${N_{\gamma}=1,2,3,\dots}$. If the intrinsic shape is axially symmetric around $z$, i.e., ${\gamma=0}$ or $\pi$, the ${\hat{B}^{{\dagger}}_{\gamma}(\gamma)}$ excitation can be combined with $\hat{B}^{\dagger}_{z}$ in such a way that the projection $K$ of angular momentum on the symmetry axis is conserved. The multiple $\gamma$ excitations with ${K=0}$ exist only for $N_{\gamma}$ even and are given by $\propto(\hat{d}^{\dagger}_{+2}\hat{d}^{\dagger}_{-2})^{N_{\gamma}/2}\hat{B}^{\dagger}(\beta,\gamma)^{N-N_{\gamma}}|0\rangle$. It can be shown, however, that the latter form of $\gamma$ excitations yields the same $\beta$ dependeces of excited potentials as the ${\gamma=0}$ or $\pi$ cuts of the potentials obtained with the aid of ${\hat{B}^{{\dagger}}_{\gamma}(\gamma)}$. This is relevant for the explanation of Fig. 5.