  We develop a formalism for computing inclusive production cross sections of
heavy quarkonia based on the nonrelativistic QCD and the potential
nonrelativistic QCD effective field theories. Our formalism applies to strongly
coupled quarkonia, which include excited charmonium and bottomonium states.
Analogously to heavy quarkonium decay processes, we express nonrelativistic QCD
long-distance matrix elements in terms of quarkonium wavefunctions at the
origin and universal gluonic correlators. Our expressions for the long-distance
matrix elements are valid up to corrections of order $1/N_c^2$. These
expressions enhance the predictive power of the nonrelativistic effective field
theory approach to inclusive production processes by reducing the number of
nonperturbative unknowns, and make possible first-principle determinations of
long-distance matrix elements once the gluonic correlators are known. Based on
this formalism, we compute the production cross sections of $P$-wave charmonia
and bottomonia at the LHC, and find good agreement with measurements.
