We perform a comprehensive study of the time evolution of heavy-quarkonium states in an expanding hot QCD medium by implementing effective field theory techniques in the framework of open quantum systems. The formalism incorporates quarkonium production and its subsequent evolution in the fireball including quarkonium dissociation and recombination. We consider a fireball with a local temperature that is much smaller than the inverse size of the quarkonium and much larger than its binding energy. The calculation is performed at an accuracy that is leading-order in the heavy-quark density expansion and next-to-leading order in the multipole expansion. Within this accuracy, for a smooth variation of the temperature and large times, the evolution equation can be written as a Lindblad equation. We solve the Lindblad equation numerically both for a weakly-coupled quark-gluon plasma and a strongly-coupled medium. As an application, we compute the nuclear modification factor for the $\Upsilon(1S)$ and $\Upsilon(2S)$ states. We also consider the case of static quarks, which can be solved analytically. Our study fulfils three essential conditions: it conserves the total number of heavy quarks, it accounts for the non-Abelian nature of QCD and it avoids classical approximations.