  For decades, the Vlasov-Darwin model has been recognized to be attractive for
particle-in-cell (PIC) kinetic plasma simulations in non-radiative
electromagnetic regimes, to avoid radiative noise issues and gain computational
efficiency. However, the Darwin model results in an elliptic set of field
equations that renders conventional explicit time integration unconditionally
unstable. Here, we explore a fully implicit PIC algorithm for the Vlasov-Darwin
model in multiple dimensions, which overcomes many difficulties of traditional
semi-implicit Darwin PIC algorithms. The finite-difference scheme for Darwin
field equations and particle equations of motion is space-time-centered,
employing particle sub-cycling and orbit-averaging. The algorithm conserves
total energy, local charge, canonical-momentum in the ignorable direction, and
preserves the Coulomb gauge exactly. An asymptotically well-posed fluid
preconditioner allows efficient use of large time steps and cell sizes, which
are determined by accuracy considerations, not stability, and can be orders of
magnitude larger than required in a standard explicit electromagnetic PIC
simulation. We demonstrate the accuracy and efficiency properties of the
algorithm with various numerical experiments in 2D-3V.
