  In this paper we address the temporal energy growth associated with numerical
approximations of the perfectly matched layer (PML) for Maxwell's equations in
first order form. In the literature, several studies have shown that a
numerical method which is stable in the absence of the PML can become unstable
when the PML is introduced. We demonstrate in this paper that this instability
can be directly related to numerical treatment of boundary conditions in the
PML. First, at the continuous level, we establish the stability of the constant
coefficient initial boundary value problem for the PML. To enable the
construction of stable numerical boundary procedures, we derive energy
estimates for the variable coefficient PML. Second, we develop a high order
accurate and stable numerical approximation for the PML using
summation--by--parts finite difference operators to approximate spatial
derivatives and weak enforcement of boundary conditions using penalties. By
constructing analogous discrete energy estimates we show discrete stability and
convergence of the numerical method. Numerical experiments verify the
theoretical results
