[en] A new approach for implicit PIC plasma simulation is described. In this method, only the deviation of the true distribution f from a local Maxwellian f_M(x,t) is represented by PIC methods. The parameters of the Maxwellian (density, mean velocity, temperature) satisfy the usual moment equations, which are coupled to Maxwell's equations and differenced implicitly. To close the moment equations, the traceless stress and heat flux are evaluated on the mesh from δf, the difference between f and f_M by advancing a collection of simulation particles and using PIC methods. A 1-D electromagnetic, two-fluid code based on this algorithm has been constructed and has demonstrated several features of the method. First, noise is greatly reduced for situations in which the plasma remains close to LTE. Second, linear numerical stability seems insensitive to methods used to evaluate the PIC related moments. Thus, for example, substepping, orbit averaging, and gyroaveraging have been demonstrated to give stable time advance schemes. Finally, nonlinear stability depends on maintaining consistency between the particle and moment data. These features are illustrated by examples of the simulation of linear waves. Extension of this method to 3-D plasma simulation is described