If the columns of the matrix are given by
$$U = \left[\matrix{u_1&u_2\ldots&u_n}\right]$$
Then the columns of the iterated Khatri-Rao product are given
$$U^{\odot M}
= \left[\matrix{u_1^{\otimes M}&u_2^{\otimes M}\ldots&u_n^{\otimes M}}\right]$$
where the iterated Kronecker product of a vector with itself has been denoted by
$$a^{\otimes M}
= \operatorname*{\bigotimes}\limits_{k=1}^M a
\;=\; a\otimes a\otimes\ldots\otimes a
$$
If you like you can use the standard basis vectors $\{e_k\}$ to rearrange things.
$$\eqalign{
u_k &= Ue_k \\
u_k^{\otimes M} &= U^{\otimes M}\,e_k^{\otimes M} \\
U^{\odot M} &= U^{\otimes M}
\left[\matrix{e_1^{\otimes M}&e_2^{\otimes M}\ldots&e_n^{\otimes M}}\right] \\
&= U^{\otimes M}\,E_M \\
}$$
This reformulation swaps the iterated Khatri product for an iterated Kronecker product, multiplied by a very large sparse binary matrix.