As you know, simulating AR(1) is to simulate the distributed error path.
Assume the bivariate errors distributed $\sim F(x),\sim F(y)$ with copula $C(u,v)$ to model their dependence.
Then the bivariate joint error distribution is given by Sklar's theorem:
$$F(x,y)=C(F(x),F(y))$$
You can simulate from this distribution using Conditional Sampling:
To obtain a realization of a bivariate Copula $C(u,v)$, one draws the first variable $u$ as
random number $\sim U(0,1)$. The second variable $v$ is generated from another
independent random number $z$ plugged into the inverse Copula $C^{-1}(z\,|u=u)$ under the first generated
(conditional) random number $u$:
- Draw $\bar{u},\bar{z}\sim U(0,1)$
- Set $\bar{v} = C_{\bar{u}}^{-1}(\bar{z})$ (quasi-inverse Copula under $\bar{u}$, or conditional $C^{-1}(t,u\,|u=\bar{u})$ )
From this you get $$(x=F(\bar{u})^{-1},y=F(\bar{v})^{-1})$$ as your two simulated errors for the AR(1) process.