The single renormalon-chain contribution to the correlator of scalar currents in QCD is calculated in the MS-scheme in the limit of a large number of fermions, N f . At n-loop order we find that in the MS-scheme the factorial growth of the perturbative coefficients due to renormalons takes over almost immediately in the euclidean region. The essential differences between the large-order growth of perturbative coefficients in the present scalar case, and in the previously-studied vector case are analysed. In the timelike region a stabilization of the corresponding perturbative series for the imaginary part, with n-loop behaviour S n /[log(s/Λ 2 )] n−1 , where S n is essentially constant for n ≤ 6, is observed. Only for n ≥ 7 does one discern the factorial growth and alternations of sign. We use the new all-orders results to scrutinize the performance of multiloop estimates , using a largeβ 0 = (11N c − 2N f )/12 approximation, the so-called "naive nonabelianization" procedure, and within the effective charges approach. The asymptotic behaviour of perturbative coefficients, in both large-N f and large-N c limits, is analysed both in the spacelike and timelike regions. A contour-improved resummation technique in the time-like region is developed. Some subtleties connected with scheme-dependence are analysed , and illustrated using results in the MS and V -schemes. The all-orders series under investigation are summed up with the help of the Borel resummation method. The results obtained are relevant to the analysis of the theoretical uncertainties in the 4-loop extractions of the running and invariant s-quark masses from QCD sum rules, and in calculations of the Higgs boson decay width into a quark-antiquark pair.