Abstract: Understanding the internal structure and dynamics of protons and neutrons, the complex many-body systems consisting of strongly-interacting quarks and gluons is at the core of exploring the visible matter universe. Gluons, which serve as mediator bosons of the strong interaction, play a key role in the nucleon’s mass and spin structures. In contrast, understanding of the gluon distributions and their role in hadron structures remains one of the most challenging but fundamental issues in nuclear and particle physics. In this talk, I will present lattice QCD calculations of matrix elements of bilocal operators composed of two gluon fields that can be used to determine the unpolarized and polarized gluon distributions. In particular, I will focus on the first lattice QCD determination of the gluon helicity parton distribution function with numerical evidence toward disfavoring negative gluon polarization in the nucleon.