We report the implementation of the relativistic equation-of-motion coupled-cluster method to calculate double ionization spectra (DI-EOMCC) of the closed-shell atomic systems. This method is employed to calculate the principal valence double ionization potential values of He and alkaline earth metal (Be, Mg, Ca, Sr and Ba) atoms. Our results are compared with the results available from the national institute of standards and technology (NIST) database and other ab initio calculations. We have achieved an accuracy of ∼ 0.1%, which is an improvement over the first principles T-matrix calculations [J. Chem. Phys. 123, 144112 (2005)]. We also present results using the second-order many-body perturbation theory and the random-phase approximation in the equation-of-motion framework and these results are compared with the DI-EOMCC results.PACS numbers: 31.15.ac, 31.15.bw, 32.10.Hq Recent advances in the experimental techniques, such as x-ray free electron laser of Linac Coherent Light sources of SLAC [1,2] and attosecond pulses [3,4], have enabled studies of multi-ionization processes. The double photo-ionization of atoms in which two electrons are ejected to continuum orbitals is a three-body quantal problem. The complex interplay between the relativistic effects and the electron correlation is of central importance in the accurate description of these processes [5]. The advent of sub-femtosecond laser technology with the generation of attosecond pulses from the vacuum ultraviolet to the extreme ultraviolet wavelength region has opened up new perspectives on the observation of the correlated electron dynamics involved in the studies of the double ionization processes [6,7]. One of the outstanding theoretical problems in these studies is to explain the simultaneous double ionization mechanisms [8], which are different than the sequential ionization events. There has been experimental progress in the direction of attosecond tracing of the correlated electron-emissions in the nonsequential double ionization processes in atomic systems which requires a suitable theory that could describe the effects of the dynamical electron corrections adequately [9]. To complement the sophisticated experimental techniques, it is desirable to have accurate theoretical methods to treat the double ionization continua. Attempts have been made using the T-matrix [10,11], delta selfconsistent-field [12], and with the four-component twoparticle propagator methods [13,14]. It is well known that not only the electron correlation, but also relaxation effects plays significant role in the accurate description of atomic states. Therefore measurements and calculations based on the lower-order many-body methods do not agree with each other [15]. The equation-of- * h.pathak@ncl.res.in motion coupled-cluster (EOMCC) method [16][17][18] provides a balanced treatment of the electron correlation and relaxation effects to determine the atomic states and also calculating differences of the energies in a direct manner. It uses a large configurational spa...