Bose-Einstein condensation is a different phase of matter, meaning that you need to go through a phase transition to reach it. Hence, some thermodynamics quantities (or their derivatives) will exhibit a discontinuity the usual one for the BEC is the heat capacity.
Since it's a phase transition, you want to find the critical value of the parameter which is driving it. <br> This is the temperature, and it is $T_c$.
The physical meaning of the BEC is a saturation effect. <br> For non-interacting bosons, the mean occupancy state $j$ is given by:
$$ f(E_j) = \frac{1}{e^{(E_j - \mu)/kT}-1} .$$
Now you see that the ground state $E_0$, the occupancy is infinity. This is because for the $E_0$ state the chemical potential $\mu$ also needs to be zero, in order to guarantee $f$ to still be positive. Physically, the chemical potential is defined as $\partial U/\partial N$, i.e. the energy added when you add one particle to the system. But if you add it to the $E=0$ state, then the extra energy is 0...
Bose-Einstein condensation begins when you saturate the excited states and start macroscopically occupying the ground state, which has infinite occupancy.
Below $T_c$, $f(E_0)$ starts blowing up so it does not make sense using the above distribution anymore, since the atoms start amassing into the ground state. <br>
So $T_c$ is extracted from when your total $N$ is equal to the number of atoms in the excited states, $N_{ex} = \int_0 ^{\infty}dE \, g(E) \, f(E)$ where $g(E)$ is the *density of states*, i.e. the number of states in a given interval $[E, E+dE]$. The sum should have been over the *states*, but I changed it to the *energy* $E$ just by introducing this density of states term.
The density of states $g(E)$ scales with the number of dimensions $d$. For a free $d$ dimensional system it goes as $g(E) \propto E^{d/2 -1}$, while for $d$ dimensional harmonic potential it scales as $g(E) \propto E^{d-1}$.
In general you can write:
$$ g(E) \propto E^{\alpha -1},$$
with $\alpha$ being the number of degrees of freedom in the system divided by 2. For free particles in $d$ dimensions, $\alpha = d/2$, and for a $d$ dimensional harmonic potential, the degrees of freedom are $2d$ ($d$ translations and $d$ oscillations) so $\alpha = d$. All agree with the above.
The integral above can be rewritten as:
$$N_{ex} = \int_0 ^{\infty}dE \, g(E) \, f(E) \propto (k T_c)^{\alpha} \int_0 ^{\infty} dx\frac{x^{\alpha-1}}{e^x - 1} $$
where I defined $x$ as $E/k T_c$.
The intregral
$$ \int_0 ^{\infty} dx \frac{x^{\alpha-1}}{e^x - 1} = \Gamma(\alpha) \zeta(\alpha), \qquad \alpha > 1$$
with $\Gamma$ being the gamma function, $\zeta$ being the Riemann zeta function.
Which gives you:
$$ k T_c \propto \frac{1}{[\Gamma(\alpha) \zeta(\alpha)]^{1/\alpha}}. $$
To have a BEC transiton, you want $T_c \neq 0$, i.e. a non-trivial solution.
In free space, $d=2,3,4$ have $\alpha = 1, 3/2, 2$:
$$
\begin{array}{ccc}
\alpha & \Gamma(\alpha) & \zeta(\alpha) \\
\hline
1/2 & \text{integral does not converge} \\
1 & 1 & \infty \\
3/2 & \sqrt{\pi}/2 & 2.612 \\
2 & 1 & \pi^2/6 \\
\dots & \dots & \dots
\end{array} $$
So in a free system with $d = 1,2$ the only solution is $T_c = 0$, but for $d>2$, $T_c$ is finite.
Hence, in 2D you may well have a lot of bosons in the ground state. But they are not forced to be there because of saturation, therefore it is not a distinct phase of matter. And the heat capacity (or whatever) will not have discontinuities.