\section{Introduction}

The radial-velocity method remains the most versatile technique for
determining the full orbital properties (e.g.~minimum mass,
eccentricity) of extrasolar planets.  More than 500 planets have been
discovered by this method, and considerable efforts are afoot to
understand the true underlying distribution of planetary properties
based on these data \citep{cumming08, shen08, otoole09a, monster,
howard10, foreverpaper, jupiters, etaearth}.  However, there are
significant biases in the measurement of planetary parameters from
sparsely sampled radial-velocity data, especially if the amplitude ($K$)
of the signal is small.  Particularly problematic is the orbital
eccentricity, a quantity which is critically important for understanding
the formation and dynamical history of planetary systems \citep{zhou07,
ford08, shen08}.  Keplerian orbit-fitting algorithms, in a blind
mathematical attempt to minimise $\chi^2$, often resort to increasing
the eccentricity of a model fit.  \citet{shen08} found this introduced
bias toward higher fitted eccentricities to be most egregious for data
with signal-to-noise ratio $K/\sigma\,\ltsimeq$3.  Similarly,
\citet{otoole09a} showed that there is a computational bias against
fitting eccentricities near zero, and that the true uncertainties in
orbital parameters can be 5-10 times larger than the formal
uncertainties emerging from standard least-squares fits.  There is,
therefore, reason to suspect that the observed eccentricity distribution
of radial-velocity detected exoplanets is biased towards higher
eccentricities.

Recently, further radial-velocity monitoring has revealed additional
planets in two systems previously thought to host a single planet
(HD\,142 and HD\,159868: Wittenmyer et al.~2012b), and the best-fit
eccentricities of the previously known planets have significantly
decreased as a result of fitting \textit{two} planetary signals.  This
was most obvious for HD\,159868, where the previously known planet, with
a period of $\sim$1180 days, had an orbit best fit with a very high
$e=0.69\pm0.02$ \citep{otoole07}.  As a result of new data, and a
two-planet solution, HD\,159868b is now best fit with a circular orbit
($e=0.01\pm$0.03).  The possibility that two nearly-circular orbits can
masquerade as a single, eccentric orbit has been explored by
\citet{ang10} and \citet{rh09}.  Motivated by these findings, we now ask
``Which eccentric single-planet systems can be better fit with two
low-eccentricity planets?''

The ambiguities in orbital eccentricity can arise from one of four
degeneracies, which we summarise here.  First, there is the degeneracy
betweeen a single planet on an eccentric orbit and two planets on
circular orbits in a 2:1 configuration \citep{ang10}.  Second is the
degeneracy between one eccentric planet and two co-orbital planets (i.e.
in a 1:1 resonance or ``Trojan pair''), as described in
\citet{laughlin02} and \citet{g12}.  Third is the degeneracy between a
single eccentric planet and a circular planet with a long-period
companion \citep{rh09}.  The fourth degeneracy is that noted above for
HD\,159868: that between a single eccentric planet and two
nearly-circular planets with poorly sampled orbital periods
\citep{142paper}.

This paper is organised as follows: Section 2 describes the data and
analysis procedures, Section 3 gives the results of our efforts to fit
two low-eccentricity planets, and details those systems which had the
most promising results.  Finally, we discuss our conclusions in Section
4.

\section{Data Analysis and Orbit Fitting}

