Let's assume I wanted to sample from a $\mathcal N(\boldsymbol 0, \boldsymbol\Sigma)$ distribution, where $$\boldsymbol\Sigma = \begin{pmatrix} 1 & \rho \\ \rho & 1\end{pmatrix}$$ for some $\rho\in(-1,1)\setminus\{0\}$. Let $\boldsymbol X\in\mathbb R^{n\times 2}$ denote the sample (of size $n$) when sampling from this distribution. Now I apply the rank statistic to $\boldsymbol X$ columnwise to transform each column to a sequence of (unique) integers ranging from $1$ to $n$. Let $\boldsymbol Y$ denote this transformed matrix. Since $\rho\neq 0$, there is a dependence (depending on $\rho$) between the columns of $\boldsymbol X$. After the transformation there should be also a correlation between the columns of $\boldsymbol Y$ (as large values in the first column correspond to large positive/negative values in the second column due to the correlation). Let $\varrho$ denote this correlation.
Is it possible to sample $\boldsymbol Y$ "directly" from a discerte bivariate distribution with correlation $\varrho$? In other words, is it possible to sample from a bivariate discrete uniform distribution with a certain correlation structure described by $\varrho$? If so, how? Unfortunately, I have no idea how to appraoch this. I came across copulas, which seemed to be promising. However, sampling discrete values from a copula seems to be tedious...
Example: $$\boldsymbol X = \begin{pmatrix} 1.452061 &-0.7210985 \\ -0.753706 & -0.4739298 \\ 1.146363 & 0.6793950 \end{pmatrix} \qquad\stackrel{\text{after the transformation}}{\leadsto}\qquad\boldsymbol Y = \begin{pmatrix} 3 & 1 \\ 1 & 2 \\ 2 & 3\end{pmatrix}$$