  In this paper, we compute the constrained QCD effective potential up to
two-loop order with finite quark mass and chemical potential. We present the
explicit calculations by using the double line notation and analytical
expressions for massless quarks are obtained in terms of the Bernoulli
polynomials or Polyakov loops. Our results explicitly show that the constrained
QCD effective potential is independent on the gauge fixing parameter. In
addition, as compared to the massless case, the constrained QCD effective
potential with massive quarks develops a completely new term which is only
absent when the background field vanishes. Furthermore, we discuss the relation
between the one- and two-loop constrained effective potential. The surprisingly
simple proportionality that exists in the pure gauge theories, however, is in
general no longer true when fermions are taken into account. On the other hand,
for high baryon density $\mu_B$ and low temperature $T$, in the massless limit,
we do also find a similar proportionality between the one- and two-loop
fermionic contributions in the constrained effective potential up to ${\cal
O}(T/\mu_B)$.
