This question boils down to: what is the likelihood of an observed ordering of a vector of independent normal variates with potentially different means?
It turns out that this question has a somewhat tractable answer based on the multivariate normal CDF; see here combinatorics - Compute probability of a particular ordering of normal random variables - Mathematics Stack Exchange
However, Stan does not have a multivariate normal CDF function, and in general it’s a hard one to compute. Some progress in this direction in Stan has been made e.g. here: Multivariate normal CDF and here Multivariate normal cdf (those posts have the same title but are different posts). Perhaps @spinkney or @martinmodrak has more to say about this?