By employing the differential equations, we compute analytically the elliptic sectors of two-loop master integrals appearing in the NNLO QCD corrections to CP-even heavy quarkonium exclusive production and decays, which turns out to be the last and toughest part in the relevant calculation. The integrals are found can be expressed as Goncharov polylogarithms and iterative integrals over elliptic functions. The master integrals may be applied to some other NNLO QCD calculations about heavy quarkonium exclusive production, like γ * γ → QQ, e + e − → γ + QQ, and H/Z 0 → γ + QQ, heavy quarkonium exclusive decays, and also the CP-even heavy quarkonium inclusive production and decays.