Entanglement spectrum of the reduced density matrix contains information beyond the von Neumann entropy and provides unique insights into exotic orders or critical behavior of quantum systems. Here we show that strongly-disordered systems in the many-body localized phase have universal power-law entanglement spectra, arising from the presence of extensively many local integrals of motion. The power-law entanglement spectrum distinguishes many-body localized systems from ergodic systems, as well as from ground states of gapped integrable models or systems in the vicinity of scale-invariant critical points. We confirm our results using large-scale exact diagonalization. In addition, we develop a matrix-product state algorithm which allows us to access the eigenstates of large systems close to the localization transition, and discuss general implications of our results for variational studies of highly excited eigenstates in many-body localized systems. Introduction.-Recently, much progress has been made towards understanding the mechanisms of ergodicity and its breakdown in an isolated quantum many-body system. Currently, two generic classes of many-body systems are known: ergodic (thermal) systems and manybody localized (MBL) systems [1][2][3][4]. An ergodic system is one that acts as a heat bath for its subsystems, and therefore thermalizes as a result of unitary evolution [5][6][7]. By contrast, in MBL systems transport of energy is quenched by disorder, via a mechanism akin to the single-particle Anderson localization [8]. Nevertheless, MBL systems do reach stationary states [9, 10], which are highly non-thermal due to the emergence of extensively many quasi-local integrals of motion (LIOMs) [11][12][13].In addition to distinct dynamical properties, ergodic and MBL systems are sharply distinguished by the microscopic nature of their eigenstates. This difference can be probed via quantum-information measures, such as entanglement entropy (EE). Given a pure quantum state ψ of a many-body system S = L ∪ R, consisting of two subsystems L and R, the EE is defined as S = − D i λ i ln λ i , where {λ i }, i = 1, ..., D, are the eigenvalues of the reduced density matrixρ R = Tr L |ψ ψ|, and D is the dimensionality of the Hilbert space of R. The EE of highlyexcited eigenstates of thermal systems, which obey the "Eigenstate Thermalization Hypothesis" [5][6][7], is known to generically scale as the number of degrees of freedom in R ("volume law"). On the other hand, in an MBL system the EE of nearly all eigenstates obeys the "area law" [11,12,14]. This weaker scaling of EE makes MBL systems reminiscent of ground states of gapped systems [15].EE, while providing a quantitative measure of entanglement in a many-body state, contains no information about how it is created or how different degrees of freedom are entangled with each other. Therefore, to gain