  A binary neutron star merger has been observed in a multi-messenger detection
of gravitational wave (GW) and electromagnetic (EM) radiation. Binary neutron
stars that merge within a Hubble time, as well as many other compact binaries,
are expected to form via common envelope evolution. Yet five decades of
research on common envelope evolution have not yet resulted in a satisfactory
understanding of the multi-spatial multi-timescale evolution for the systems
that lead to compact binaries. In this paper, we report on the first successful
simulations of common envelope ejection leading to binary neutron star
formation in 3D hydrodynamics. We simulate the dynamical inspiral phase of the
interaction between a 12$M_\odot$ red supergiant and a 1.4$M_\odot$ neutron
star for different initial separations and initial conditions. For all of our
simulations, we find complete envelope ejection and final orbital separations
of $a_{\rm f} \approx 1.3$-$5.1 R_\odot$ depending on the simulation and
criterion, leading to binary neutron stars that can merge within a Hubble time.
We find $\alpha_{\rm CE}$-equivalent efficiencies of $\approx 0.1$-$2.7$
depending on the simulation and criterion, but this may be specific for these
extended progenitors. We fully resolve the core of the star to $\lesssim 0.005
R_\odot$ and our 3D hydrodynamics simulations are informed by an adjusted 1D
analytic energy formalism and a 2D kinematics study in order to overcome the
prohibitive computational cost of simulating these systems. The framework we
develop in this paper can be used to simulate a wide variety of interactions
between stars, from stellar mergers to common envelope episodes leading to GW
sources.
