Including higher-order modes in a quadrupolar eccentric numerical relativity surrogate using universal eccentric modulation functions

$\mathtt{\text{gwNRHME}}$ is a framework that converts multi-modal (i.e., containing several spherical harmonic modes) quasi-circular waveforms into their eccentric counterparts, provided the quadrupolar eccentric mode is known, by exploiting universal eccentric modulation functions. Leveraging this framework, we combine the quasi-circular NR surrogate model $\mathtt{\text{NRHybSur3dq8}}$ with the quadrupolar, non-spinning, eccentric surrogate $\mathtt{\text{NRSurE\_q4NoSpin\_22}}$ to construct a multi-modal, non-spinning, eccentric model, denoted as \model{}, which includes nine modes: $(2,\{1,2\})$, $(3,\{1,2,3\})$, $(4,\{2,3,4\})$, and $(5,5)$. When compared against 156 eccentric SXS NR waveforms, \model{} achieves median frequency-domain mismatches (computed using the Advanced LIGO design sensitivity) of $\sim 9\times {10}^{-5}$, with a standard deviation of $\sim 2\times {10}^{-4}$. To demonstrate the modularity of the framework, we further combine $\mathtt{\text{NRSurE\_q4NoSpin\_22}}$ with effective-one-body (EOB) models $\mathtt{\text{SEOBNRv5HM}}$ and $\mathtt{\text{TEOBResumS-Dali}}$ in their non-spinning limits, yielding eccentric waveforms with median mismatches of $\sim 2\times {10}^{-4}$ and $\sim {10}^{-3}$, respectively, with standard deviation of $\sim 2\times {10}^{-3}$ and $\sim 2\times {10}^{-2}$ respectively. Finally, we provide both a surrogate model, $\mathtt{\text{gwEccEvolve\_q4NoSpin\_Sur}}$, and an analytical model, $\mathtt{\text{gwEccEvNSv2}}$, for the eccentricity evolution up to $2M$ before merger, based on eccentricity definitions derived from the universal modulation functions. The $\mathtt{\text{gwNRHME}}$ framework is publicly available through the $\mathtt{\text{gwModels}}$ package, and the resulting waveform models will be released via the $\mathtt{\text{gwsurrogate}}$ package.