erf_arg = -((param['rho_star'] - rho_bar +
                 (rho_bar - rho_0) * np.exp(-n / (tau_eff * f))) /
                (np.sqrt(sigma_rho_sq * (1.0 -
                                         np.exp(-2.0 * n / (tau_eff * f))))))