In this paper, we develop a provably energy stable discontinuous Galerkin spectral element method (DGSEM) approximation of the perfectly matched layer (PML) for the three and two space dimensional (3D and 2D) linear acoustic wave equations, in first order form, subject to well-posed linear boundary conditions. First, using the well-known complex coordinate stretching, we derive an efficient un-split modal PML for the 3D acoustic wave equation, truncating a cuboidal computational domain. Second, we prove asymptotic stability of the continuous PML by deriving energy estimates in the Laplace space, for the 3D PML in a heterogeneous acoustic medium, assuming piece-wise constant PML damping. Third, we develop a DGSEM for the wave equation using physically motivated numerical flux, with penalty weights, which are compatible with all well-posed, internal and external, boundary conditions. When the PML damping vanishes, by construction, our choice of penalty parameters yield an upwind scheme and a discrete energy estimate analogous to the continuous energy estimate. Fourth, to ensure numerical stability of the discretization when PML damping is present, it is necessary to systematically extend the numerical numerical fluxes, and the inter-element and boundary procedures, to the PML auxiliary differential equations. This is critical for deriving discrete energy estimates analogous to the continuous energy estimates. Finally, we propose a procedure to compute PML damping coefficients such that the PML error converges to zero, at the optimal convergence rate of the underlying numerical method. Numerical solutions are evolved in time using the high order Taylor-type time stepping scheme of the same order of accuracy of the spatial discretization. By combining the DGSEM spatial approximation with the high order Taylor-type time stepping scheme and the accuracy of the PML we obtain an arbitrarily accurate wave propagation solver in the time domain. Numerical experiments are presented in 2D and 3D corroborating the theoretical results.solutions everywhere. The well-posedness and temporal stability analysis of the PML has been considered extensively in the literature [11,10,9,17]. For general systems, there is no guarantee that all solutions decay with time. In [9], however, the geometric stability condition was introduced to characterize the temporal stability of initial value problems for PMLs. If this condition is not satisfied, then there are modes of high spatial frequencies with temporally growing amplitudes. This result have been extended to PML initial boundary value problems (IBVPs) [1,11]. Even when the geometric stability condition is satisfied, however, numerical experiments have also shown that the PML can be unstable [32,33]. For models that satisfy the geometric stability condition, like the acoustic wave equation, recent results [1,2] have revealed the impact of numerical boundary procedures on the stability of discrete PMLs, using high order summation-by-parts (SBP) finite difference method.By the results in [11...