First, let me say I'm not sure what is meant by a BEC with $T\gt 0$. Condensation is a finite temperature phenomenon, which occurs due to the presence of pair-wise interactions (generally attractive, but pairing can happen even for repulsive potential) in a many-body system. For instance, in a superconductor below some critical temperature $T_c$, electrons with opposite momenta and spin (s-wave pairing) pair up to form a bound state called a Cooper pair.
The ground state of the unpaired electron gas for $T \gt T_c$ is characterized by the Fermi energy $E_F$. After condensation, the many-body system has a new ground state at energy $E_{bcs} = E_F - \Delta $, where $\Delta \sim k_B T_c$ is the binding energy of a Cooper pair. $\Delta$ is also known as the gap.
For $T\lt T_c$ all the electrons are not paired up due to thermal fluctuations. However, the number of unpaired electrons as a fraction of the total number of electrons (the condensate fraction) goes as $N_{free}/N_{all} = 1- (T/T_c)^\alpha$, where $\alpha\gt 0$. The number of free electrons drops rapidly as $T$ is decreased below $T_c$. In lab setups, BEC's generally undergo some form of evaporative cooling to get rid of particles with energies greater than $\Delta$. At this point the condensate can be treated as a gas of interacting (quasi)particles (cooper pairs) with an approximate hard-core repulsion.
So the gas, before and after condensate formation, is always at finite temperature! This is reflected, for instance, in the dependance of the condensate fraction on $T$ as mentioned above.
The mean-field solutions for low-energy excitations of the condensate are given by the Gross-Pitaevskii equation(GPE):
$$ \left( - \frac{\hbar^2}{2m}\frac{\partial^2}{\partial r^2} + V(r) + \frac{4\pi\hbar^2 a_s}{m} |\psi(r)|^2 \right) \psi(r) = \mu \psi(r) $$
where $a_s$ is the scattering length for the hard-core boson interaction, with $a \lt 0$ for an attractive interaction and $a \gt 0$ for a repulsive interaction.
Presumably one should be able to construct a canonical ensemble with solutions $\psi(k)$ ($k$ being a momentum label of the above equation), but this is by no means obvious because of the non-linearity represented by the $|\psi(r)|^2$ term. Here a "zero temperature" state would correspond to a perfect BEC with no inhomogeneities, i.e. the vacuum solution of the GPE. However, the entire system is at some finite temperature $T \lt T_c$ as noted above. The resulting thermal fluctuations will manifest in the form of inhomogeneities in the condensate, the exact form of which will be determined by the solutions of the GPE.
Of course, the GPE's regime of validity is that of dilute bose gases ($l_p \gg a_s$ - the average interparticle separation $l_p$ is much greater than the scattering length). For strong coupling I do not know of any similar analytical formalism. If I had to take a wild guess I'd say that the strong-coupling regime could be made analytically tractable by mapping it to a dual gravitational system, but that's another story altogether.
As $l_p$ approaches $a_s$ from above, the GPE breaks down and it will have singular solutions for any given $T$ and these are likely the singularities that you are referring to.
Reference: The single best reference I can suggest is Fetter and Walecka's book on many-body physics. I'm sure you can find more compact sources with a little effort. But generally the brief explanations leave one wanting for a comprehensive approach such as the one F&W provides.
This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivs 3.0 Unported License.
To determine the upper limit on chemical potential for a gas of $\mathcal N$ bosons, look at the form of the Bose distribution in the grand canonical ensemble with $\langle N \rangle = \mathcal N$. When using the GCE, it's easiest to work at chemical potential $\mu$ and to then choose $\mu(\mathcal N)$ so that $\langle N\rangle(\mu)=\mathcal N$. Each state $s$ has average occupancy
$$
\langle n_s\rangle=\frac{\sum_{n\geq 0} ne^{-\beta n(\epsilon_s-\mu)}}{\sum_{n\geq 0} e^{-\beta n(\epsilon_s-\mu)}}=\frac{1}{\Xi_s}\frac{\partial}{\partial(\beta\mu)}\Xi_s,\quad \Xi_s=\frac{1}{1-e^{-\beta(\epsilon_s-\mu)}},\\
=-\partial_{(\beta\mu)}\log(1-e^{-\beta\epsilon_s+(\beta\mu)})=\frac{e^{\beta\mu}}{1-e^{-\beta(\epsilon_s-\mu)}}.
$$
This is finite as long as $\mu<\epsilon_s$. In order for $\langle N\rangle=\sum_s \langle n_s\rangle$ to be finite we need $\mu<\min_s \epsilon_s=\epsilon_0$. Hence, for any system of bosons where $N$ is conserved we have $\mu<\epsilon_0$. It is conventional to set $\epsilon_0=0$ for simplicity, but you can have systems with $\epsilon_0\neq 0$. As you implied by your final question, in these systems the critical value of $\mu$ is $\epsilon_0$ in the thermodynamic limit, with $N, V, E\rightarrow \infty$ and $\mu, p, T$ held constant. Of course, if the system has no BEC phase then as $T\rightarrow 0$, the chemical potential $\mu$ never exceeds some value $\mu_\max<\epsilon_0$.
Best Answer
It feels like you are going to fast in your way to think. The difference between the chemical potential and the ground state energy is clear :
The ground-state energy of your system is here $\epsilon_\textbf{k}=0$ correponding to the ground-state $|\textbf{k=0}\rangle$, which is macroscopically occupied in a BEC (i.e. $N\sim N_0$).
The chemical potential $\mu$ is the increment in energy that you give to your condensate by adding one particle. It is generally defined using free energy $F$ : $$\mu=\left(\frac{\partial F}{\partial N}\right)_{V,\;T}$$
If you want to evaluate $\mu$, considere the ground-state occupation : $$N_0=\frac{1}{e^{-\beta\mu}-1}$$
If you want $N_0$ to be the biggest possible (macroscopic occupation), you must fulfill the condition $e^{-\beta\mu}\rightarrow 1$, in other words $\mu\rightarrow 0$.
For larger resolution on $\mu$, one can do a Taylor expansion on the expression of $N_0$ : $$N_0\underset{\mu\rightarrow 0}{\sim}-\frac{1}{\beta\mu}$$ Then, $$\mu\simeq -\frac{1}{N_0\beta}=-\frac{k_\mathrm{B} T}{N_0}=-\frac{k_\mathrm{B} T}{N-N_\mathrm{exc}}$$