The recent detection of gravitational waves from merging neutron star events has opened a new window on the many unknown aspects of their internal dynamics. A key role in this context is played by the transition from baryon to quark matter described in the neutron star equation of state (EoS). In particular, the binary pulsar observation of heavy neutron stars requires appropriately stiff dense matter in order to counter gravitational collapse, at variance with the predictions of many phenomenological quark models. On the other side, the LIGO observations favor a softer EoS therefore providing a lower bound to the equation stiffness. We introduce a quantum chromodynamics (QCD) description of the neutron star's high baryon density regime where the pressure and energy density distributions are directly obtained from the matrix elements of the QCD energy momentum tensor. Recent ab initio calculations allow us to evaluate the energy-momentum tensor in a model independent way including both quark and gluon degrees of freedom. Our approach is a first effort to replace quark models and effective gluon interactions with a first principles, fully QCD-based description. Most importantly, the QCD energy momentum tensor matrix elements are connected to the Mellin moments of the generalized parton distributions which can be measured in deeply virtual exclusive scattering experiments. As a consequence, we establish a connection between observables from high energy experiments and from the analysis of gravitational wave events. Both can be used to mutually constrain the respective sets of data. In particular, the emerging QCD-based picture is consistent with the GW170817 neutron star merger event once we allow a first-order phase transition from a low-density nuclear matter EoS to the newly-constructed high-density quark-gluon one.