We selected all radial-velocity detected single-planet systems with
publicly available data and published eccentricities $e>0.3$.  This can
be considered a high-eccentricity subset, as the mean eccentricity for
the current population of confirmed planets is 0.22\footnote{Planet data
obtained from the Exoplanet Orbit Database at http://exoplanets.org}.
The mean uncertainty on the published eccentricities is 0.05, though we
note that the uncertainties arising from least-squares fit can be
underestimated \citep{otoole09a}.  In addition, we excluded any
transiting planets (e.g.~HD\,17156b and HD\,80606b), because for these
cases, the transit should be simultaneously fit with the radial
velocities, a task which is beyond the scope of this paper.  After
applying these selection criteria, 82 stars remained.  A summary of the
data used here is given in Table~\ref{datasummary}.  All previously
unpublished AAT data used in this work are now given in the Appendix,
Tables~\ref{2039vels}--\ref{216437vels}.

To facilitate the comparison of the two-planet models with the single
eccentric-planet model, we first re-fit all available radial-velocity
data with a single planet (with no restrictions on $e$).  For those
stars with data from multiple sources, this approach ensures consistent
treatment by using the same fitting procedure for all stars.  These
results are referred to as Method 1 as given in Table~\ref{allfits}.
Then, for each system, we fit all available data using a genetic
algorithm, which has proven useful in previous work where the system
parameters are extremely uncertain or data are sparse \citep{NNSer,
HUAqr, 47205paper, tinney11, cochran07}.  The genetic algorithm used
here has the advantage that the range of allowed parameter space can be
restricted: in this work, we wish to fit the data with two
\textit{low-eccentricity} planets.  We thus direct the fitting process
to model two Keplerian orbits with $e<0.2$.  We note in passing that
applying this procedure to the data published in \citet{otoole07} yields
a 2-planet system essentially identical to that presented in
\citet{142paper}.  The best-fit set of parameters\footnote{While genetic
algorithms are commonly lauded as ways to find a \textit{global} best
fit, we note that their effectiveness depends on the choice of input
parameters such as mutation rates and the degree to which $\chi^2$ is
allowed to increase between generations. Hence, the ``global'' solutions
found by our approach cannot necessarily be guaranteed to be the
absolute best-fit for the the very complicated, multi-modal parameter
space of the second planet.} resulting from 10,000 iterations of the
genetic algorithm is then used as initial input for the
\textit{Systemic} Console \citep{mes09}.  We then use \textit{Systemic}
to perform Keplerian model fits to the data, again requiring $e<0.2$ for
both planets.  These results are referred to as Method 2 as given in
Table~\ref{allfits}.

The 12 systems for which a 2-planet fit gave a physically meaningful
result (i.e.~no crossing orbits) with $\chi^2$ similar to the one-planet
fit were subjected to a more detailed fitting process.  We used the
Runge-Kutta algorithm within \textit{Systemic} to perform a dynamical
fit which accounted for the gravitational interactions between the
modelled planets.  The model systems were then integrated for $10^5$
years as a basic stability check.  The results of this analysis, those
nine systems which remained stable for $10^5$ yr, are given in
Table~\ref{promising}.  For those systems which proved stable in the
initial check, we produced detailed dynamical maps of a broad range of
parameter space about the best-fit orbits.  For this final step in the
dynamical feasibility testing, we turned to the \textit{Mercury} N-body
integrator \citep{chambers99}.  Following our previous work in dynamical
mapping of extrasolar planetary systems \citep{24Sex, NNSer, HWVir,
155358, 204313, 142paper, HUAqr}, we sampled a 3$\sigma$ region of
four-dimensional parameter space: semimajor axis ($a$), eccentricity
($e$), mean anomaly ($M$), and argument of periastron ($\omega$).  Due
to computational limitations and the large uncertainties involved for
these speculative 2-planet systems, we chose coarse grids: 21x21x5x5 in
($a, e, M, \omega$), respectively.  As in previous work, we held the
best-constrained planet (planet 1 as given in Table~\ref{promising})
fixed and altered the initial parameters of the second planet.  We ran
each simulation for $10^8$ yr, or until the system destabilised (via
ejection or collision).

\section{Results}

We give the results of all orbit fits, including the reduced $\chi^2$
and rms velocity scatter, in Table~\ref{allfits}.  Many of the attempted
2-planet fits resulted in two Keplerians at nearly identical periods.
These systems may have a physically plausible solution with a slightly
worse $\chi^2$, but testing such possibilities is beyond the scope of
this paper.  Hence, there may exist additional ``good'' 2-planet
solutions which have been missed by our approach.  The Keplerian fitting
methods used here incorporate no physics: they are simply seeking a
lowest-$\chi^2$ solution regardless of the physicality of the resulting
system parameters.  A 1:1 resonant configuration is dynamically
possible, as evidenced by the abundance of such ``Trojan'' objects in
our own Solar system \citep{levison97, jonti1, jonti2}.  However, the
radial-velocity signature of a 1:1 configuration is extremely difficult
to disentangle \citep{laughlin02,g12}.  Some extrasolar planetary
systems have been proposed to host dynamically stable planets in 1:1
configurations \citep{goz06, cress09, schwarz09, funk12}.  However,
owing to the difficulty of maintaining dynamically stable
configurations, we consider such cases to be beyond the scope of this
work.  In this section, we will focus on those systems where a 2-planet
fit produced a ``plausible'' result (Table~\ref{promising}) with a
$\chi^2$ and rms similar to, or better than, the single-eccentric-planet
fit as given in Table~\ref{allfits}.

\textit{HD 3651} -- The best 2-planet fit resulted in a 2:1 period
commensurability, with the second signal at $P=31.08$ days.  This is
quite close to the monthly observing window, so the fit results in large
uncertainties in phase ($\omega$ and mean anomaly).  The $\chi^2$ and
rms are similar to, but slightly higher than, the one-planet fit.
Detailed dynamical simulations (Figure~\ref{dynamics1}, panel a) clearly
show the 2:1 resonance as a vertical strip of stability throughout the
range of allowed eccentricities, and the best fit for a second planet
places it comfortably within an extremely stable region.

\textit{HD 7449} -- The two-planet model results in two giant planets,
with minimum masses of 1.2 and 0.4 \mbox{M$_{\rm Jup}$}\ (Table~\ref{promising}).  The
dynamical tests (Figure~\ref{dynamics1}, panel b) show that the best-fit
eccentricity for the second (innermost) planet is on the edge of a
stable region.  As the eccentricity of the second planet is
increased beyond $e=0.4$, the system stability is quickly degraded.

\textit{HD 52265} -- As with HD\,3651, this system also gave a 2:1
configuration, with the second planet at half the period of the known
planet, and a slightly worse $\chi^2$ and rms.  Figure~\ref{dynamics1}
(panel c) shows that the candidate planet is well within a broad stable
region.

\textit{HD 65216} -- A second planet can be fit, with $P=152.6\pm$0.5
days and a mass of 0.17$\pm$0.03 \mbox{M$_{\rm Jup}$}.  The initial dynamical check in
\textit{Systemic} showed this system to be stable, and further dynamical
mapping (Figure~\ref{dynamics1}, panel d) shows that the entire 3-sigma
$a-e$ parameter space is stable.

\textit{HD 85390} -- Here, fitting a second planet improves the $\chi^2$
(12.2 to 4.3) and rms (2.3\,\mbox{m\,s$^{-1}$} to 1.4\,\mbox{m\,s$^{-1}$}).  The eccentricities of
both planets are then consistent with zero.  The second planet would be
a Jupiter analog (e.g.~Wittenmyer et al.~2011a), with $P\sim$10 yr and a
mass of 0.20$\pm$0.02\mbox{M$_{\rm Jup}$}.  The detailed dynamical map
(Figure~\ref{dynamics2}, panel a) shows a broad swath of stability for
all orbits with $a\gtsimeq$2\,AU.

\textit{HD 89744} -- The possible second planet has a period of 85.2
days, very close to a 3:1 commensurability with the known 258-day
planet.  The dynamics of the system (Figure~\ref{dynamics2}, panel b)
show that all permissible orbits are stable, despite the relatively high
mass of the candidate planet (3.2 \mbox{M$_{\rm Jup}$}).

\textit{HD 92788} -- Figure~\ref{dynamics2} (panel c) shows this planet
candidate is almost certainly trapped in the 2:1 resonance with the
known 326-day planet.  This resonance would allow relatively large
eccentricities (and even crossing orbits) to be long-term stable.

\textit{HD 117618} -- Adding a second planet improves the $\chi^2$ (8.1
to 5.7) and rms (5.6\mbox{m\,s$^{-1}$}\ to 4.6\mbox{m\,s$^{-1}$}).  With a period of 318$\pm$2 days,
the second signal is far enough from one year to allay fears of
aliasing.  A dynamical map for this Saturn-mass candidate is shown in
Figure~\ref{dynamics2} (panel d) -- again, there is a vast region of
stability across the 3$\sigma$ parameter space.

\textit{GJ 649} -- As the two proposed planets are very widely
separated, and of low mass, we did not perform the dynamical testing for
this system.  The planet candidates are separated by 21.8 mutual Hill
radii, which is certainly sufficient dynamical room for any interaction
to be negligible \citep{chambers96}.

\textit{HD 192310 (=GJ 785)} -- While this system does not appear in
Table~\ref{promising}, we note that the 2-planet fit shown in
Table~\ref{allfits} suggests a Neptune-mass second planet with period
629$\pm$64 days.  This is broadly consistent (within $\sim\,1.5\sigma$)
with the claim of a 0.07\mbox{M$_{\rm Jup}$}\ planet by \citet{pepe11}, with
$P=526\pm$9 days.

\citet{ang10} have also approached this problem of disentangling single
eccentric planets from near-circular two-planet systems, specifically
considering the case of the 2:1 resonance.  That work showed that about
35\% of known single-planet systems were indistiguishable from 2:1
resonant solutions.  By comparison, 9/82=11\% of the systems examined in
the present work resulted in a 2:1 configuration (Table~2).  We also
find that a further 24/82=29\% of systems tested resulted in co-orbital
(1:1) configurations.  \citet{hollis12} also performed an extensive and
self-consistent Bayesian re-analysis of available radial-velocity data
for 94 exoplanet systems.  As part of their analysis, they attempted
two-planet fits to the known single-planet systems (Hollis et al.~2012,
Table~A2).  Unlike this work, their two-planet fits were not restricted
to low eccentricities.  For example, \citet{hollis12} find a two-planet
solution for HD\,3651 where $(P_{1},e_{1})=(62.25,0.60)$ and
$(P_{2},e_{2})=(295,0.32)$.  By contrast, our result in
Table~\ref{promising} gives a 2:1 configuration with $e\ltsimeq\,0.06$
for both planets.

Ultimately, there is no substitute for sampling density: given infinite
observational resources, one would ideally observe every system as often
as possible \citep{monster}.  Such techniques have been successful in
identifying low-mass planets \citep{16417paper, 61vir, alphacen} and in
clarifying the true orbital period and mass of eccentric planets where
large velocity excursions occur on short timescales \citep{37605,
endl06}.  It is, of course, more practical to optimally plan one's
observations to confirm suspected planet candidates \citep{scheduling}.
All of the systems discussed above withstood detailed dynamical
scrutiny, as shown in Figures~\ref{dynamics1} and \ref{dynamics2}.  For
those potential 2-planet systems, we now ask how one may observationally
discern between the one- and two-planet models.  At present, the two
models are typically close to each other in goodness-of-fit; an example
is shown in Figure~\ref{comparenow} for HD\,65216.  In the future,
however, the two possible models may diverge.  In
Figures~\ref{orbits1}--\ref{orbits5}, we overplot the two models for
each system, to give a qualitative estimate of the optimal times to
observe them.  When they diverge sufficiently, well-timed observations
of reasonable precision ($\sim$2-3 \mbox{m\,s$^{-1}$}) can distinguish between the 1-
and 2-planet models.  We note, however, that the 2-planet fits often
have large phase gaps (and hence relatively large uncertainties in $P$
and $T_0$) which can shift the model curves shown here.  For this
reason, we leave more detailed quantitative simulations of detectability
(e.g. Wittenmyer et al. 2013) for future work.

\section{Summary and Conclusions}

We have examined 82 known moderately eccentric single-planet systems,
reanalysing the available radial-velocity data to test the hypothesis
that some may actually be low-eccentricity two-planet systems.  We have
identified nine particularly promising candidate systems, and performed
detailed dynamical stability simulations of the candidate planets.  All
of the systems proved to be dynamically stable on timescales of at least
$10^8$ yr.  We have also given model orbits to show qualitatively when
the one- and two-planet solutions diverge enough to be distinguishable
by future well-timed radial-velocity observations.  Our results suggest
that at least 11\% of apparently single-planet systems may in fact host
two low-eccentricity planets, a figure likely to rise as more
observations are obtained.  The difference between the one- and
two-planet solutions should typically be observable within the next 3
years, adding weight to the case for continued observations of these
systems.

\acknowledgements

We gratefully acknowledge the UK and Australian government support of
the Anglo-Australian Telescope through their PPARC, STFC and DIISR
funding; STFC grant PP/C000552/1, ARC Grant DP0774000 and travel support
from the Australian Astronomical Observatory.  The work was supported by
iVEC through the use of advanced computing resources located at the
Murdoch University, in Western Australia.  RW is grateful to Nanjing
University for the support of his visit to Nanjing.  S.~Wang's visit to
UNSW was funded by Faculty Research Grant PS27262.

This research has made use of NASA's Astrophysics Data System (ADS), and
the SIMBAD database, operated at CDS, Strasbourg, France.  This research
has made use of the Exoplanet Orbit Database and the Exoplanet Data
Explorer at exoplanets.org \citep{wright11}.


\section{Introduction}
A class of interaction potentials in quantum mechanics characterized by the
Hamiltonian,
\begin{equation}
\mathrm{H}=\mathrm{p}^2+\mathrm{r}^2+\lambda |\mathrm{r}|^{-\alpha}
\equiv \mathrm{H_0}+\lambda |\mathrm{r}|^{-\alpha}, \ \
\mathrm{r} \in [0, \infty]
\end{equation}
where $p=-i\ \partial /\partial \mathrm{r}$,
have found widespread applications in many areas of atomic, molecular, nuclear physics
and are often referred as the spiked harmonic oscillators (SHO). Here $\mathrm{H_0}$
formally denotes
the simple harmonic oscillator Hamiltonian; the coupling parameter $\lambda$
determines the strength of the perturbative potential and the positive constant
$\alpha$ represents the type of singularity at the origin. The higher the value of
$\lambda$, the higher the singularity. There has been an upsurge of interest [1-20]
to calculate the SHO eigenvalues over
the past three decades and it continues to grow. An interesting feature of
this potential is that once the perturbation $\lambda |\mathrm{r}|^{-\alpha}$
is switched on, complete turn-off is impossible; vestigial effects of the interaction
persists leading to the so called ``Klauder phenomenon" [1,2]. From a purely mathematical
viewpoint, on the other hand, this poses considerable challenges to some of the
well-established and widely used mathematical theories. For example, the commonly used
Rayleigh-Schr\"odinger perturbation series diverges [3] according to the relation
$\mathrm{n} \geq \frac{1}{\alpha-2}$, where $\mathrm{n}$ is the order of the
perturbation term. Consequently, a singular
perturbation theory was to be specially devised to treat these potentials. These
potentials also exhibit the phenomenon of {\it supersingularity} [1] in the region
of $\alpha \geq 5/2,$ i.e., every matrix element of the potential is infinite.
The numerical solutions of the pertinent Schr\"odinger equations are notoriously
difficult as well; especially those involving the finite-difference (FD) schemes
and often require special care.

Several analytical (both variational and perturbative) methodologies [3-20] are
available for the exact and approximate calculation of these systems.
For example, the modified (nonpower) perturbation series [3] to finite order
for the ground-state eigenenergies valid for small values of $\lambda$ and arbitrary
values of $\alpha$, a large coupling perturbative expansion [5] for the approximate
estimates of the same for large positive values of $\lambda$, the weak coupling expansion
expressions of the nonsingular ($\alpha < 5/2$) SHO through
the resummation technique for $\alpha = 1/2, 1, 3/2\ $ [6] and for $\alpha = 2$ [10],
the exact and approximate (variational) solutions [8] for some particular values of the
parameters in the interaction potential, a modified WKB treatment [11], etc. Besides,
the upper and lower bounds of ground and excited states [12-16] of the SHO as well as
the generalized SHO, the analytical
pseudoperturbation shifted-$\ell$ expansion technique (PSLET) [19,20] have also been
developed. The extensions to N-dimensions are reported lately [14,20]. On the numerical
side, the FD methods [21,22] through Richardson extrapolation, integration
of the Schr\"odinger equation [23] using a Lanczos grid method for the cases of
$\alpha = 4,6$ for small values of the ($\lambda \leq 0.01$), the analytic continuation
method [24] for ground and excited states have been reported.

Despite all these attempts, a general prescription which can accurately and reliably
calculate the bound states of these potentials in a uniform and simple way for a general
set of potential parameters with the promise of furnishing ground and excited states
with equal ease, would be highly desirable. This is because
physically meaningful and good accuracy results are obtainable only by some of these
methods. Additionally some of these methods can give
satisfactory results for a certain type of parameters while perform rather poorly in
other cases. Much attention has been paid to the ground states; excited states are
reported less frequently and definitively, presumably because of the greater challenges
compared to the ground states. Moreover, much work has been devoted to the eigenvalues;
only few results are available for the {\it eigenfunctions}[15]. Also some
of these methodologies are often fraught with rather tedious and cumbersome
mathematical complexities. This work proposes a simple methodology to study these
systems by using the GPS scheme which has shown considerable promise for a variety
of atomic and molecular processes including both static and dynamic situations in recent
years (see, for example [25-29] and the references therein). This
formalism helps alleviate some of the well-known discomfitures of the
FD schemes widely used and discussed in the literature [5,21,22,30],
e.g., the necessity of significantly larger spatial grid points to deal with
the singularity at the origin. The GPS method essentially works in a {\it nonuniform}
and {\it optimal} spatial grid; thus a much smaller number of points suffices to
achieve good accuracy. However its applicability has been so far restricted to the
cases of Coulomb singularities; no attempts have been made to deal with the other
singularities characterizing many physical systems. The objective of this Letter is
two-fold: (a) to extend the regions of applicability and judge the performance of it
on the SHOs, (b) to calculate accurately the bound-state spectra of these systems.
Comparison with the literature data has been made wherever possible.

The organization of the article is as follows. Section II presents an overview of the
basic formalism. Section III makes a discussion of the results while a few
concluding remarks are made in section IV.

\section{The GPS formalism for the SHO}\label{sec:method}
In this section, we present the GPS formalism used to solve the radial
eigenvalue problem with the SHO potentials. A detailed account of the GPS method
can be found in the references [25-29].

The time-independent radial Schr\"odinger equation to be solved can be written in the
usual way (atomic units employed unless otherwise mentioned),
\begin{equation}
\left[-\frac{1}{2} \ \frac{\mathrm{d^2}}{\mathrm{dr^2}} + \frac{\ell (\ell+1)}
{2\mathrm{r^2}} + \mathrm{v(r)} \right]
\psi_{n,\ell}(\mathrm{r})=E_{n,\ell}\ \psi_{n,\ell}(\mathrm{r})
\end{equation}
where v(r) is the SHO potential given by,
\begin{equation}
\mathrm{v(r)}=[\mathrm{r}^2+\lambda/\mathrm{r}^{\alpha}]/2.
\end{equation}
The 1/2 factor is introduced here only for easy comparison
with the literature and $\ell$ signifies the usual angular momentum quantum number. The
GPS formalism facilitates the use of a denser mesh at small r regions and a relatively
coarser mesh at the large r regions while preserving the similar accuracy at both the
regions.

The key step in this formalism is to approximate a function $f(x)$ defined in the
interval $x \in [-1,1]$ by an N-th order polynomial $f_N(x)$ {\it exactly} at the
discrete collocation points $x_j$ as in the following,
\begin{equation}
f(x) \cong f_N(x) = \sum_{j=0}^{N} f(x_j)\ g_j(x),
\end{equation}
\begin{equation}
f_N(x_j) = f(x_j).
\end{equation}
Within the Legendre pseudospectral method that the current work uses, $x_0=-1$,
$x_N=1$, and $x_j (j=1,\ldots,N-1)$ can be determined from the roots of the first
derivative of the Legendre polynomial $P_N(x)$ with respect to $x$, i.e.,
\begin{equation}
P'_N(x_j) = 0.
\end{equation}
The $g_j(x)$s in Eq.~(4) are the cardinal functions expressed as,
\begin{equation}
g_j(x) = -\frac{1}{N(N+1)P_N(x_j)}\ \  \frac{(1-x^2)\ P'_N(x)}{x-x_j},
\end{equation}
satisfying the relation $g_j(x_{j'}) = \delta_{j'j}$.
At this stage, we use a transformation $r=r(x)$ to map the semi-infinite domain
$r \in [0, \infty]$ onto the finite domain $x \in [-1,1]$. One can make use of the
following algebraic nonlinear mapping,
\begin{equation}
r=r(x)=L\ \ \frac{1+x}{1-x+\alpha},
\end{equation}
where L and $\alpha=2L/r_{max}$ are the mapping parameters. Finally introduction of the
following relation,
\begin{equation}
\psi(r(x))=\sqrt{r'(x)} f(x)
\end{equation}
in conjunction with a symmetrization procedure gives the following transformed Hamiltonian
\begin{equation}
\hat{H}(x)= -\frac{1}{2} \ \frac{1}{r'(x)}\ \frac{d^2}{dx^2} \ \frac{1}{r'(x)}
+ \mathrm{v}(r(x))+v_m(x),
\end{equation}
The advantage of this is that one ends up with a \emph {symmetric} matrix eigenvalue
problem which can be solved readily and efficiently to give accurate eigenvalues
and eigenfunctions by using standard routines. It may be noted that $v_m(x)=0$ for the
above transformation leading to the following set of discretized coupled equations,
\begin{widetext}
\begin{equation}
\sum_{j=0}^N \left[ -\frac{1}{2} D^{(2)}_{j'j} + \delta_{j'j} \ v(r(x_j))
+\delta_{j'j}\ v_m(r(x_j))\right] A_j = EA_{j'},\ \ \ \ j=1,\ldots,N-1,
\end{equation}
\end{widetext}
where
\begin{equation}
A_j  = \left[ r'(x_j)\right]^{1/2} \psi(r(x_j))\ \left[ P_N(x_j)\right]^{-1}.
\end{equation}
and the symmetrized second derivatives $D^{(2)}_{j'j}$ of the cardinal functions are
given in [26]. Thorough checks are made on the variation of the energies with respect
to the mapping parameters for large ranges of the interaction parameters in the
potential available in the literature. After a series
of such calculations, a choice has been made at the point where the results changed
negligibly with any variation. In this way, a consistent and uniform set of parameters
($r_{max}=200,$ $\alpha=25$ and $N=300$) has been used.

\section{Results and discussion}
\subsection{The charged harmonic oscillator, $\alpha=1$}
Before considering the general case of relatively stronger spikes, {\it viz.,}
$\alpha \neq 1$, it is worthwhile to study the simpler special case of $\alpha=1$.
This does not exhibit supersingularity and the Hamiltonian takes the
simplified confined Coulomb potential type form. It has been pointed out [31] that
this possesses an infinite set of {\em elementary} solutions.
Table~\ref{tab:table1} displays such elementary solutions calculated by the present
method along with the exact analytical results. Note that $E=3/2$ is a trivial
solution corresponding to $\lambda=0$, i.e., the unperturbed Hamiltonian. The
other $\lambda$s in this table are taken from the solutions of the polynomial
equation [31]. It may be noted that all the calculated results in this table and
throughout the article are {\em truncated} and therefore all the digits in the
reported numbers should be considered as correct. It is seen that
for all values of $\lambda$, our results match excellently up to a 12 digit accuracy
with the exact values.

\begingroup
\squeezetable
\begin{table}
\caption {\label{tab:table1}Some elementary solutions (in a.u.) of the SHO
with $\alpha=1$ for several values of $\lambda$ corresponding to the ground
state.}
\begin{ruledtabular}
\begin{tabular}{ccc}
$\lambda$ & \multicolumn{2} {c} {Energy} \\
\cline{2-3}
     & This work & Exact\footnote{Ref. [31]. These results have
been halved to take care of a 2 factor.}  \\  \hline
 0                          & 1.49999999999  &  1.5   \\
 2                          & 2.49999999999  &  2.5   \\
 $\sqrt{20} $               & 3.50000000000  &  3.5   \\
 $(30+6\sqrt{17})^{1/2}$    & 4.49999999999  &  4.5   \\
 $(70+6\sqrt{57})^{1/2}$    & 5.49999999999  &  5.5   \\
 14.450001026966            & 6.49999999999  &  6.5   \\
 18.503131410003            & 7.49999999999  &  7.5   \\
\end{tabular}
\end{ruledtabular}
\end{table}
\endgroup

\begingroup
\squeezetable
\begin{table}
\caption {\label{tab:table2}Calculated ground-state energies E (in a.u.) of the
SHO with $\alpha=1$ for several values of $\lambda.$}
\begin{ruledtabular}
\begin{tabular}{cccccl}
$\lambda$ & \multicolumn{2}{c}{Energy}& $\lambda$ & \multicolumn{2}{c}{Energy} \\
\cline{2-3} \cline{5-6}
         & This work & Literature\footnotemark[1]  &      & This work &
Literature\footnotemark[1]     \\  \hline
$-0.001$ & 1.49943577146 & 1.49943577146  & 0.001 & 1.50056415064 & 1.5005641506   \\
$-0.005$ & 1.49717807794 &                & 0.005 & 1.50281997477 &                \\
$-0.01$  & 1.49435420563 & 1.49435420565  & 0.01  & 1.50563800525 & 1.50563800525  \\
$-0.05$  & 1.47169265799 &                & 0.05  & 1.52811261097 &                \\
$-0.1$   & 1.44318757957 & 1.4431875796   & 0.1   & 1.55603345324 & 1.55603345325  \\
$-0.5$   & 1.20765362342 &                & 0.5   & 1.77283783394 &                \\
$-1$     & 0.892602739638 & 0.89260273965 & 1     & 2.02893850398 & 2.0289385040   \\
$-5$  & $-2.90807895034$ & $-2.90807895035$ & 5   & 3.69201586294 & 3.69201586295  \\
$-10$ & $-12.4404995301$ & $-12.44049953015$ & 10 & 5.28874176968 & 5.288741697    \\
$-50$ & $-312.497600033$ &                 & 50   & 13.7025706824 &                \\
$-100$ & $-1249.99940000$ &                & 100  & 21.2314590573 &                \\
\end{tabular}
\end{ruledtabular}
\footnotetext[1]{Ref. [31]. The quoted results are halved to
take care of a 2 factor.}
\end{table}
\endgroup

Next we report in Table~\ref{tab:table2} the ground-state energies for a large range
of (+)ve and (-)ve $\lambda$s (left and right sides of the table respectively)
along with the available literature data. One can envisage three distinct regions
in this case depending on the values of $\lambda$, {\em viz.,} (a) the Coulomb region,
corresponding to large (-)ve $\lambda,$ (b) the strong-coupling region having large
(+)ve $\lambda,$ and (c) the weak-coupling region having small (both (+)ve and
(-)ve) $\lambda$. The perturbation expressions corresponding to regions (a) and
(c) are obtained through an amalgamation of the hypervirial
and Hellmann-Feynman theorem [31]. For some of the (-)ve and (+)ve $\lambda$s
ground states are examined by the Renormalization as well as the direct
numerical integration methods [31]. Also for $\lambda \leq -1$ and $\lambda \geq 1$,
the Coulomb series and strong coupling series solutions are available [31].
Good agreement is observed for $\lambda=-10$ and $\lambda=10$ involving these methods;
for other $\lambda$s, they vary significantly from each other. Here, the numerical
results are quoted for comparison. No results were available for $\lambda = \pm 0.005,
\pm 0.05, \pm 0.5, \pm 50, \pm 100)$. It is seen that the current results
are in excellent agreement with theirs. At this point mention may be made of one of
the uncomfortable features in some of the available methodologies, {\em viz.,} the
presence of the unphysical roots, e.g., in the Riccati-Pad\'e method for the small
$\lambda$s of these potentials [31]. However, no such solutions have been found
in the present calculations. In some instances, very slight differences
are observed in our results from the literature data. Furthermore, in
table~\ref{tab:table3}, we present the calculated first three states
corresponding to $\ell=0,1,2,3$ for these systems. The ground states are repeated
for the sake of completeness. Again a wide range of both positive and negative
$\lambda$ values are chosen. No results are available for these states to our
knowledge and we hope that they could be useful in future calculations.

\begingroup
\squeezetable
\begin{table}
\caption {\label{tab:table3} Excited state energies (in a.u.) of the charged harmonic
oscillator for several positive and negative values of $\lambda$. First three states
are presented corresponding to $\ell=0,1,2,3.$}
\begin{ruledtabular}
\begin{tabular}{cllll}
 $\lambda$  & $\ell=0$  &   $\ell=1$      &  $\ell=2$     &   $\ell=3$     \\ \hline
 $-0.001$ & 1.49943577146   & 2.49962386468    & 3.49969909505    &  4.49974208263  \\
         & 3.49952982655   & 4.49966148110    & 5.49972058921    &  6.49975641177   \\
         & 5.49958155072   & 6.49968700564    & 7.49973670957    &  8.49976780986   \\
$-0.1$   & 1.44318757957   & 2.46229789284    & 3.46987169094    &  4.47418745667   \\
         & 3.45282982176   & 4.46609746678    & 5.47203352713    &  6.47562590482   \\
         & 5.45807015200   & 6.46866692538    & 7.47365242360    &  8.47676916627   \\
$-10  $  &$-12.4404995301$ &$-2.62119802134$  & 0.004574720607   & 1.67473695981    \\
         &$-2.41723883317$ & 0.551233459914   & 2.41797686878    & 3.89431134350    \\
         & 0.869699218970  & 3.039428778956   & 4.66996941607    & 6.05143674320    \\
$-100$   &$-1249.99940000$ &$-312.494000168$  &$-138.863694966$  &$-78.0530829912$  \\
         &$-312.491600236$ &$-138.852898182$  &$-78.0243305403$  &$-49.7760242855$  \\
         &$-138.847499606$ &$-78.0051588012$  &$-49.7312779249$  &$-34.2071382286$  \\
0.001    & 1.50056415063   & 2.50037611746    & 3.50030089728    & 4.50025791310    \\
         & 3.50047014252   & 4.50033850860    & 5.50027940556    & 6.50024358500    \\
         & 5.50041843194   & 6.50031298742    & 7.50026328651    & 8.50023218755    \\
  0.1    & 1.55603345324   & 2.53752389333    & 3.53005208746    & 4.52577056764    \\
         & 3.54686142702   & 4.53380037959    & 5.52791522734    & 6.52434325783    \\
         & 5.54175768729   & 6.53126513099    & 7.52631014389    & 8.52320695048    \\
  10     & 5.28874176968   & 5.63241238009    & 6.19962012502    & 6.89697621559    \\
         & 7.07543947857   & 7.46149523127    & 8.06916090729    & 8.79599713087    \\
         & 8.89811648443   & 9.32030795656    & 9.96020685568    & 10.7102538225    \\
  100    & 21.2314590573   & 21.3064355531    & 21.4546955286    & 21.6730703270    \\
         & 22.9756496882   & 23.0536647120    & 23.2077138583    & 23.4341264159    \\
         & 24.7309414007   & 24.8119454465    & 24.9716696244    & 25.2059214990    \\
\end{tabular}
\end{ruledtabular}
\end{table}
\endgroup

\begingroup
\squeezetable
\begin{table}
\caption {\label{tab:table4}Calculated ground-state energies E (in a.u.) of the SHO
with $\alpha=4$ and 6 for several values of $\lambda.$ The
literature results are divided by a 2 factor.}
\begin{ruledtabular}
\begin{tabular}{cllll}
$\lambda$ & \multicolumn{2}{c}{Energy ($\alpha=4$)} &
\multicolumn{2}{c}{Energy ($\alpha=6$)} \\
\cline{2-3} \cline{4-5}
      & This work   & Literature   & This work    & Literature \\   \hline
0.001 & 1.53438158545 & 1.53438158545\footnotemark[1], 1.534385\footnotemark[2] &
1.63992791296 & 1.63992791296\footnotemark[1]  \\
0.005 & 1.57417615416 & 1.574176155\footnotemark[3],1.574175\footnotemark[4],    &
1.71144209213 & 1.71144208\footnotemark[3],1.71144\footnotemark[4],    \\
&     & 1.574195\footnotemark[5] &  & 1.71151\footnotemark[5] \\
0.01  & 1.60253374753 & 1.60253374753\footnotemark[1], 1.60254\footnotemark[2],  &
1.75272613799 & 1.75272613799\footnotemark[1],1.752726195\footnotemark[3], \\
      &   & 1.602533745\footnotemark[3],1.602535\footnotemark[4],   &    &
1.752725\footnotemark[4],1.75287\footnotemark[5],       \\
  &       & 1.602635\footnotemark[5], 1.602535\footnotemark[6]    &
  &  1.7527265\footnotemark[6]       \\
0.05  & 1.71258069752 &                        & 1.88277010302 &          \\
0.1 & 1.78777599560 & 1.78777599560\footnotemark[1],
     1.787785\footnotemark[2] & 1.95783261264 &         \\
   &  & 1.787775\footnotemark[6]   &    &               \\
0.5   & 2.06529243634 &                        & 2.19395453013 &          \\
1     & 2.24708899168 & 2.24708899168\footnotemark[1],
     2.24709\footnotemark[2]$^,$\footnotemark[6]  & 2.32996998478 &
2.32996998478\footnotemark[1],2.329970\footnotemark[6]  \\
5  & 2.89222177088 & 2.89222\footnotemark[6]
   & 2.75657950709 & 2.7565795\footnotemark[6]               \\
10    & 3.30331125601 & 3.30331125601\footnotemark[1],
3.30331\footnotemark[2]$^,$\footnotemark[6]  & 3.00160451444 &
3.00160451444\footnotemark[1],3.0016045\footnotemark[6], \\
   &  & 3.3033112560\footnotemark[7]   &    &  3.00160451\footnotemark[7]  \\
50    & 4.73277787167 &                        & 3.76776072255 &         \\
100   & 5.63254021587 & 5.63254021587\footnotemark[1], 5.63254\footnotemark[2],  &
4.20667914031 & 4.2066791403\footnotemark[7] \\
    &   & 5.6325402\footnotemark[7]  &      &      \\
500   & 8.73793385806 &                        & 5.57607711626 &              \\
1000  & 10.6847312660 & 10.6847312660\footnotemark[1], 10.68473\footnotemark[2], &
6.35930853290 &    \\
   &   &  10.684731265\footnotemark[7]   &    &    \\
\end{tabular}
\end{ruledtabular}
\footnotetext[1] {Ref. [24].}
\footnotetext[2] {Ref. [31].}
\footnotetext[3] {Ref. [23].}
\footnotetext[4] {Ref. [21].}
\footnotetext[5] {Ref. [1].}
\footnotetext[6] {Ref. [8].}
\footnotetext[7] {Ref. [33].}
\end{table}
\endgroup

\subsection{$\alpha \neq 1$}
Now results are presented for $\alpha \neq 1$. Here we focus on the $\alpha$ values
4 and 6; however,  the present scheme has been thoroughly checked to reproduce the
results of similar accuracy and reliability for other values of $\alpha$
available in the literature. In table~\ref{tab:table4}, ground state energies are
tabulated for these two cases ($\alpha=4$ in the left and $\alpha=6$ in the right),
for small and large $\lambda$s. Two new $\lambda$ values
are introduced here (500 and 1000) in addition to those employed in
table~\ref{tab:table2}. Both of these $\alpha$ values can lead to supersingularity and
have been investigated by many workers.  The present results are seen to be in good
agreement with the accurate analytic continuation results [24]. These results are
available for $\lambda= 0.001, 0.01, 1, 10$ for both $\alpha=4,6$, while
$\lambda=0.1, 100, 1000$ for $\alpha=4$ only. Various other results are also
available for the smaller $\lambda$s (0.005, 0.01) [23,21,1,33]
and our results show good agreement with these. Direct integration results
[8] for $\lambda=0.01, 0.1, 1, 5, 10$ are also presented for comparison. It may be
noted that the current results surpass in accuracy to all others except [24].

\begingroup
\squeezetable
\begin{table}
\caption {\label{tab:table5}Calculated $\ell \ne 0$ state energies (in a.u.) of the
SHO with $\alpha=4$ (top) and 6 (bottom) for several values of $\lambda$s. The
literature results have been divided by a 2 factor.}
\begin{ruledtabular}
\begin{tabular}{cccccc}
$\ell$ & $\lambda=0.001$ & $\lambda=0.01$ & $\lambda=0.1$ & $\lambda=1$ &
$\lambda=10$ \\ \hline
 3 & 4.50005713956 & 4.50057109970 & 4.50568201308 & 4.55432930375 & 4.91961566042  \\
   & 4.50005713956\footnotemark[1] & 4.50057109970\footnotemark[1] &
     4.50568201309\footnotemark[1] & 4.55432930376\footnotemark[1] &                \\
 4 & 5.50003174537 & 5.50031739444 & 5.50316804961 & 5.53112085969 & 5.77200022575  \\
   & 5.50003174537\footnotemark[1] & 5.50031739444\footnotemark[1] &
     5.50316804961\footnotemark[1] & 5.53112085969\footnotemark[1] &                \\
 5 & 6.50002020182 & 6.50020200030 & 6.50201821626 & 6.52000759152 & 6.68566506197  \\
   & 6.50002020182\footnotemark[1] & 6.50020200030\footnotemark[1] &
     6.50201821626\footnotemark[1] & 6.52000759153\footnotemark[1] &                \\
10 & 11.5000050125 & 11.5000501247 & 11.5005011983 & 11.5050070693 & 11.5495902896  \\
20 & 21.5000012507 & 21.5000125077 & 21.5001250765 & 21.5012506189 & 21.5124915772  \\
30 & 31.5000005556 & 31.5000055570 & 31.5000555707 & 31.5005556887 & 31.5055549875  \\
40 & 41.5000003125 & 41.5000031254 & 41.5000312547 & 41.5003125438 & 41.5031249904  \\
50 & 51.5000001999 & 51.5000020001 & 51.5000200019 & 51.5002000183 & 51.5020000374  \\
   &               &               &               & 51.5002000183\footnotemark[1] & \\
\hline
5  & 6.50000577192 & 6.50005771227 & 6.50057643602 & 6.50570148892 & 6.55258902874  \\
10 & 11.5000005897 & 11.5000058970 & 11.5000589687 & 11.5005894939 & 11.5058757159  \\
20 & 21.5000000675 & 21.5000006760 & 21.5000067609 & 21.5000676086 & 21.5006759799  \\
30 & 31.5000000194 & 31.5000001949 & 31.5000019498 & 31.5000194985 & 31.5001949796  \\
40 & 46.5000000080 & 41.5000000811 & 41.5000008117 & 41.5000081181 & 41.5000811806  \\
50 & 51.5000000040 & 51.5000000411 & 51.5000004123 & 51.5000041240 & 51.5000412410  \\
\end{tabular}
\end{ruledtabular}
\footnotetext[1] {Ref. [16].}
\end{table}
\endgroup

Now table~\ref{tab:table5} gives the results for low- and high-$\ell$ states for a wide
range of $\lambda$s (0.001, 0.01, 0.1,1,10) for both $\alpha=4,6$. Upper and
lower bounds as well as the numerical eigenenergies for $\ell=3,4,5$ have been studied
recently by [16] for $\alpha=4$ for first four $\lambda$s. Our results match almost
completely with theirs except very slight discrepancies in three instances at the
last digit (our results are lower by $10^{-11}$). Also, the eigenvalues of
$\ell=5,10,20 \cdots ,50$ for all the mentioned values of $\alpha$
and $\lambda$ are given as a test of this method for the very high excited states.
The present result is in complete agreement with the lone available result of
$\ell=50$ (for $\alpha=4, \lambda=1$). Next table~\ref{tab:table6} gives results
for the first 10 eigenvalues of the SHO with the parameters $\alpha=6, \lambda=10$.
We have considered $\ell=0,1,2,3,4$ and no results could be found for these states.

\begingroup
\squeezetable
\begin{table}
\caption{\label{tab:table6}The first 10 eigenvalues (in a.u.) for $\ell=0,1,2,3,4,$
of the SHO. The parameters are: $\alpha=6$ and $\lambda=10.$}
\begin{ruledtabular}
\begin{tabular}{ccccc}
$\ell=0$      &  $\ell=1$      & $\ell=2$       & $\ell=3$       & $\ell=4$  \\  \hline
 3.00160451444 & 3.32389487858 & 3.91806927392 & 4.70345973112 & 5.60034573904  \\
 5.38666914834 & 5.64859182812 & 6.14941053052 & 6.84505383654 & 7.67818492569  \\
 7.66493489996 & 7.89510635043 & 8.34183646535 & 8.97697716408 & 9.75876666401  \\
 9.88940298242 & 10.0989808012 & 10.5089815278 & 11.0996720737 & 11.8395020670  \\
 12.0805100442 & 12.2752417097 & 12.6580500545 & 13.2142245371 & 13.9191144855  \\
14.2485304976 & 14.4318512146 & 14.7933959369 & 15.3217108425 & 15.9969892830  \\
16.3994450470 & 16.5736092280 & 16.9178854110 & 17.4230621225 & 18.0728546082  \\
18.5370785045 & 18.7036642068 & 19.0335218023 & 19.5190555671 & 20.1466196426  \\
20.6640433547 & 20.8242086533 & 21.1417690049 & 21.6103351695 & 22.2182898578  \\
22.7822131369 & 22.9368391875 & 23.2437342203 & 22.6974360875 & 24.2879213038  \\
\end{tabular}
\end{ruledtabular}
\end{table}
\endgroup

\begingroup
\squeezetable
\begin{table}
\caption {\label{tab:table7}Calculated expectation values (in a.u.) for the SHO
for some selected values of $\alpha\ $ and $\lambda\ $. The first three states
corresponding to $\ell=0,1,2\ $ are presented.}
\begin{ruledtabular}
\begin{tabular} {ccccc}
$\alpha$ &$\lambda$ &$\ell$ &$\langle r^{-1}\rangle$ & $\langle r \rangle $ \\ \hline
 1   & 10      &   0     & 0.579335567 & 1.88860444    \\
     &         &         & 0.572186022 & 2.20351385    \\
     &         &         & 0.562374825 & 2.49562513    \\
 4   & 10      &   0     & 0.546623313 & 1.93946889    \\
     &         &         & 0.483512472 & 2.38064959    \\
     &         &         & 0.443751364 & 2.74044976    \\
 6   & 10      &   0     & 0.558986259 & 1.89176957    \\
     &         &         & 0.477245223 & 2.38688068    \\
     &         &         & 0.431013450 & 2.77224304    \\
\end{tabular}
\end{ruledtabular}
\end{table}
\endgroup

As a test on the quality of the eigenfunctions, In table~\ref{tab:table7}, we present
some of the calculated expectation values $\langle r^{-1} \rangle$ and
$\langle r \rangle $ for $\alpha=1,4$ and 6. The parameter
$\lambda$ is kept fixed at 10 in all these cases and the first three states are
reported for $\ell=0,1,2$. No results could be found for any of these
values in the literature. Finally, figure 1 depicts the radial probability
distribution functions for the first three states of $\ell=0,1,2$ along
with the potential ($\alpha=6, \lambda=10$). As expected they show the
requisite number of nodes in these plots.

\begin{figure}
\begin{minipage}[b]{0.40\textwidth}
\centering
\includegraphics[scale=0.30]{fig1d.eps}
\end{minipage
\\
\begin{minipage}[b]{0.40\textwidth}
\centering
\includegraphics[scale=0.30]{fig1c.eps}
\end{minipage
\\
\begin{minipage}[b]{0.40\textwidth}
\centering
\includegraphics[scale=0.30]{fig1b.eps}
\end{minipage
\\
\begin{minipage}[b]{0.40\textwidth}
\centering
\includegraphics[scale=0.30]{fig1a.eps}
\end{minipage

\caption{The radial probability distribution function, $|rR_{n\ell}|^2$ (in a. u.)  for
the first three states corresponding to $\ell=0,1,2$ of the spiked oscillator potential
having the parameters, $\alpha=6,\ \lambda=10.\ $ (a) The potential, (b) the ground
state, (c) the first excited state and (d) the second excited state.}\label{fig:fig1}
\end{figure}

\section{Conclusion}
The GPS formalism is shown to deliver {\em accurate} and {\em reliable} results for
the eigenvalues, expectation values and the radial densities of the SHOs. The
simplicity and viability of the method is demonstrated by calculating the low and
high excited states of these potentials for weak and strong values of the interaction
parameter in the potential. Excellent agreement with the literature data is observed
in all cases. Some states are reported here for the first time. Finally the approach
may be as well equally successful and useful for other singularities (e. g., the
Hulth\'en, Yukawa, Hellman potentials etc.) in quantum mechanics. Work in this direction is
under progress.

\begin{acknowledgments}
I gratefully acknowledge the hospitality provided by the Department of Chemistry,
University of New Brunswick, Fredericton, Canada.
\end{acknowledgments}


\section{Introduction}

This paper addresses spontaneous symmetry breaking in non-relativistic
quantum systems of finite size. Strictly speaking, spontaneous
symmetry breaking can only occur in infinite systems. Then the ground
state exhibits a lesser degree of symmetry than the Hamiltonian
itself, i.e., the ground state is invariant under a symmetry group
${\cal H}$ that is a proper subgroup of the Hamiltonian's symmetry
group ${\cal G}$. The low-energy excitations are strongly constrained
by symmetry and given in terms of (weakly interacting) Nambu-Goldstone
modes. These can be calculated within an effective field theory (EFT)
that is solely based on the pattern of symmetry breaking. From a
technical point of view, the EFT is a nonlinear $\sigma$ model with
fields that parameterize the coset space ${\cal G}/{\cal
  H}$~\cite{weinberg1968,coleman1969,callan1969}. Examples for
spontaneous symmetry breaking are the breaking of spin-rotational
symmetry in a ferromagnet, the breaking of translational symmetry in a
crystal lattice, and the breaking of chiral symmetry in quantum
chromodynamics (QCD). In these examples, the Nambu-Goldstone modes are
magnons, phonons, and pions, respectively, and EFTs have been
developed for all these
cases~\cite{Gasser1984,Leutwyler1994,Leutwyler1996,roman1999,hofmann1999,Baer2004,kampfer2005},
see Refs.~\cite{weinbergbook,brauner2010} for reviews.

In finite systems, the ground state exhibits the full symmetry of the
Hamiltonian, and spontaneous symmetry breaking becomes evident in
symmetry-unrestricted mean-field
calculations~\cite{aberg1990,Nazarewicz1994,Yannouleas2007}. It is
then a major effort (and complication!) to restore the symmetry with
the help of projection techniques. The expressions ``obscured symmetry
breaking''~\cite{Koma1994} or ``emergent symmetry
breaking''~\cite{Yannouleas2007} (which we adopt here) have been
proposed for such systems. There are two distinct cases for which
emergent symmetry breaking plays a role. First, numerical simulations
of infinite physical systems are usually limited to a finite volume,
and it is then important to understand the finite-size corrections.
Some rigorous results are known in this
case~\cite{Horsch1988,Koma1994}, and finite-size corrections to
partition functions and thermodynamical observables have been worked
out within EFTs for simulations of QCD on finite
lattices~\cite{Leutwyler1987,Gasser1988}, and for spin
systems~\cite{Hasenfratz1993}. Genuinely finite systems constitute the
second and probably physically most interesting case. Prominent
examples are the emergence of superfluidity in trapped ultracold Bose
gases~\cite{matsumoto2002,Enomoto2006}, pairing in atomic nuclei (both
breaking a $U(1)$ phase symmetry), and non-spherical shapes of
molecules and atomic nuclei (both breaking $O(3)$ rotational symmetry
in the limit of infinite system size). Here, the techniques for
constructing EFTs for spontaneous symmetry breaking need to be
modified, and the interest is in spectra and transitions rather than
in thermodynamical observables.

The EFT for a finite system with emergent symmetry breaking is, of
course, related to the EFT for the corresponding infinite system with
spontaneous symmetry breaking. The symmetry must be realized
nonlinearly, and the Nambu-Goldstone fields parameterize the coset
space ${\cal G}/{\cal H}$~\cite{weinberg1968,coleman1969,callan1969}.
In the infinite system, the proper Nambu-Goldstone fields depend on
space and time and exhibit fluctuations of small amplitudes and long
wavelengths. A purely time-dependent (and spatially constant) mode is
forbidden because it would relate states of inequivalent Hilbert
spaces. In the finite system, however, this zero mode, i.e., the
spatially constant mode of the Nambu-Goldstone field, must be singled
out and treated separately. That mode undergoes large-amplitude
fluctuations and upon quantization restores the symmetry. In the
finite system, the small-amplitude fluctuations, i.e., the proper
Nambu-Goldstone modes with nontrivial temporal and spatial dependence,
must likewise be quantized. The theoretical implementation of this
program is not trivial and is demonstrated for two interesting and
important cases. We first consider as an example the emergent breaking
of a $U(1)$ phase symmetry in finite superfluids such as ultracold
bosonic atom gases or atomic nuclei. In this case, the proper
Nambu-Goldstone modes and the global phase rotations do not couple in
leading order, and both modes have the same energy scale. This
facilitates the description. Second, we consider the emergent breaking
of $SO(3)$ symmetry to its $SO(2)$ subgroup. This case describes the
low-energy physics of nonspherical objects with axial symmetry such as
linear molecules and deformed atomic nuclei. The case is more
complicated and interesting due to the interactions between global
rotations and proper Nambu-Goldstone modes, and the energy scale of
the rotational mode differs from the energy scale of the proper
Nambu-Goldstone modes. Our detailed presentation of these two problems
makes clear how to develop EFTs for systems with emergent symmetry
breaking in general.

In this paper we construct EFTs for finite systems with emergent
symmetry breaking, and we focus particularly on deformed atomic
nuclei. Such nuclei are traditionally described by generalized
collective
models~\cite{bohr_1952,bohrmottelson_1953,eisenberg,bmbook,rowe} or
the interacting boson model~\cite{arima1975,iachello}.  For more
microscopic approaches to collective motion, we refer the reader to
Refs.~\cite{aberg1990,frauendorf2001,niksic2011}.  For deformed
nuclei, the presented EFT generalizes the simpler construction of an
effective theory proposed recently~\cite{papenbrock2011}. Based on
symmetry principles alone, our model-independent approach re-derives
some of the well-known results for collective nuclear models.  We
expect that extensions of the EFT approach could be useful in
addressing well-known and long-standing limitations of the collective
models such as, e.g., the significant overprediction of
electromagnetic interband transitions, see
Refs.~\cite{garrett2001,rowe} for recent reviews of this problem.

Our procedure is patterned after the case of the infinite
ferromagnet~\cite{Leutwyler1994,roman1999, kampfer2005}.  We
generalize the expression for unitary transformations in the coset
space by including purely time-dependent variables. These account for
the dynamics of the finite system. The resulting generators define the
Nambu-Goldstone modes as well as the zero modes as classical fields.
From these we construct the building blocks of the effective
Lagrangian $L$ using arguments of invariance and energy scaling.
Quantization of the Hamiltonian obtained by a Legendre transformation
of $L$ then determines the spectrum.

\section{Emergent breaking of U(1) phase symmetry}
\label{emer}

Superfluids can be viewed as breaking $U(1)$ phase symmetry. Examples
are infinite Bose-Einstein condensates (BEC) or the paired states of a
BCS superconductor. In their mean-field description, these systems
exhibit a coherent phase at the expense of a well-defined particle
number. In finite systems the particle number is, of course, a good
quantum number. While the results we derive in this Section are well
know (see, e.g. Ref.~\cite{wen}), their derivation exhibits novel
aspects and paves the way for the description of deformed nuclei
within an EFT.

In finite superfluids the low-lying excitations are governed by two
energy scales. These are the chemical potential $\mu$ (the energy
needed to add a single boson to the system), and the energy $\Omega$
of long wave-length excitations. These scales are different for
noninteracting and for interacting systems. We consider harmonically
trapped bosons as an example. For noninteracting bosons, the proper
thermodynamic limit is defined~\cite{dalfovo1999} by keeping the
product of particle number $N$ and the third power of the trap
frequency constant. That frequency, in turn, defines $\mu$ for the
condensate. Hence $\mu \sim \Omega \sim N^{-1/3}$ scale similarly. As
$N \to \infty$, the ground states of systems with different particle
numbers $N$ become quasi degenerate, superpositions of such states
(i.e., states with a constant phase and undefined $N$) describe the
superfluid, and the $U(1)$ phase symmetry is broken in the
thermodynamic limit.  For finite uniform systems, both energies scale
as $\mu\sim\Omega\sim N^{-2/3}$~\cite{huang1987}, and the arguments
apply likewise.

For interacting Bose gases of volume $V$, the chemical potential $\mu$
typically approaches in the thermodynamic limit the nonzero value $\mu
\sim g N/V$ where $g$ measures the strength of the
interaction~\cite{dalfovo1999}. Similarly, for BCS superconductors,
adding a pair roughly costs an energy of the order of the Fermi
energy. The latter becomes constant in the thermodynamic limit. In
both cases the breaking of $U(1)$ symmetry can be understood in the
framework of the grand canonical ensemble. The system is coupled to a
particle reservoir with external chemical potential $\mu_{\rm ext}$,
and the Hamiltonian $H-\mu_{\rm ext}N$ is minimized. Adjusting
$\mu_{\rm ext}$ such that $\mu_{\rm ext}\approx\mu$ for $N \gg 1$
introduces a quasi degeneracy between states of different particle
numbers, and a superposition of these states then breaks $U(1)$
symmetry. However, the canonical and the grand canonical ensemble
yield different results. For nonextensive quantities the differences
decrease as $N^{-1/2}$ for $N \to \infty$. It is only within this
uncertainty that an isolated finite system can be viewed as equivalent
to a finite system coupled to a particle reservoir. Technically, the
introduction of a finite external chemical potential $\mu_{\rm ext}$
breaks time-reversal invariance, and the resulting effective theory
differs from the case $\mu_{\rm ext} = 0$. As we will see, the
latter can be recovered from the former by simply setting $\mu_{\rm
  ext} = 0$ in the leading-order equations we derive below. For
nonzero $\mu_{\rm ext}$, the low-energy scales of interest are $\mu
-\mu_{\rm ext}$ and $\Omega$, and we assume that both are similar in
size.

In the case of a broken $U(1)$ symmetry we have ${\cal G} = U(1)$ and
${\cal H} = 1$. The Nambu-Goldstone fields parameterize the coset
${\cal G}/{\cal H}\sim {\cal G}$, which is the group itself. The
fields induce local phase transformations, and in a finite system the
relevant operator is
\begin{equation}
U(\alpha,\beta) = e^{i\alpha(t)}e^{iV^{1/2}\beta(t,\vec{x})} \ .
\label{uni}
  \end{equation}
Here, $\beta$ is the Nambu-Goldstone field, the angle $\alpha$ is the
zero mode that needs to be singled out, and $V$ is the volume. For a
proper Nambu-Goldstone field we have $\int_V {\rm d}^3x \beta = 0$.
Following Refs.~\cite{weinbergbook,brauner2010}, we build the
invariants of our theory from the derivatives ($\nu = x, y, z$)
\begin{eqnarray}
-i U^{-1} \partial_t U &=& \dot\alpha +V^{1/2} \dot\beta \ , \\
-i U^{-1} \partial_\nu U &=& V^{1/2} \partial_\nu \beta
  \end{eqnarray}
which are in the Lie algebra of ${\cal G}$. Here the dot denotes the
time derivative. Under a global phase transformation with angle
$\gamma$, the operator $U$ becomes $e^{i\gamma}U(\alpha,\beta) =
U(\alpha+\gamma,\beta)$. Thus, $\dot\alpha$, $\dot\beta$, and
$\partial_\nu \beta$ are invariant under global phase transformations.
Note also that $\beta$ is a truly ``intrinsic'' degree of freedom
because it is unaffected by a global phase transformation.

We use these invariants to construct the leading-order terms in the
effective Lagrangian, taking account of energy scales. The scale
associated with the $\alpha$ degree of freedom is $\mu - \mu_{\rm
    ext}$, that associated with $\beta$ is
$\Omega$. Assuming also invariance under rotations, we have the
invariants
\begin{eqnarray}
L_0 &\equiv&  C_1\mu_{\rm ext}\int\limits_V {\rm d}^3x \dot\alpha
=  C_1V\mu_{\rm ext}{\dot\alpha} \ , \nonumber \\
L_1 &\equiv& {C_1\over 2} \int\limits_V {\rm d}^3x \dot\alpha^2 =
{C_1V\over 2}{\dot\alpha}^2 \ , \nonumber\\
L_2 &\equiv& {C_2 \over 2} \int\limits_V {\rm d}^3x \dot\beta^2 \ , \
L_3 \equiv {D \over 2} \int\limits_V {\rm d}^3x \left(\nabla \beta
\right)^2 \ .
  \end{eqnarray}
Here, $C_1$, $C_2$, and $D$
are constants that have to be determined from low-energy data. The
Lagrangian is
\begin{equation}
\label{eff}
L = {C_1 V \over 2}{\dot\alpha}^2 + C_1 V \mu_{\rm ext} \dot\alpha +
\int\limits_V {\rm d}^3x \left( {C_2\over 2} \dot\beta^2 -{D\over 2}
\left(\nabla\beta\right)^2 \right) \ .
  \end{equation}
The conserved quantity corresponding to invariance under global phase
transformations is the particle number
\begin{equation}
\label{nparticle}
N\equiv p_\alpha \equiv {\partial L\over \partial\dot\alpha} \ .
  \end{equation}
This conserved quantity can be derived via the Noether theorem, and
$p_\alpha$ is the canonical momentum of $\alpha$.

We expand $\beta(t,\vec{x}) = \sum_j \beta_j \phi_j(\vec{x})$ in a set
of orthonormal complex functions $\phi_j(\vec{x})$, $j=1, 2, \ldots$
with $\int_V {\rm d}^3 x \phi_j(\vec{x}) = 0$ (absence of zero modes
for the Nambu-Goldstone field). As an example, we choose the
eigenfunctions of a free particle in a spherical cavity of volume $V$
with von Neumann boundary conditions. The Lagrangian becomes
\begin{equation}
L = {C_1 V\over 2} {\dot\alpha}^2 + C_1 V \mu_{\rm ext}\dot\alpha +
\sum_{j>0}\left({C_2\over 2}{\dot\beta}_j^2 - {D k_j^2\over 2}\beta_j^2
\right) \ .
  \end{equation}
Here, $k_j^2$ denotes the squared momentum of the spherical wave
$\phi_j(\vec{x})$. A Legendre transformation with $p_{\beta j} \equiv
\partial L / \partial {\dot\beta}_j$ and Eq.~(\ref{nparticle}) yield
the Hamiltonian
\begin{equation}
H = {\left(p_\alpha - C_1V\mu_{\rm ext}\right)^2\over 2C_1V} + \sum_{j>0}
\left({p_{\beta j}^2\over 2C_2}  + {D k_j^2 \over 2}\beta_j^2 \right) \ .
  \end{equation}
We quantize $H$ by putting $p_\alpha = -i\partial_\alpha$, and
$p_{\beta j} = -i \partial_{\beta_j}$. Then $p_\alpha e^{iN\alpha} = N
e^{iN\alpha}$. The intrinsic degrees of freedom $\beta_j$ yield
harmonic-oscillator spectra with energies
\begin{equation}
\label{evib}
\omega_j \equiv k_j \left({D\over C_2}\right)^{1/2} \ .
  \end{equation}
In the long-wave-length limit we have $k_j \sim V^{-1/3}$, and the
measurement of the low-energy collective excitations of the superfluid
determines the ratio $D/C_2$. The amplitude of the Nambu-Goldstone
modes (i.e., the oscillator length) is
\begin{equation}
\label{osclen}
l_j\equiv \left(C_2Dk_j^2\right)^{-1/4} \ ,
  \end{equation}
and this dimensionless quantity is assumed to be small, $l_j \ll 1$.

For fixed particle number $N$, the $\alpha$-dependent part of the
Hamiltonian yields the energy $E_N \equiv \left(N-C_1V\mu_{\rm ext}
\right)^2 / (2C_1V)$. For $N \gg 1$ the energy difference of a system
with $N + 1$ and one with $N$ particles is
\begin{equation}
\label{chempot}
E_{N+1}-E_N \approx {N \over C_1 V} - \mu_{\rm ext} \ .
  \end{equation}
We note that $E_{N+1}-E_N = \mu-\mu_{\rm ext}$ on physical
  grounds.  The emergent breaking of $U(1)$ requires
  $E_{N+1}-E_N\approx 0$, and Eq.~(\ref{chempot}) relates the constant
  $C_1$ to the chemical potential and to the density of the
  system. We see that the chemical
potential and the frequencies and amplitudes of the quantized
collective vibrations determine the low-energy constants of the EFT.
In addition to the collective vibrations with frequencies~(\ref{evib})
one finds approximately equidistant levels with spacing
  $\mu - \mu_{\rm ext}$ belonging to superfluids with different
particle numbers. In superfluid atomic nuclei, these harmonic
excitations belonging to different numbers of paired nucleons are
known as pairing vibrations~\cite{Bes1966}. In summary we have shown
that in finite superfluids, the spectra of systems with different
particle numbers are related to each other, and this is a
model-independent prediction of the EFT.

We turn to higher-order corrections
and establish our power counting.  The energy scales used in the
construction of $L$ in Eq.~(\ref{eff}) are (i) the scale $\mu -
\mu_{\rm ext}$ associated with a change in particle number (we have
assumed $\dot\alpha \sim \mu-\mu_{\rm ext}$) and (ii) the scale
$\Omega$ of the collective vibrations. We have assumed that both
scales are of similar size, $\mu - \mu_{\rm ext} \sim \Omega$.
Moreover, in the low-energy domain the amplitudes of the
Nambu-Goldstone modes $\beta_j$ given by Eq.~(\ref{osclen}) obey $l_j
\ll 1$. Thus, we have
\begin{eqnarray}
\label{scalingu1}
\beta_j \sim l_j \ , \ {\dot\beta}_j \sim \Omega l_j \ , C_2 \sim
\Omega^{-1} l_j^{-2} \ , \ p_{\beta j} \sim l_j^{-1} \ .
  \end{eqnarray}
In the construction of terms of higher order, we consider only powers
of single derivatives because higher-order derivatives of the
Nambu-Goldstone field $\beta$ and of the zero mode $\alpha$ can be
eliminated via perturbative field
redefinitions~\cite{damour1991,grosseknetter1994}. The EFT has a
breakdown scale $\Lambda$ with $\Lambda \gg \Omega, |\mu_{\rm ext}
-\mu|$. At this scale, $\dot\alpha$ (and $\dot\beta$, $\partial_\nu
\beta$, $\beta$) are a factor $\sqrt{\Lambda/|\mu_{\rm ext}-\mu|}$
(and $\sqrt{\Lambda/\Omega}$, respectively) larger than in the
low-energy domain. At the breakdown scale, a Lagrangian term of the
form $C_{pqr}{\dot\alpha}^p{\dot\beta}^q \beta^r$ must scale as
$\Lambda$ when written in terms of such velocities and fields. This
determines the scaling of $C_{pqr}$.  When evaluated for velocities
and fields in the low-energy domain, that term yields a contribution
to $L$ of order $\Lambda (|\mu_{\rm ext} - \mu| / \Lambda)^{p/2}
(\Omega/\Lambda)^{q/2+r/2}$. Similar considerations apply for terms
containing powers of the spatial derivatives $\partial_\nu\beta$. For
energies below the breakdown scale, all these higher-order corrections
are perturbatively small. Thus, our procedure yields a perturbative
expansion in powers of $|\mu_{\rm ext}-\mu|/\Lambda$ and
$\Omega/\Lambda$.

\section{Emergent breaking of rotational symmetry}

For finite objects with an axially symmetric ground-state deformation,
the emergent symmetry breaking is from the rotation group to the
subgroup of axial symmetry. It is, therefore, useful to recall the
breaking of spin-rotational symmetry $O(3)$ to $O(2)$ in infinite
ferromagnets~\cite{Leutwyler1994,roman1999, Baer2004}. In the ground
state, all spins point in the same direction, violating the
spin-rotational symmetry of the Hamiltonian. Ground states with
macroscopically different spin directions have zero overlap and define
inequivalent Hilbert spaces. As a consequence, the low-lying spectrum
of the ferromagnet is dominated by Nambu-Goldstone modes, i.e., spin
waves of long wave length that locally induce small rotations of the
aligned spins. For a ferromagnet of finite size~\cite{Hasenfratz1993},
the formerly inequivalent Hilbert spaces are connected by
non-vanishing tunneling matrix elements, the ground states belonging
to different spin directions have nonzero overlap, and there exist
states that are superpositions of these ground states. Such states
are, for instance, the Wigner $D$-functions for rotational motion.
Physically that implies that the ferromagnet may rotate about an axis
perpendicular to the direction of the aligned spins. An analogous
situation occurs in linear molecules and in axially symmetric deformed
even-even nuclei. In the limit of infinitely large size, the low-lying
parts of the spectra of these systems would be determined by
Nambu-Goldstone modes. But finite linear molecules or finite axially
symmetric nuclei may rotate about an axis perpendicular to the
symmetry axis. In all three cases, rotational motion occurs as a
consequence of emergent symmetry breaking due to the finite size of
the system. Rotational motion is, therefore, distinctly different from
the Nambu-Goldstone modes. It plays the same role as the zero mode
$\alpha$ in Eq.~(\ref{uni}) for the emergent breaking of $U(1)$ phase
symmetry.

For a deformed nucleus with $A \gg 1$ nucleons, we can quantify these
statements. The linear extension of the system is $\propto A^{1/3}$.
The moment of inertia is proportional to mass $\times$ (length)$^2
\propto A^{5/3}$. With increasing $A$, the frequency of rotational
motion tends to zero like $A^{- 5/3}$, faster than the wave number
$\propto A^{- 1/3}$ of the massless modes (analogues of the
Nambu-Goldstone modes in the infinite system). Hence, there exists a
regime of $A$ values where the scale $\xi$ for rotational motion is
small compared to the scale $\omega$ for vibrational motion and where
$\omega$ in turn is small compared to the breakdown scale $\Lambda$ of
the EFT defined by pair-breaking excitations~\cite{dean2003}. In the
regime $\xi \ll \omega \ll \Lambda$, it is meaningful to consider
rotational motion as a small (in energy) correction to the Goldstone
theorem. That is the regime we study here. For instance, in rare-earth
nuclei we have $\xi \approx 80$~keV from the lowest energy spacing in
a rotational band, and $\omega \approx 1$~MeV from the lowest
``vibrational'' band head, while $\Lambda \approx 2-3$~MeV is the cost
of pair breaking. The condition $\omega \ll \Lambda$ is met only
marginally. Corresponding considerations apply to linear molecules
such as CO$_2$. Here the scale of rotational energies is $\xi
\approx$~1~cm$^{-1}$, that of vibrational energies is $\omega
\approx$~500~cm$^{-1}$, while the breakdown scale
$\Lambda\approx$~10,000~cm$^{-1}$ is defined by electronic
excitations. The conditions $\xi \ll \omega \ll \Lambda$ are very well
fulfilled.

Axially symmetric deformed even-even nuclei consist of nucleons, and
linear molecules consist of nuclei. The locations of these
constituents have body-fixed spherical coordinates $r, \theta, \phi$.
Vibrations of the nucleus/molecule about an axis perpendicular to the
symmetry axis act locally on the constituents. The resulting
dislocations are assumed to have small amplitude. We expect that
Nambu-Goldstone modes related to the coordinate $r$ have higher
frequencies than those due to $\theta$ and $\phi$. We, therefore,
confine attention to the latter variables although our approach can be
straightforwardly generalized. In addition to these small-amplitude
vibrations we also consider global rotations of the entire
nucleus/molecule.

Our effective theory is universal and applies both to linear molecules
and to deformed nuclei. Apart from the magnitude of the low-energy
constants, the key difference between both systems is in the symmetry
properties of the deformed ground--state wave functions as these
define the symmetry properties of the admissible low--lying
Nambu-Goldstone modes~\cite{weinbergbook}. We assume that molecular
and nuclear ground-state wave functions are axially symmetric about
the body-fixed $z'$-axis, invariant under time-reversal, and have
positive parity. As a consequence of pairing (superfluidity),
even-even nuclei differ from linear molecules in that their intrinsic
ground states are also invariant under rotations by $\pi$ about any
axis perpendicular to the symmetry axis, i.e., possess positive ${\cal
  R}$ parity~\cite{bohr1958,bmbook}. Hence, low-energetic intrinsic
excitations in nuclei must also have positive ${\cal R}$ parity.

\subsection{Dynamical Variables and Power Counting}

As done for the case of the
ferromagnet~\cite{Leutwyler1994,roman1999,kampfer2005} and in
Sect.~\ref{emer}, we consider the Nambu-Goldstone modes as classical
fields that are later quantized.  We prefer to write these fields in
the space-fixed (rather than the body-fixed) coordinate system because
here the commutation relations of the three generators $J_x, J_y, J_z$
of infinitesimal rotations about the space-fixed $x, y, z$ axes,
respectively, are of standard form. The molecular/nuclear ground state
is invariant under $SO(2)$ rotations about the body-fixed $z'$-axis
while $SO(3)$ symmetry is broken by the deformation. Therefore, the
Nambu-Goldstone modes lie in the coset space $SO(3) /
SO(2)$~\cite{coleman1969,callan1969,brauner2010}. The modes depend on
the angles $\theta, \phi$ defined above and on time $t$, and are
generated by a unitary transformation $U$. As in Eq.~(\ref{uni}) we
parameterize the matrix $U$ in product form,
\begin{eqnarray}
\label{1}
U &=& g(\alpha, \beta) u(x,y) \ , \nonumber\\
g(\alpha, \beta) &=& \exp \left\{ -i \alpha(t) J_z \right\}
\exp\left\{ -i \beta(t) J_y \right\} \ , \nonumber\\
u(x,y) &=& \exp\left\{-i x(\theta, \phi, t) J_x
-i y(\theta, \phi, t) J_y \right\} \ .
  \end{eqnarray}
The purely time-dependent variables $\alpha(t)$ and $\beta(t)$ are the
zero modes. They describe global rotations of the finite system and
are factored out~\cite{Leutwyler1987}. They are not Nambu-Goldstone
modes but upon quantization generate rotational
bands~\cite{chandrasekharan2008,papenbrock2011}. The fields $x(\theta,
\phi, t)$ and $y(\theta, \phi, t)$ with $|x|, |y| \ll 1$ generate
small-amplitude vibrations of the constituents. These depend
non-trivially on $\theta$ and $\phi$ so that
\begin{equation}
\label{2}
\int {\rm d} \Omega \ x(\theta, \phi, t) = 0 = \int {\rm d} \Omega \
y(\theta, \phi, t) \ .
  \end{equation}
Here ${\rm d} \Omega$ is the surface element of the three-dimensional
unit sphere. In an infinite system $x(\theta, \phi, t)$ and $y(\theta,
\phi, t)$ would be genuine Nambu-Goldstone modes.  Eqs.~(\ref{1}) and
(\ref{2}) define the dynamical variables of the system. Eq.~(\ref{1})
may look like a rather special ansatz but actually follows from the
most general form of $U$.

Further progress hinges on the identification of the energy scales
$\xi$, $\omega$, and $\Lambda$ defined above. The ranges of the
variables $\alpha$ and $\beta$ being of order unity, the ratios
$\dot{\alpha} / \alpha$ and $\dot{\beta} / \beta$ are governed by the
energy scale $\xi$ of rotational motion. We have $|x|, |y| \ll 1$,
indicating that the amplitudes of the Nambu-Goldstone fields are small.
Then $|\dot{x}|\sim \omega |x|$ and $|\dot{y}| \sim \omega|y|$.  We
are going to show that power counting based upon the inequalities $\xi
\ll \omega \ll \Lambda$ together with the symmetry requirements
formulated above uniquely determine the leading-order part of the
Hamiltonian, except for a small number of constants that have to be
determined by fits to data.

Our EFT is characterized by two breakdown scales. The first scale is
set by $\Lambda \gg \omega$ and marks the appearance of neglected
degrees of freedom. In deformed nuclei these are single-particle
degrees of freedom or pair-breaking effects and in linear molecules,
electronic excitations. The second scale is related to large-amplitude
excitations of the Nambu-Goldstone fields. That scale is reached when
these excitations are so large that they practically restore the
spherical symmetry of the intrinsically deformed object, or when the
energy due to the zero-mode velocities $\dot\alpha$ and $\dot\beta$
reaches the vibrational scale $\omega$. This second scale is set by
$\omega^2/\xi$. For well-deformed nuclei, that scale considerably
exceeds $\Lambda$, while both scales are similar in size for linear
molecules. We now discuss both breakdown scales separately.

At energies below $\Lambda$ the neglected degrees of freedom cause the
appearance of higher-order terms in the effective Lagrangian of the
EFT. Such terms involve powers of the leading-order fields and
velocities and, possibly, also higher derivatives. The latter can be
eliminated via perturbative field
redefinitions~\cite{damour1991,grosseknetter1994} and are not
considered here. At the breakdown scale (where the amplitudes and
velocities of the Nambu-Goldstone fields are a factor
$(\Lambda/\omega)^{1/2}$ larger than in the low-energy domain), a term
in the effective Lagrangian with $\tau$ velocities and $n$ powers of
$x$ or $y$ is of order $\Lambda$. At the low-energy scale $\omega$
that term yields a contribution of order $\omega (\omega /
\Lambda)^{(n+\tau)/2-1}$. Terms like that give rise to small corrections
and, at each order, are finite in number. Similar considerations apply
to the spatial derivatives.

The scale for the breaking of emergent symmetry is reached when the
amplitudes, velocities and higher derivatives are a factor $(\omega /
\xi)^{1/2}$ larger than in the low-energy domain. The arguments of the
previous paragraph can essentially be repeated by replacing the scale
$\Lambda$ by $\omega^2/\xi$, and the ratio $\Lambda / \omega$ by
$\omega / \xi$. At the breakdown scale the contribution of a term in
the effective Lagrangian which contains $\sigma$ velocities is of
order $\omega$. At the low-energy scale $\xi$ that term scales as $\xi
(\xi / \omega)^{\sigma/2-1}$.

These arguments establish our power counting. As a result, our EFT
provides an expansion in the two small parameters $\omega / \Lambda$
and $\xi / \omega$, and there is a finite number of terms for each
power of these parameters.

\subsection{Effective Lagrangian}

The effective Lagrangian is built from invariants. These are
constructed from elements $a_\mu^x, a_\mu^y, a_\mu^z$ defined by
\begin{eqnarray}
\label{3}
U^{-1} i \partial_\mu U &=& a_\mu^x J_x + a_\mu^y J_y + a_\mu^z J_z \ .
  \end{eqnarray}
The symbol $\partial_\mu$ with $\mu = 1, 2, 3$ stands for the partial
derivatives with respect to the angles $\theta, \phi$ and time $t$.
Explicit expressions for these elements are obtained from
Eqs.~(\ref{1}) and (\ref{3}) in terms of a series expansion in powers
of $x$, $y$, and their partial derivatives where only leading-order
terms are kept. We use the Baker--Campbell--Hausdorff
expansion and obtain
\begin{eqnarray}
a^x_t &=& \dot{x} +{y \over 6} (x \dot{y} - y \dot{x}) -\dot{\alpha}
\sin \beta - y \dot{\alpha} \cos \beta + \ldots\ , \nonumber \\
a^y_t &=& \dot{y} -{x \over 6} (x \dot{y} - y \dot{x}) + \dot{\beta}
+ x \dot{\alpha} \cos \beta + \ldots \ , \\
\label{3a}
a^z_t &=& -{1 \over 2} (x \dot{y} - y \dot{x}) + \dot{\alpha} \cos
\beta - y \dot{\alpha} \sin \beta - x \dot{\beta} + \ldots \ , \nonumber
  \end{eqnarray}
and
\begin{eqnarray}
a^x_\nu &=& \partial_\nu{x} +{y \over 6} (x \partial_\nu{y} - y \partial_\nu{x}) + \ldots\ , \nonumber \\
a^y_\nu &=& \partial_\nu{y} -{x \over 6} (x \partial_\nu{y} - y \partial_\nu{x}) + \ldots \ , \\
\label{3b}
a^z_\nu &=& -{1 \over 2} (x \partial_\nu{y} - y \partial_\nu{x}) + \ldots \ , \nonumber
  \end{eqnarray}
To define the invariants we calculate the changes induced on the
variables $\alpha, \beta, x, y$ by infinitesimal rotations $r$ of $U$
about angles $\delta \chi_k$ around the space-fixed $k = x, y, z$
axes.  With $h(\gamma) = \exp \{ - i \gamma J_z \}$ we have from
Eq.~(\ref{1}) that
\begin{equation}
\label{rot1}
r g(\alpha, \beta) = g(\alpha', \beta') h(\gamma') \ .
  \end{equation}
Here, the primed angles depend on the angles of the rotation $r$ and
the angles $\alpha$, $\beta$.  The right-hand side of Eq. (\ref{rot1})
has the form of a general rotation with Euler angles ($\alpha',
\beta',\gamma'$), with explicit expressions given in
Ref.~\cite{papenbrock2011}. Thus,
\begin{eqnarray}
r U
&=&g(\alpha', \beta') h(\gamma') u \nonumber\\
&=&g(\alpha', \beta') \ [ h(\gamma') u h^\dag(\gamma') ]\ h(\gamma') \ .
  \end{eqnarray}
As a result we find that the angles $\alpha$ and $\beta$ transform
nonlinearly as the azimuthal and polar angle of the two-sphere,
respectively, while $x$ and $y$ transform linearly as the $x$ and $y$
components of a vector under rotations around the $z$ axis. In other
words, under a rotation, the nucleus as a whole changes its
orientation and undergoes a rotation around its symmetry axis. This
transformation behavior under rotations confirms that $\alpha$ and
$\beta$ describe the global orientation of the axially symmetric
nucleus, while $x$ and $y$ are ``intrinsic'' degrees of freedom.  Thus
any combination of $x$ and $y$ that formally exhibits axial symmetry
is indeed fully invariant under rotations. For example, $x^2 + y^2$ is
invariant under rotations, and the four quantities $x, y, \dot{x},
\dot{y}$ are transformed into linear combinations of $x', y',
\dot{x}', \dot{y}'$. These transformation properties are
characteristically different from the ones for an infinite system
where in Eq.~(\ref{1}) we would have $g(\alpha, \beta) = 1$.

Time-reversal invariance requires that invariants involving time
derivatives must contain even powers of $a^x_t, a^y_t, a^z_t$. The
lowest-order invariants obtained from Eqs.~(\ref{3a}) are
\begin{eqnarray}
\label{4}
{\cal L}_{1a} &=& \dot{\beta}^2 + \dot{\alpha}^2 \sin^2 \beta \ ,
\nonumber \\
{\cal L}_{1b} &=& \dot{x}^2 + \dot{y}^2 + 2(x \dot{y} - y \dot{x})
\dot{\alpha} \cos \beta \ , \nonumber \\
{\cal L}_{1c} &=& (x \dot{y} - y \dot{x})^2 \ , \nonumber \\
{\cal L}_{1d} &=& (x^2 +y^2)[\dot{x}^2 + \dot{y}^2 + 2(x \dot{y} - y
\dot{x}) \dot{\alpha} \cos \beta] \ .
  \end{eqnarray}
We note that the invariant ${\cal L}_{1a}$ is essentially the
Lagrangian of a rotor, and that the Lagrangian density ${\cal L}_{1b}$
couples global rotations to the Nambu-Goldstone modes.  The invariant
${\cal L}_{1 c}$ is related to the angular momentum of the
Nambu-Goldstone modes, see Eq.~(\ref{K}).  The invariant ${\cal
  L}_{1d}$ is obtained by multiplying ${\cal L}_{1b}$ with the
invariant $(x^2 + y^2)$ and is of the same order as ${\cal L}_{1c}$.
As for the invariants constructed from $a^x_\nu$, $a^y_\nu$ $a^z_\nu$
with $\nu = \theta$ or $\nu = \phi$, we use that for fixed $\nu$ and
$\nu'$, the forms $(a^x_\nu)^2 + (a^y_\nu)^2$, $a^z_\nu$ and $a^z_\nu
a^z_{\nu'}$ are invariant. Admissible linear combinations of these
expressions are defined by the requirement of axial symmetry.
Suppressing terms of higher order than $x^4$ and multiplying with the
additional invariant $(x^2 + y^2)$, we find the invariants
\begin{eqnarray}
\label{5}
{\cal L}_{2a} &=& (\vec{\bf L} x)^2 + (\vec{\bf L} y)^2 \ , \nonumber \\
{\cal L}_{2b}&=& ({\bf L}_z x)^2 + ({\bf L}_z y)^2 \ , \nonumber \\
{\cal L}_{2c} &=& (x \vec{\bf L} y - y \vec{\bf L} x)^2 \ , \nonumber \\
{\cal L}_{2d} &=& (x^2 + y^2) \left( (\vec{\bf L} x)^2 + (\vec{\bf L}
y)^2 \right) \ .
  \end{eqnarray}
Here $\vec{\bf L}$ (${\bf L}_z)$ is the vector (the $z$-component)
operator of orbital angular momentum, respectively, written in terms
of $\theta$ and $\phi$~\cite{varshalovich1988}. The occurrence of the
term ${\cal L}_{2 b}$ reflects the fact that we impose only axial
rather than rotational symmetry on ${\cal L}$. The Lagrangian $L$ is
given by
\begin{eqnarray}
\label{6}
L &=& L_1 + L_2 \nonumber\\
&=& \sum_{i = a, b, c, d} \int {\rm d} \Omega \ \bigg(
\frac{C_i}{2} {\cal L}_{1 i} - \frac{D_i}{2} {\cal L}_{2 i} \bigg) \ .
  \end{eqnarray}
Here $C_i$ and $D_i$ with $i = a, b, c, d$ are constants that are
determined by low-energy data.

We expand the real variable $x$ as
\begin{eqnarray}
\label{7}
x = \sum_{\lambda = 2}^\infty \sum_{\mu = -\lambda}^\lambda
x_{\lambda \mu} Z_{\lambda \mu}
  \end{eqnarray}
and correspondingly for $y$, $\dot{x}$, $\dot{y}$. Aside from
normalization constants, the real orthonormal functions $Z_{\lambda
  \mu}$ are equal to the real part (for $\mu\ge 0$) and imaginary part
(for $\mu<0$) of the spherical harmonics $Y_{\lambda \mu}$. The
coefficients $x_{\lambda\mu}$ are real. Terms with $\lambda = 0$ and
$\lambda = 1$ are excluded since $\lambda = 0$ violates Eq.~(\ref{2})
and describes global rotations while $\lambda = 1$ describes
translations in space~\cite{bmbook}. We insert the
expansions~(\ref{7}) into Eq.~(\ref{6}) and use the resulting
expression for $L$ to define the real canonical momenta
\begin{eqnarray}
\label{8}
p_\beta &=& \frac{\partial L}{\partial \dot{\beta}} \ , \ p_\alpha
= \frac{\partial L}{\partial \dot{\alpha}} \ , \nonumber\\
p^x_{\lambda \mu} &=& \frac{\partial L}{\partial \dot{x}_{\lambda \mu}} \ , \
p^y_{\lambda \mu} = \frac{\partial L}{\partial \dot{y}_{\lambda \mu}} \ .
  \end{eqnarray}
From Noether's theorem we obtain explicit expressions for the three
components $I_x, I_y, I_z$ of angular momentum. These are
\begin{eqnarray}
I_x &=& - p_\beta \sin \alpha - p_\alpha \cot \beta \cos \alpha +
{\cos \alpha \over \sin \beta} \ K \ , \\
I_y &=& p_\beta \cos \alpha - p_\alpha \cot \beta \sin \alpha + {\sin
\alpha \over \sin \beta} \ K \ , \\
I_z &=& p_\alpha \ .
  \end{eqnarray}
Here
\begin{eqnarray}
\label{K}
K = \int {\rm d} \Omega \ (x p_y - y p_x )
  \end{eqnarray}
is the angular momentum of the two-dimensional oscillators that
describe the intrinsic vibrations. The square of the total angular
momentum is
\begin{equation}
\label{9}
I^2 = p_\beta^2 +{1 \over \sin^2\beta} \left(p_\alpha^2 -2Kp_
\alpha \cos \beta + K^2 \right) \ .
  \end{equation}
The terms in Eqs.~(\ref{9}) obtained by putting $K = 0$ can be shown
to be equal to the square of the total angular momentum of the rigid
rotor.

\subsection{Effective Hamiltonian and Quantization}

We use Eqs.~(\ref{8}) and the standard Legendre transformation to
transform the effective Lagrangian $L$ into the effective Hamiltonian
$H$. The scales of the coefficients $C_i$ and $D_i$ (and the terms
that are kept in $H$) are determined by assuming $C_a {\cal L}_{1a}
\sim \xi$, $C_b {\cal L}_{1b} \sim \omega$, $C_c {\cal L}_{1c}, C_d
{\cal L}_{1d} \sim |x|^2 \omega$, $p_\beta, p_\alpha \sim 1$, and
$p_x, p_y \sim |x|^{-1}$. Relevant parts of $H$ are given
below.

For the rigid-rotor part of $H$, quantization is achieved by
symmetrization with respect to $\alpha, \beta$ and by putting $p_\beta
= -i (\sin \beta)^{- 1/2} \partial_\beta (\sin \beta)^{1/2}$,
$p_\alpha = -i \partial_\alpha$. This is the usual quantization
  for a particle on the sphere. For the remaining canonical momenta
we have
\begin{equation}
\label{10}
p^x_{\lambda \mu} = - i \frac{\partial}{\partial x_{\lambda \mu}} \ , \
p^y_{\lambda \mu} = - i \frac{\partial}{\partial y_{\lambda \mu}} \ .
  \end{equation}
Substitution of these expressions into Eq.~(\ref{9}) yields the
quantized form of the square $\hat{I}^2$ of the operator of total
angular momentum. The operator $\hat{K}$ is given by $\hat{K} =
\sum_{\lambda\mu} \left(x_{\lambda\mu} p^y_{\lambda\mu} -
y_{\lambda\mu} p^x_{\lambda \mu} \right)$. The three components
$\hat{I}_x$, $\hat{I}_y$, $\hat{I}_z$ of the quantized angular
momentum can be shown to obey the standard commutation
relations. Moreover, every component commutes with $\hat{K}$ and with
the quantized Hamiltonian $\hat{H}$. A complete set of commuting
operators is, therefore, $\hat{I}^2, \hat{I}_z, \hat{K}, \hat{H}$.

\subsection{Spectra}
The leading-order (${\cal O}(\omega)$) contribution
to $\hat{H}$ is given by
\begin{eqnarray*}
\hat{H}_\omega = \sum_{\lambda \mu} \left( {(p^x_{\lambda
\mu})^2 + (p^y_{\lambda \mu})^2 \over 2C_b} + {C_b\over 2} \omega^2_{\lambda \mu}
\left( x_{\lambda \mu}^2 + y_{\lambda \mu}^2 \right) \right)
  \end{eqnarray*}
and describes a set of uncoupled harmonic oscillators with frequencies
$\omega_{\lambda \mu} = [(\lambda ( \lambda + 1) D_a + \mu^2 D_{b}) /
  C_b]^{1/2}$. We combine the degrees of freedom $x_{\lambda \mu}$ and
$y_{\lambda \mu}$ into a two-dimensional $SO(2)$ symmetric harmonic
oscillator with quantum numbers $n_{\lambda \mu} = 0, 1, 2, \ldots$
and $k_{\lambda \mu} = 0, \pm 1, \pm 2, \ldots$. In units of
$\omega_{\lambda \mu}$ the energies are $2n_{\lambda \mu} +
|k_{\lambda \mu}| +1$. The operator $\hat{K}$ has eigenvalues $K =
\sum_{\lambda\ge 2} \sum_{\mu = - \lambda}^\lambda k_{\lambda \mu}$.

Next-order corrections to the Hamiltonian $\hat{H}_\omega$ are either
of order ${\cal O}(\xi)$ or of order ${\cal O}(x^2\omega)$.
The former couple rotations to vibrations. The latter add
anharmonicities to the harmonic vibrations and thereby lift the
degeneracies. We confine ourselves here to the former. The Hamiltonian
is
\begin{equation}
\label{12}
\hat{H}_{\omega,\xi}= \hat{H}_{\omega} + {\hat{I}^2 - \hat{K}^2
\over 2C_a} \ ,
  \end{equation}
with $\hat{I}^2$ given by the quantized form of Eq.~(\ref{9}). The
eigenfunctions of $(\hat{I}^2 - \hat{K}^2)$ are Wigner $D$-functions
$D_{M,K}^I(\alpha,\beta,0)$ that depend on total integer spin $I$ and
its projections $-I \le M, K \le
I$~\cite{varshalovich1988,papenbrock2011}. The eigenvalues of
$\hat{I}^2$ are $I (I + 1)$ with $I \ge |K|$. We see that each
vibrational state of the leading-order Hamiltonian $\hat{H}_\omega$
becomes a band head with spin $|K|$, and the spectrum exhibits a
rotational band on top of each band head.  At this order in the EFT,
all rotational bands have the same moment of inertia $C_a$.
Differences in the moments of inertia are higher-order
effects~\cite{Zhang2013}.

The ground state has quantum numbers $n_{\lambda \mu} = 0$,
$k_{\lambda \mu} = 0$ (this implies $K = 0$) and spin $I = 0$. It has
positive parity and, in the case of nuclei, positive ${\cal R}$
parity. For nuclei, this limits the rotational states in the
ground-state band to even values of $I$.

We turn to excited states. Here nuclei and linear molecules differ.
For $D_b > 0$ the lowest single vibrational excitation corresponds to
the mode $(x_{2 0}, y_{2 0})$. The fields $x$ and $y$ have positive
parity, in keeping with the parity of the axial vectors $J_x$ and
$J_y$ in Eqs.~(\ref{1}). The lowest excitation has $|K| = 1$, negative
intrinsic ${\cal R}$ parity and, thus, values of $I=1,2,3,\ldots$. For
linear molecules, states with $|K| = 1$ are indeed the lowest-lying
vibrations~\cite{Herzberg1945}.  In contradistinction, in nuclei such
states are excluded because paired states have positive ${\cal R}$
parity.  Pair breaking (i.e., generation of states with odd ${\cal R}$
parity) happens only at the breakdown scale $\Lambda$ of our EFT.
Thus, pairing excludes low-lying magnetic dipole
excitations~\cite{heyde2010,bentz2011} and more generally any vibrational band
head with odd spin and positive parity in the low-energy regime.  This
essential element provides the only difference in the low-energy
spectra of molecules and of deformed nuclei. Indeed, there are no
low-lying $I^\pi = 1^+$ states in deformed even-even
nuclei~\cite{davidson1981, aprahamian2006}.

As is well known from data~\cite{davidson1981, aprahamian2006}, the
ground states of deformed even-even nuclei consistently have quantum
numbers $I = 0 = K$ and positive parity.  Low-lying vibrational states
have $K = 0$ and $|K| = 2$ and positive parity. In rare-earth nuclei,
both band heads have an excitation energy of about $1$ MeV. In the
present approach, these states are generated by local rotations around
an axis that is perpendicular to the symmetry axis of the ground
state.  However, the wavefunction corresponding to $K=0$ is symmetric
under any $SO(2)$ rotation of $x$ and $y$ and must thus be viewed as an axially
symmetric excitation that corresponds to the $\beta$ excitation of the
geometrical model~\cite{bmbook}. The $|K|=2$ wave functions exhibit no
symmetry under exchange of $x$ and $y$ and thus break the axial
symmetry. Thus, they correspond to the $\gamma$ mode of the
geometrical model.

Finally, we note that the present approach can be extended to EFTs for
other cases such as general molecules. In this case, the ground state fully
breaks $SO(3)$ to the point group ${\cal P}$ of the molecule, and the
coset space is thus $SO(3)/{\cal P}$. This case is technically simpler than
the breaking from $SO(3)$ to $SO(2)$ because the coset essentially
has a group structure. The coset can be parameterized by three Euler angles,
and the zero modes indeed describe the orientation of the molecule in space
while the Nambu-Goldstone fields describe the intrinsic vibrations.

\section{Summary}
In summary, we have shown how to develop effective field theories for
finite systems with emergent symmetry breaking. We have applied this
approach to two types of systems. (i) Superfluids like infinite
Bose-Einstein condensates (BEC) or the paired states of a BCS
superconductor that break $U(1)$ phase symmetry. (ii) Systems with
non-spherical ground states such as molecules and atomic nuclei that
break rotational symmetry. In both cases, symmetry arguments alone
yield the universal features of the low-lying excitations. In case (i)
these are vibrations. We also relate the spectra of systems with
different particle numbers. In case (ii) these are vibrations that are
the heads of rotational bands. The moment of inertia is a fit
parameter and will, in general, be different for different physical
systems (for instance, superfluid and normal systems). Nuclei and
molecules differ in that the ground states of even-even nuclei are
paired. This accounts for the absence in nuclei of low-lying band
heads with odd spin and positive parity.

In contrast to phenomenological approaches, and except for a small
number of constants, our approach yields an explicit expression for
the Hamiltonian in leading order, and a systematic procedure to
generate terms of higher order. It may, thus, be a useful starting
point for the analysis of spectra and electromagnetic transitions.  It
is textbook knowledge that the traditional collective models of
deformed nuclei~\cite{bmbook,eisenberg,iachello} overpredict
transitions between the $\beta$-band ($\gamma$ band) and the
ground-state band by a factor of about ten (four)~\cite{rowe}.
Furthermore, the interpretation of low-lying collective $0^+$ states
as $\beta$ band heads has been put into question by the available data
on $E2$ transitions~\cite{garrett2001}. This makes it very interesting
to study electromagnetic couplings within the EFT approach.

\begin{acknowledgments}
The authors thank N. Pietralla and A. Richter for discussions.
This work has been supported by the U.S.  Department of Energy under
grant Nos. DE-FG02-96ER40963 (University of Tennessee) and
DE-AC05-00OR22725 with UT-Battelle, LLC (Oak Ridge National
Laboratory), and by the Alexander-von-Humboldt Foundation.
\end{acknowledgments}