  We discuss a gauge-invariant prescription to take the mean-field
approximation self-consistently in the PNJL model (Nambu--Jona-Lasinio model
with the Polyakov loop). We first address the problem of non-vanishing color
density in normal quark matter, which is an artifact arising from gauge-fixed
treatment of the Polyakov loop mean-fields. We then confirm that the gauge
average incorporated in our prescription resolves this problem and ensures
color neutrality. We point out that the proposed method has an advantage in
computing the expectation value of any function of the Polyakov loop matrix. We
discuss the Casimir scaling as an immediate application of the method.
