  We develop the multi-layer multi-configuration time-dependent Hartree method
for bosons (ML-MCTDHB), a variational numerically exact ab-initio method for
studying the quantum dynamics and stationary properties of bosonic systems.
ML-MCTDHB takes advantage of the permutation symmetry of identical bosons,
which allows for investigations of the quantum dynamics from few to many-body
systems. Moreover, the multi-layer feature enables ML-MCTDHB to describe mixed
bosonic systems consisting of arbitrary many species. Multi-dimensional as well
as mixed-dimensional systems can be accurately and efficiently simulated via
the multi-layer expansion scheme. We provide a detailed account of the
underlying theory and the corresponding implementation. We also demonstrate the
superior performance by applying the method to the tunneling dynamics of
bosonic ensembles in a one-dimensional double well potential, where a
single-species bosonic ensemble of various correlation strengths and a weakly
interacting two-species bosonic ensemble are considered.
