Calculating mean of correlation between correlated random variables

Question

I am constructing correlated random variables according to the equations given in in this post: How does the formula for generating correlated random variables work?

I understand this equation, although what I don't understand is if I run a bunch of simulations where I calculate the correlation numerically between $Y$ and $X_1$ and calculate the mean of this correlation, it seems to be systematically reduced as compared to the true value $\rho$. It seems, specifically, that there needs to be a "correction factor" to the mean of the correlation of $(N+1)/N$, where $N$ is the number of sampled values. This was just found as a guess, it could also be wrong but seems to fit. Here's the python code to reproduce this issue:

import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import ttest_ind


N = 50000
Npoints_list = [2, 100]

corr_values = [0.5]


list_of_ab = []
list_of_cd = []


rand_b_list = [1.0]
std_scale_factor = 1
C = 0
delta = 0

for Npoints in Npoints_list:
    corr_values_ac = []
    for Delta in rand_b_list:
        list_of_ab = []
        list_of_cd = []
        for corr in corr_values:
            p_values_ab = []
            p_values_cd = []
            for ii in range(N):
                a = (np.random.randn((Npoints))*std_scale_factor+C)

                c =corr*a+np.sqrt(1-corr**2)*(C*(1-corr)/np.sqrt(1-corr**2)+np.random.randn((Npoints))*std_scale_factor)
    
                corr_ac = np.corrcoef(a,c)
        
                corr_values_ac.append(corr_ac[1,0])
                
    print(f'mean of corr of ac for {Npoints} points is {np.mean(corr_values_ac)}')
    print(f'mean of corr of ac for {Npoints} points after correction is {np.mean(corr_values_ac)*(Npoints+1)/Npoints}')