The Evolution of Massive Stellar Multiplicity in the Field I. Numerical simulations, long-term evolution and final outcomes

We investigate how the multiplicity of binary, triple and quadruple star systems changes as the systems evolve from the zero-age main-sequence to the Hubble time. We find the change in multiplicity fractions over time for each data set, identify the number of changes to the orbital configuration and the dominant underlying physical mechanism responsible for each configuration change. Finally, we identify key properties of the binaries which survive the evolution. We use the stellar evolution population synthesis code Multiple Stellar Evolution (MSE) to follow the evolution of $3 \times 10^4$ of each 1+1 binaries, 2+1 triples, 3+1 quadruples and 2+2 quadruples. The coupled stellar and orbital evolution are computed each iteration. The systems are assumed to be isolated and to have formed in situ. We generate data sets for two different black hole natal kick mean velocity distributions (sigma = 10 km/s and sigma = 50 km/s and with and without the inclusion of stellar fly-bys. Our fiducial model has a mean black hole natal kick velocity if sigma = 10 km/s and includes stellar fly-bys. Each system has at least one star with an initial mass larger than 10 solar masses. All data will be publicly available. We find that at the end of the evolution the large majority of systems are single stars in every data set (> 85%). As the number of objects in the initial system increases, so too does the final non-single system fraction. The single fractions of final systems in our fiducial model are 87.8 $\pm$ 0.2 % for the 2 + 2s, 88.8 $\pm$ 0.3 % for the 3 + 1s, 92.3 $\pm$ 0.2 % for the 2 + 1s and 98.9 $\pm$ 0.3 % for the 1 + 1s.