Abstract. We present an analysis of stationary discrete shock profiles for a discontinuous Galerkin method approximating scalar nonlinear hyperbolic conservation laws with a convex flux. Using the Godunov method for the numerical flux, we characterize the steady state solutions for arbitrary approximation orders and show that they are oscillatory only in one mesh cell and are parametrized by the shock strength and its relative position in the cell. In the particular case of the inviscid Burgers equation, we derive analytical solutions of the numerical scheme and predict their oscillations up to fourth-order of accuracy. Moreover, a linear stability analysis shows that these profiles may become unstable at points where the Godunov flux is not differentiable. Theoretical and numerical investigations show that these results can be extended to other numerical fluxes. In particular, shock profiles are found to vanish exponentially fast from the shock position for some class of monotone numerical fluxes and the oscillatory and unstable characters of their solutions present strong similarities with that of the Godunov method.Key words. discontinuous Galerkin method, discrete shock profile, scalar conservation laws, convex flux, inviscid Burgers equation, linear stability, spectral viscosity AMS subject classifications. 65N30, 65N121. Introduction. Discontinuous Galerkin (DG) methods are high-order finite element discretizations which were introduced in the early 1970s for the numerical simulation of the first-order hyperbolic neutron transport equation [20,30]. In recent years, these methods have become very popular for the solution of nonlinear convection dominated flow problems [7,8,17]. These methods allow high-order of accuracy and locality, which make them well suited to parallel computing, hp-refinement, hpmultigrid, unstructured meshes, application of boundary conditions, etc.However, the DG method suffers from spurious oscillations in the vicinity of discontinuities that develop in solutions of hyperbolic systems of conservation laws. These oscillations are due to the Gibbs phenomenon [9] and may cause the solution to become locally nonphysical leading to robustness issues for the computation. Quadrature rules are usually used to compute integrals in the discretization of the equations. Properties of the DG method therefore depend on the local structure of the numerical solution at faces and within elements of the mesh. The control of such oscillations at a reasonable cost while keeping accuracy, robustness and stability is essential for the efficiency of the DG method and remains a challenge. Strategies have been proposed such as limiters [7,8], non-oscillatory reconstructions [1,29], hp-adaptation [12], shock capturing techniques [27,10], etc. The latter methods aim at adding artificial viscosity to spread the structure of the discontinuity so that it can be resolved at the discrete level. For a DG method with polynomials of degree p > 0 and cells of size h, the resolution scale is h/p thus meaning that the m...