We compute the static self-energy of SU(3) gauge theory in four spacetime dimensions to order α 20 in the strong coupling constant α. We employ lattice regularization to enable a numerical simulation within the framework of stochastic perturbation theory. We find perfect agreement with the factorial growth of high order coefficients predicted by the conjectured renormalon picture based on the operator product expansion.PACS numbers: 11.15. Bt,12.38.Cy,12.38.Bx,11.10.Jj,12.39.Hg Little is known about properties of quantum field theories from first principles. This is particularly so for asymptotically free gauge theories such as quantum gluodynamics. One of the most salient features of this theory is the confinement of charged objects. Yet this property has not been proven, and the best evidence comes from the linearly rising static potential at large distances obtained in lattice simulations. Another expected property is the asymptotic nature of perturbative weak coupling expansions. In four dimensional non-Abelian gauge theories one particular pattern of asymptotic divergence should be determined by the structure of the operator product expansion (OPE). It is usually named renormalon [1] or, more specifically, infrared renormalon. Its existence has also not been proven but only tested assuming the dominance of β 0 -terms, which amounts to an effective Abelianization of the theory, or in the two dimensional O(N ) model [2], where it is suppressed by powers of 1/N . Moreover, the possible non-existence or irrelevance of renormalons in Quantum Chromodynamics has been suggested in several papers, see, e.g. [3,4] and references therein. This has motivated dedicated high order perturbative expansions of the plaquette, e.g. [5][6][7][8], in lattice regularization, with conflicting conclusions. Powers as high as α 20 were achieved in the most recent simulation [9]. However, the expected asymptotic behaviour was not seen. If confirmed, this non-observation would cast doubt on the well-accepted lore of the OPE and renormalon physics (see [10] for a comprehensive review), and would significantly affect the phenomenological analysis of data from high energy physics experiments on the decay of heavy hadrons, heavy quark masses, the running coupling parameter, parton distributions, etc.. Therefore, this issue should be clarified unambiguously.In this letter we present compelling numerical evidence that the expected renormalons indeed exist not only in models but in real gluodynamics. We also argue why previous analyses based on the plaquette have failed to detect them. The vital and new ingredients of our study are as follows. (a) We consider a perturbative series whose leading renormalon is dictated by a dimension d = 1 operator, rather than by the d = 4 plaquette. (b) Using a higher order integrator and employing twisted boundary conditions, among other improvements, we are able to obtain results of unprecedented precision on an extensive set of spacetime volumes. (c) We carefully extrapolate to the infinite volume limi...