Improving control variate for variance reduction

I have tried stock price as control variate for my monte carlo simulation, and I am trying to reduce the variance of my estimated price for European Put option. And the code look like this:

path <- function(rho,T, S0, r, kappa, theta,xi,v0, Maturity ) { 
  
  dt = Maturity / T
  COV <- matrix(c(1,rho,rho,1), nrow = 2,ncol = 2)
  W <- rmvnorm(T, , COV)
  W_v <- W[,1]
  W_s <- W[,2]
  
  Vt <- rep(0,T)
  Vt[1] <- v0
  St <- rep(0,T)
  St[1] <- S0
  
  for (i in 2:T){
    Vt[i] <- Vt[i-1] + kappa* (theta - Vt[i-1])*dt + xi* sqrt(Vt[i-1]) * W_v[i-1] * sqrt(dt)
    St[i] <- St[i-1]*exp((r-0.5*Vt[i-1]) *dt + sqrt(dt*Vt[i-1]) * W_s[i-1])
    
  }
  St
}
set.seed(214)
N<-2500
Stm <- replicate(100000, path(0.5,N,45,-0.02,10,0.09,0.2,0.09,2))
payoff <- pmax (40 - Stm[N,], 0 ) *exp(-0.02 * 2)
Hest_MC1 <- mean(payoff)
S_aver <- mean(Stm[N,])

ST <- 45* exp(-0.02 *2 )

hest_mc <- Hest_MC1 + -reg$coefficients[[2]]*(S_aver - ST)
hest_mc
reg <- lm(payoff~Stm[N,])

reg$coefficients


sd_old <- sqrt(var(payoff)/length(payoff))
sd_reduced <- sqrt(var(payoff)/length(payoff) + reg$coefficients[[2]]^2 *var(Stm[N,])/length(Stm[N,])
                         + 2*( -reg$coefficients[[2]] )*cov(Stm[N,],payoff)/length(payoff))
# result
        sd_old sd_reduced
[1,] 0.06926151 0.04898853

As you can see, the result only improved marginally after 100,000 Simulation. Do anyone has any Idea how to continue improve on this? I did try the European Put (plain) as control variate, however, the payoff is very weakly correlated with the payoff of the option I tried to reduce variance on.

Any comment would be appreciated.