In our manuscript, we present a detailed multi-domain modeling approach for the numerical study of free-running harmonic frequency comb (HFC) formation in terahertz quantum cascade lasers (QCLs). We investigate the influence of the chosen eigenstate basis on the gain spectrum and present self-consistent simulation results of stable HFC operation in a double metal terahertz QCL. In our simulations, the studied QCL gain medium shows self-starting harmonic mode-locking for different bias and waveguide configurations, resulting in a mode spacing of up to eight times the cavity round trip frequency. Furthermore, we characterize the spectral time evolution of the coherent HFC formation process, yielding the spontaneous build-up of a dense multimode state which is gradually transferred into a broad and clear HFC state.