Low correlation problem in control variate method

I have been trying to use control variate method to reduce the variance of my Monte Carlo Simulation, however, the model is suffering from low correlation problem which makes the control variate method does not work, how do I deal with this.

I am thinking that could I generate two set of "random" number that artifically correlates with each other, but I am not sure what this will do to the expected value of my monte carlo simulation, Thanks for any comment.

option_monte <- function(sigma,s_0,r,k, nSim,tau, type, q){
set.seed(204)

Z <- rnorm(nSim, 0,1)
W_T <- sqrt(tau) *Z
S_T <- s_0*exp((r -q- 0.5*sigma^2)*tau + sigma * W_T)

if (type == "call"){
sim_call_payoff <- exp(-r*tau)*pmax(S_T-k,0)
option_price <- mean(sim_call_payoff)
variance <-   var(sim_call_payoff)/length(sim_call_payoff)
x <- rep(NA,length(sim_call_payoff))
x <- sim_call_payoff
return(c(option_price, variance,x))
}

if (type == "put") {
sim_put_payoff <- exp(-r*tau)*pmax(k-S_T,0)
option_price <-  mean(sim_put_payoff )
variance <-   var(sim_put_payoff)/length(sim_put_payoff)
x <- rep(NA,length(sim_put_payoff))
x <- sim_put_payoff
return(c(option_price, variance,x))
}
}

#control variate

MC_put1 <- option_monte(0.3,45,0.02,40,1000,2,"put", 0.04)
MC_call <- option_monte(0.3,45,0.02,40,1000,2,"call", 0.04)
BS_call <- BS(45,40,0.04,0.3,0.02,2,"call")
BS_put  <- BS_put<- BS(45,40,0.04,0.3,0.02,2,"put")
MC_Put <- MC_put1[1] + (BS_call - MC_call[1])
sd_new1 <- sqrt(MC_put1[2] + MC_call[2] - 2*cov(MC_put1[3:length(MC_put1)],MC_call[3:length(MC_call)])/length(MC_call))

#Result
var_reduced_put_price Monte_Carlo_put BS_analytic_put Var_monte Var_reduced_Var

4.732958        5.092364        5.273973  0.2127453       0.5008263

As you can see, the "reduced" standard error is 0.5 which is larger than unreduced (0.2127). I think this is due to the two independent random variables that I used for my control variates. However, I have no good solution for this.

$$`$$