The energy-energy correlation (EEC) function in $e^+e^-$ annihilation is currently the only QCD event shape observable for which we know the full analytic result at the next-to-leading order (NLO). In this work we calculate the EEC observable for gluon initiated Higgs decay analytically at NLO in the Higgs Effective Field Theory (HEFT) framework and provide the full results expressed in terms of classical polylogarithms, including the asymptotic behavior in the collinear and back-to-back limits. This observable can be, in principle, measured at the future $e^+e^-$ colliders such as CEPC, ILC, FCC-ee or CLIC. It provides an interesting opportunity to simultaneously probe our understanding of the strong and Higgs sectors and can be used for the determinations of the strong coupling.