Initial rotational quantum state of a molecule in an ensemble with Boltzmann distribution

Question

The dynamic of rotational states of an ensemble of molecules in temperature T is the purpose of my simulations. Thermal distribution of molecules is defined by Boltzmann distribution and spacing between the rotational states is small enough that for a temperature like 300K the states from J=0 to J=15 are occupied. A laser pulse interacts with the molecules and this brings them to a series of what is called as "Rotational revivals". Hence after time "t" the state of molecule is:

$$\big|\psi(t)\big\rangle = \sum_{JM} a_{JM}(t)\big|J,M\big\rangle$$
where "J" is the rotational quantum number, M is magnetic quantum number and $|J\rangle$ is the orthogonal kets of rotation space.

I am trying to numerically simulate the process so I need the initial states and My question is: What are the initial expansion coefficient, $a_J(t=0)$, for each J?. I guess that should be $a_J(t=0)=\sqrt{P_J}$ where $P_J$ is the probability of finding the molecule in state J according to Boltzmann distribution. This also holds the normalization condition that $\sum_J |a_J(t=0)|^2=1.$

I am aware of the density matrix but in many papers dealing with the alignment of molecules no density matrix is used to formulate the dynamics of molecules while interacting with the laser pulse.