  The evolution of Bose-Einstein condensates is amply described by the
time-dependent Gross-Pitaevskii mean-field theory which assumes all bosons to
reside in a single time-dependent one-particle state throughout the propagation
process. In this work, we go beyond mean-field and develop an essentially-exact
many-body theory for the propagation of the time-dependent Schr\"odinger
equation of $N$ interacting identical bosons. In our theory, the time-dependent
many-boson wavefunction is written as a sum of permanents assembled from
orthogonal one-particle functions, or orbitals, where {\it both} the expansion
coefficients {\it and} the permanents (orbitals) themselves are {\it
time-dependent} and fully determined according to a standard time-dependent
variational principle. By employing either the usual Lagrangian formulation or
the Dirac-Frenkel variational principle we arrive at two sets of coupled
equations-of-motion, one for the orbitals and one for the expansion
coefficients. The first set comprises of first-order differential equations in
time and non-linear integro-differential equations in position space, whereas
the second set consists of first-order differential equations with
time-dependent coefficients. We call our theory multi-configurational
time-dependent Hartree for bosons, or MCTDHB($M$), where $M$ specifies the
number of time-dependent orbitals used to construct the permanents. Numerical
implementation of the theory is reported and illustrative numerical examples of
many-body dynamics of trapped Bose-Einstein condensates are provided and
discussed.
