  In QCD with $N_c$ colors the anti-symmetric valence quark color space singlet
state $\sim \epsilon_{i_1\cdots i_{N_c}} |i_1,\cdots, i_{N_c}\rangle$ of the
proton corresponds to the reduced density matrix $\rho_{ij}=(1/N_c)
\delta_{ij}$ for a single color degree of freedom. Its degenerate spectrum of
eigenvalues, $\lambda_i=1/N_c$, the purity $\mathrm{tr}~\rho^2 = 1/N_c$, and
the von~Neumann entropy $S_\mathrm{vN}=\log(N_c)$ all indicate maximal
entanglement of color. On the other hand, for $N_c\to\infty$ the spatial wave
function of the proton factorizes into valence quark wave functions determined
by a mean field (E. Witten, Nucl. Phys. B 160 (1979) p. 57) where there is no
entanglement of spatial degrees of freedom. A model calculation at $N_c=3$
using a simple three quark model light-front wave function by Brodsky and
Schlumpf, predicts percent level entanglement of spatial degrees of freedom.
Using light-cone perturbation theory we also derive the density matrix
associated with the four parton $|qqqg\rangle$ Fock state. Tracing out the
quarks, we construct the reduced density matrix for the degrees of freedom of
the gluon, which encodes its entanglement with the sources. Our expressions
provide the dependence of the density matrix on the soft cutoff $x$ for the
gluon light-cone momentum, and on the collinear and ultraviolet regulators.
Numerical results obtained in a simple approximation indicate stronger
entanglement for the gluon (with $x_g < \langle x_q\rangle$) than for quarks in
the three quark Fock state.
