You could define the data generating process as follows:
\begin{equation}
x_{1,t} = \phi_{11} x_{1,t-1} + \cdots \phi_{1p} x_{1,t-p} + \epsilon_{1,t} \,, \quad \hbox{NID}(0, \sigma^2_1) \\
x_{2,t} = \phi_{21} x_{2,t-1} + \cdots \phi_{2p} x_{2,t-p} + \epsilon_{2,t} \,, \quad \hbox{NID}(0, \sigma^2_2) \\
\hbox{Cov}(\epsilon_{1,t}, \epsilon_{2,s}) = \sigma \, \hbox{ if } t=s \hbox{ and zero otherwise.}
\end{equation}
The R functions MASS::mvrnorm or mvtnorm::rmvnorm can be used to generate draws from the multivariate Gaussian distribution.
The code below generates several series from the model defined above with $p=1$ and stores in object rhos
the correlation between each pair of series. The covariance matrix of the disturbance terms is defined with values $\sigma^2_1=2$, $\sigma^2_2=3$ and $\sigma=0.8$.
require("mvtnorm")
set.seed(123)
sigma <- diag(c(2,3))
sigma[1,2] <- sigma[2,1] <- 0.8
n <- 200
niter <- 1000
rhos <- rep(NA, niter)
for (i in seq_len(niter))
{
eps <- mvtnorm::rmvnorm(n = n, mean = rep(0, 2), sigma = sigma)
x1 <- arima.sim(n = n, model = list(ar = c(0.7)), innov = eps[,1])
x2 <- arima.sim(n = n, model = list(ar = c(-0.4)), innov = eps[,2])
rhos[i] <- cor(x1, x2)
}
summary(rhos)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 0.00732 0.13480 0.16990 0.16690 0.19990 0.31280
It can be checked that the correlation between $x_{1,t}$ and $x_{2,t}$ when both series follow an AR(1) process is given by:
\begin{equation}
\rho = \frac{\sigma/(1 - \phi_{11}\phi_{21})}{\sqrt{\frac{\sigma^2_1}{1-\phi_{11}^2}}\sqrt{\frac{\sigma^2_2}{1-\phi_{21}^2}}}
\end{equation}
For the parameter values of the example we have:
sd1 <- sqrt(2 / (1 - 0.7^2))
sd2 <- sqrt(3 / (1 - 0.4^2))
(0.8/(1 + (0.7*0.4))) / (sd1 * sd2)
# [1] 0.1670049
which is in agreement with the average value obtained for the correlation in the small simulation exercise.