  In this paper, we introduce and discuss an exactly energy-conserving
Particle-in-Cell method for arbitrary curvilinear coordinates. The flexibility
provided by curvilinear coordinates enables the study of plasmas in
complex-shaped domains by aligning the grid to the given geometry, or by
focusing grid resolution on regions of interest without overresolving the
surrounding, potentially uninteresting domain. We have achieved this through
the introduction of the metric tensor, the Jacobian matrix, and contravariant
operators combined with an energy-conserving fully implicit solver. We
demonstrate the method's capabilities using a Python implementation to study
several one- and two-dimensional test cases: the electrostatic two-stream
instability, the electromagnetic Weibel instability, and the geomagnetic
environment modeling (GEM) reconnection challenge. The test results confirm the
capability of our new method to reproduce theoretical expectations (e.g.
instability growth rates) and the corresponding results obtained with a
Cartesian uniform grid when using curvilinear grids. Simultaneously, we show
that the method conserves energy to machine precision in all cases.
