Supermassive black hole binary systems (SMBHs) are the most promising sources of gravitational waves for the LISA mission. However, past theoretical models and numerical simulations showed that the hardening of the binary would stall at parsec-scale separations (Final Parsec Problem, FPP) and lead to merger timescales greater than the Hubble time. This project uses large-scale, high-precision N-body simulations of a SMBH merger, with a merger-induced triaxiality of the merger remnant. The simulations show an efficient, constant hardening rate of the SMBH binary, thus successfully resolving the FPP. In this talk, I will present a novel, fully GPU-parallelized hybrid integration approach which combines direct summation with the self-consistent field (SCF) method in order to efficiently represent the galaxy merger remnant potential, while conserving computational time. Using this approach, I study the interaction and subsequent ejection of interacting stars by the SMBH binary, focusing on analyzing the high energy and orbital element changes of these stars.