Imaginary time propagation , Gross-Pitaevskii equation When deriving the Gross-Pitaevskii equation we already assume that all the atoms are in the ground state. 
Based on this assumption, we minimize the energy functional and get the Gross-Pitaevskii equation(GPE).
GPE gives us the density of the condensate.
What I don't understand is that in imaginary time propagation method , we write the wave function of the condensate(which gives us density) as a superposition of eigen functions corresponding to different energy eigen values.
It is said that when we propagate in imaginary time we get the ground state of the condensate.
But GPE is already a equation for the ground state wave function since all the atoms are in the ground state. So what does it mean to find the ground state of the BEC using imaginary time method?
According to my understanding , all solutions should have same energies since all the atoms are in the ground state.
 A: The non-interacting GPE is just the Schrödinger equation. Is the Schrödinger equation only satisfied by ground states? No.
The energy functional (or action, in general) minimisation actually does not care about the eigenenergies of the solutions to the GPE. For the energy functional (action) to be minimised such that $\delta E = 0$ ($\delta S=0$), you need the solution $\phi$ to follow a certain path.  That path is given by the equation of motion.
Given an action $S = \int L(\mathbf{q}, \dot{\mathbf{q}})\mathrm{d}t$, you can derive the Euler-Lagrande equations that will result in the equation of motion for the field(s):
$$\delta S = 0 \quad \Rightarrow \quad \frac{\partial L}{\partial \mathbf{q}}-\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial L}{\partial \mathbf{\dot{q}}} =0, $$
so any ($\mathbf{q}, \dot{\mathbf{q}}$) solution that satisfies the Euler-Lagrange equation (no matter its energy) will give that the RHS = 0 and hence $\delta S= 0$.
The reasoning is essentially the same for the energy functional and the GPE.
If your solution $\phi$ satisfies the GPE, then it minimises the energy functional $\delta E$. However, there might be several $\phi_i$ that satisfy the GPE, and they all give you a stationary functional. It's the same as the solution to $\partial_x (\sin(x)) = 0$: it gives many solutions, each at increasing $x$, but they all give you a minimum of the function.
Imaginary time evolution then finds the solution to the specific GPE (i.e. with a specified potential $V$ and interaction energy $U_0$) with the lowest energy eigenvalue. Provided your initial guess has a finite overlap with the true ground state.
