Sigma <- matrix(c(1, 0.5, 0.5, 1), nrow = 2)
X <- MASS:::mvrnorm(200000, mu = rep(0, 2), Sigma = Sigma, empirical = T)
cov(X)
     [,1] [,2]
[1,]  1.0  0.5
[2,]  0.5  1.0

P <- pnorm(X, lower.tail = FALSE)
cov(P)
           [,1]       [,2]
[1,] 0.08329436 0.04028184
[2,] 0.04028184 0.08335959

qt.scaled <- function (p, df, mean = 0, sd = 1, ncp, lower.tail = F, log.p = FALSE)
{
mean + sd * stats::qt(p, df, ncp = ncp, log.p = log.p)
}

Z1 <- qt.scaled(P[, 1], df = 5, mean = 5, sd = 1)
Z2 <- qt.scaled(P[, 2], df = 10, mean = 5, sd = 2)
Z <- cbind(Z1, Z2)
cov(Z)
        Z1       Z2
Z1 1.678528 1.421597
Z2 1.421597 5.000570