  We analyse inclusive photoproduction of polarized $^3P_1$ quarkonium in the
framework of QCD. To separate nonperturbative and perturbative parts in the
density matrix of the produced quarkonium we use a method , which is equivalent
to the diagramatic expansion widely used in analysing deeply inelastic
scatterings. A systematic expansion in a small velocity $v$, with which a heavy
quark moves inside the quarkonium in its rest frame, is performed for the
nonperturbative parts, and they are expressed as matrix elements in
nonrelativistic QCD. At the leading order of $v$ there are four matrix elements
representing nonperturbative physics. The perturbative parts are calculated at
the leading order of coupling constants. Some numerical results, especially,
numerical results for HERA are given.
