  The exact ground state of the reduced BCS Hamiltonian is investigated
numerically for large system sizes and compared with the BCS ansatz. A
"canonical'' order parameter is found to be equal to the largest eigenvalue of
Yang's reduced density matrix in the thermodynamic limit. Moreover, the
limiting values of the exact analysis agree with those obtained for the BCS
ground state. Exact results for the ground state energy, level occupations and
a pseudospin-pseudospin correlation function are also found to converge to the
BCS values already for relatively small system sizes. However, discrepancies
persist for a pair-pair correlation function, for inter-level correlations of
occupancies and for the fidelity susceptibility, even for large system sizes
where these quantities have visibly converged to well-defined limits. Our
results indicate that there exist non-perturbative corrections to the BCS
predictions in the thermodynamic limit.
