We obtain rotational (\(\Omega \ne 0\) ) states by integrating
Eq. (REF ) until we get convergence; given our initial data, the
system contains \(N\)  bosons. For \(\Omega = 0\)  and \(V_\theta =0\) , the solution
of the ARGLPE (REF ) converges, at large times, to ground states
\(\psi _0\)  that are spherically symmetric, compact objects of radius \(R\) . This
radius can be estimated by using a variational ansatz [1]} or it can
be computed numerically at both zero temperature (\(T =0\) ) [2]}, [3]} and finite temperature (\(T >0\) ) [4]}. We start from the \(T=0\)  state of
Ref. [4]}; and then we increase the value of \(\Omega \)  in steps
of 1; at each such increase in the value of \(\Omega \) , we use the converged
ARGLPE state, from the previous value of \(\Omega \) , as the initial data. In the
inset of Fig. REF  we plot the ARGLPE-converged values of
\(E-\Omega J_z\)  versus \(\Omega \) ; from this plot it is apparent that the
non-rotating state, with no vortices, loses its stability around \(\Omega \simeq 14\) , to a state with 4 vortices. We show scans in which \(\Omega \)  increases
(red lines) and decreases (green lines). The differences between these scans
indicate that this system exhibits hysteresis: It has no vortices along the
horizontal part of the red line and 4 vortices along the green ones; and it
goes from the 4-vortex branch back to the 0-vortex branch at \(\Omega \simeq 9\) . In Fig. REF  we give an illustrative plot of
\(J_z\)  versus time \(t\)  for the GPPE evolution of the 4-vortex state at
\(\Omega =14\) ; we see that \(J_z\)  is not conserved
perfectly, but it decreases slightly (\(\simeq 1\%\)  per turn). \(J_z\)  is
not conserved exactly because rotational symmetry holds only approximately,
given our PBCs; the closer we are to the center of our
computational domain the better is this rotational symmetry. We show, in Figs.
REF (a) and (b), respectively, isosurface plots of the boson
density (top view) and of \((\nabla \times (\rho v))^2\)  (side view) for the
converged solution of ARGLPE, with \({\mathcal {N}}=256\) , \(G=800\) , \(g=80 \) , and
\(\Omega = 60\) . Clearly, our ARGLPE-based algorithm yields
rotating, self-gravitating Bose-Einstein condensates (BECs), with vortices in
the GPPE; this has not been possible hitherto For example, such
states were not obtained in Ref. [6]}, probably because
pure GPPE evolution was used, without the benefit of ARGLPE-based
minimization..
<FIGURE>