  We suggest a framework based on the rainbow approximation with effective
parameters adjusted to lattice data. The analytic structure of the gluon and
ghost propagators of QCD in Landau gauge is analyzed by means of numerical
solutions of the coupled system of truncated Dyson-Schwinger equations. We find
that the gluon and ghost dressing functions are singular in complex Euclidean
space with singularities as isolated pairwise conjugated poles. These poles
hamper solving numerically the Bethe-Salpeter equation for glueballs as bound
states of two interacting dressed gluons. Nevertheless, we argue that, by
knowing the position of the poles and their residues, a reliable algorithm for
numerical solving the Bethe-Salpeter equation can be established.
