Simulate a dataset with effect-size estimates correlated to a given dataset in R

Simulate a dataset with effect-size estimates correlated to a given
dataset in R

I have a dataset which is associated with an outcome with a specific mean
and standard deviation, and I would like to simulate a second dataset that
has betas (effect estimates) which are correlated with the betas from the
first dataset at a specific correlation coefficient rho. I can simulate a
bivariate distribution with correlated errors, but I don't think this will
do the job. What I am looking for is to have a specific correlation
between effect estimates. So for example if I have:
x1 = rnorm(100, 2, 0.5)
y <- rnorm(100, mean=0, sd=1)
beta1= lm(y ~ x1)$coefficients[[2]]
I am trying to find a way to create a beta2 which is correlated with beta1
at a specific rho. Can someone please help me?