Using a copula to generate correlated data draws; attenuation bias in result correlation

Question

I would like to generate correlated data draws across two unrelated distributions (or across two empirical datasets that are from unknown distributions). The correct way to do this is to use a copula. I'm not a copula expert.

Mechanically, my steps are as follows:

1. Do multivariate normal draws with a fixed correlated structure.
2. Map the resulting draws to their quantiles.
3. Use these quantiles on the CDF (or eCDF) of the desired target distribution
4. Take a rank correlation of the resulting data to verify the correlation structure is preserved. (Simulate this many times to ensure that finite sample bias is zero)

I wrote a bunch of custom code to do this for my use case, only to discover that I had a slight (3-5%) attenuation bias. Suppose I was drawing for 2 variables with correlation rho=0.6, I'd get data correlated at 0.58. Suppose rho=0.2, I'd get 0.19. Suppose rho=-0.2, I'd get -0.19. The magnitude of the bias is multiplicative and setting rho=1 or rho=0 worked the way you'd expect (1/0 observed correlation)

Figuring there was something wonky with my code, I tried using the copula package in R to do it parametrically:

# Load copula library
library(copula)

# Normal copula, rho=0.6
copula = normalCopula(0.6, dim = 2, dispstr="ex")

# Output distributions: uniform, normal, parameterized
decl = mvdc(copula = copula,
     margins = c("unif", "norm"),
     paramMargins = list(list(min=1, max=1000),
                         list(mean=75, sd=8)))

# Run 1,000 simulations of 1,000 observations
results = replicate(1000, 
          rMvdc(1000, decl))

# Get the average correlation coefficient across all 1,000 sims
summary(t(apply(results, 3, cor, method="spearman")))

The result is 0.58 correlation.

Before writing this question, I did a search on whether this is a known property of copula. This article suggests it is the case -- under "Using Rank Correlation Coefficients":

The correlation parameter, $\rho$, of the underlying bivariate normal determines the dependence between X1 and X2 in this construction. However, the linear correlation of X1 and X2 is not $\rho$. For example, in the original lognormal case, a closed form for that correlation is: $cor(X1, X2) = \frac{e^{\rho\sigma^2} - 1}{e^{\sigma^2} - 1}$ which is strictly less than $\rho$, unless $\rho$ is exactly 1. In more general cases such as the Gamma/t construction, the linear correlation between X1 and X2 is difficult or impossible to express in terms of $\rho$, but simulations show that the same effect happens.

Well, sugar. I'm in the "more general case" and I would like a closed form for the correlation so I can develop a correction to fix this. My end goal is writing a method in a package for simulating particular data in R (think wakefield but a little different). It's not the end of the world if this is a known property that I cannot work around, but I would at least like to understand a little bit better what's going on to be able to document the behaviour.

Can anyone help intuit what's going on here or point me to a discussion of this phenomenon or, even better, let me know what kind of correction I can apply to resolve this given arbitrary marginal distributions and/or empirical distributions?

Thanks in advance.