We start by analyzing the iterative self-consistent solution of the reduced Dyson-like screening equation from which the set of frequency-dependent screened polarizabilities are obtained and subsequently used in the solution of the Casimir-Polder integral in Eq.( ). Before proceeding with the numerical analysis it is first important to understand the adequate number of quadrature points (i.e. frequencies) to be used in the GL quadrature. While a too low number of quadrature points may lead to numerical errors, a too large one, on the other hand, would lead to an unnecessary increase of the computational cost as well as memory requirements and a balance must be therefore found. The analysis reported in the Supporting Information shows that a ten-points GL quadrature leads to errors in the order of 10 -6 and this number of choice will be embraced throughout the discussion.