A recently developed implicit method for solving the set of coupled particle and field equations arising in particle-in-cell plasma simulation is described in detail. This implicit integration scheme is motivated by the desire to study efficiently low-frequency, long-wavelength plasma phenomena using a large time step. In particular, this method allows the use of a time step which is larger than the electron plasma period, when electron plasma oscillations are not of interest, and provides selective damping of the distorted remnant of the electron plasma oscillation. The implicit scheme presented here uses particle data directly without introducing fluid moment equations as an intermediary between the field and particle equations. In an electrostatic model, the essence of our scheme is a linearization of the charge density at the advanced time about an explicit approximate density and the computation of the incremental correction to the charge density that is linear in the advanced field. We are led to an elliptic field equation whose coefficients depend directly on particle data accumulated on the spatial grid in the form of an effective linear susceptibility. Prediction and iterative refinement of the solution of the implicit equations, and spatial difference representations of the equations are given. Residual restrictions on time step are described. It is demonstrated that convergence is superior when spatial diferencing and filtering are done in a consistent manner.