Nonlocal gluon condensates are vacuum expectations of the product of gluon field strength tensors. Short-distance expansions of two-, three-, and four-gluon condensates are presented up to dimension-8 local operators. We propose a method for calculating the Wilson coefficients based on the presented expansions and the Feynman diagram technique in the background field approach. The method is demonstrated using the glueball current correlators as examples. Methodological aspects of the background field approach are discussed in relation to glueball studies within QCD sum rules. We confirm the results for Operator Product Expansion (OPE) of the two-gluon 0 ±+ glueball current correlators and calculate additional contributions coming from dimension-6 four-quark condensate and dimension-8 mixed quark-gluon condensates. The OPEs used in QCD sum rules for three-gluon 0 ±+ glueballs are revisited up to dimension-6 order.product 0| : G µν (x)[x; y]G αβ (y) : |0 , where G µν is the gluon field strength tensor and [x; y] is a Wilson line, see details in Sec. 2.1.The methodological aspects of QCD SR are often considered for the case of quark current correlators, see [38], while the glueball currents correlators have particularities left unattended. Therefore, applications of the proposed approach are given by the examples of the glueball currents correlators used in QCD SRs studies [39][40][41][42].Glueballs are an exotic state that contains only gluons and no valence quarks. This type of state has candidates among observations [43,44] and is included in the running and projected large-scale experiments: Belle (Japan), BESIII (Beijing, China), LHC (CERN), GlueX (JLAB,USA), NICA (Dubna, Russia), HIAF (China), and FAIR (GSI, Germany). Electorn-Ion-Colliders have potential for glueball state observations [45]. The reviews of glueball physics can be found in [43,46]. There are many applications of QCD SR to glueball states [7,[39][40][41][42].Here we shortly discuss some of them. The first QCD SR study [39] of glueballs considered a pseudo-scalar 0 −+ state with an obtained mass of ∼ 1 GeV, where two-gluon current was used. Later QCD SR was applied [40] to a scalar 0 ++ glueball state, where the glueball mass was estimated to be ∼ 0.7 GeV. The OPE used in SRs for the scalar and pseudoscalar glueballs was extended [61][62][63]75] by including the direct instanton contribution, the two loop radiative corrections [76] to the perturbative term, and the one loop radiative correction [77] to dimension-4 term (dimension of the considered condensate is four in energy units). Three-gluon glueballs were first studied in [41] for a 0 ++ -state. Then, QCD SRs for three-gluon glueballs was extended [42] to the 0 −+ scalar, vector and tensor states [58].Here we study the technical aspects of the calculation of OPE's coefficients up to dimension-8 order and suggest a new way to organize and perform calculations using NLCs [35][36][37]. Using the algorithms formulated in [34] for calculating higher power corrections, we extend the results of [3...