Hamiltonian Monte Carlo (HMC) is a popular sampling method in Bayesian inference. Recently, Heng & Jacob (2019) studied Metropolis HMC with couplings for unbiased Monte Carlo estimation, establishing a generic parallelizable sampling scheme. However, in practice a different HMC method, multinomial HMC, is considered as the go-to method, e.g. as part of the no-U-turn sampler (NUTS). In multinomial HMC, proposed states are not limited to the end-points as in Metropolis HMC; instead points along the entire trajectory can be proposed. In this paper, we establish couplings for multinomial HMC, based on optimal transport for the multinomial sampling in the transition. We prove an upper bound for the meeting time -- the time it takes for the coupled chains to meet -- using analysis based on local contractivity. We evaluate our methods using three targets: a high-dimensional Gaussian, Bayesian logistic regression and log-Gaussian Cox point processes. Compared to Heng & Jacob (2019), we find that coupled multinomial HMC is more robust to the choice of step sizes and trajectory lengths and generally attains a smaller meeting time with a better variance-bias trade-off. The robustness allows re-use of existing HMC adaptation methods and, together with other improvements, paves the way for a wider and more practical use of coupled HMC methods.