Matrix expression for $\operatorname{vec}(X^{\top}X)$?

I believe $B$ takes the form $$ B = \left(I_n \otimes \left(\operatorname{vec}\left(I_m\right)\right)^{\top} \otimes I_n\right)\left( I_{mn} \otimes K \right), $$ where $K$ is the 'commutation matrix.' That is, $K$ is the matrix such that $K \operatorname{vec}(X) = \operatorname{vec}\left(X^{\top}\right).$

Letting $$ V = \left(I_{mn} \otimes K \right) \operatorname{vec}\left(\operatorname{vec}\left(X\right)\left(\operatorname{vec}\left(X\right)\right)^{\top}\right) = \operatorname{vec}\left(\operatorname{vec}\left(X^{\top}\right)\left(\operatorname{vec}\left(X\right)\right)^{\top}\right), $$ then the $i,j$ element of $\operatorname{vec}\left(X^{\top}X\right)$ should be: $$ \operatorname{vec}\left(X^{\top}X\right)_{i,j} = \sum_k X_{k,i} X_{k,j} = V_{i+mk+mn(k+mj)}, $$ where we index from zero, as suggested by @Chrystomath. The term $\left(I_n \otimes \left(\operatorname{vec}\left(I_m\right)\right)^{\top} \otimes I_n\right)$ then captures that indexing and the repeated $k$.

I have confirmed this relationship with some R code:

library(matrixcalc)
m <- 5
n <- 3
set.seed(1234)
X <- matrix(rnorm(m*n),nrow=m)
LHS <- vec(t(X) %*% X)
VM  <- vec(vec(X) %*% t(vec(X)))
B   <- (diag(n) %x% t(vec(diag(m))) %x% diag(n)) %*% (diag(m*n) %x% commutation.matrix(m,n))
RHS <- B %*% VM
max(abs(RHS - LHS))

[1] 0