\section{Introduction}
Despite very different natures, gauge theories and gravity have deep connections; one of the oldest and the most prominent example is the double copy structure \cite{Kawai:1985xq, Bern:2008qj, Bern:2010ue}. Originally it was discovered from Kawai-Lewellen-Tye (KLT) relations \cite{Kawai:1985xq} in string theory, and a modern realization of double copy has relied on the duality between color and kinematics for gauge theory amplitudes, where the Bern-Carrasco-Johansson (BCJ) kinematic numerators satisfy the same Jacobi relations as the color factors. The duality and double copy have led to tremendous progress in the study of amplitudes both in gauge theory and gravity (see~\cite{Bern:2019prr, Bern:2022wqg, Adamo:2022dcm} and references therein). More recently, the authors of~\cite{Cheung:2021zvb} have revealed a so-called covariant color-kinematics (CCK) duality for a large class of theories including Yang-Mills theory (YM) and its coupling to bi-adjoint $\phi^3$ (YMS). As a consequence, the duality implies new, closed-form expression for BCJ numerators of all tree-level amplitudes in YMS and YM theory.
Previous works on BCJ numerators and kinematic algebras include~\cite{Mafra:2011kj,Bargheer:2012gv,CHY3,He:2015wgf,Fu:2017uzt,Teng:2017tbo, Du:2017kpo, He:2017spx, Edison:2020ehu,He:2021lro, Monteiro:2011pc,Monteiro:2013rya,Cheung:2016prv,Chen:2019ywi,Edison:2022jln} and references therein.

On the other hand, recent years have seen progress on revealing new geometric/combinatorial structures underlying scattering amplitudes {\it e.g.} from the (all-loop) amplituhedron of supersymmetric Yang-Mills ~\cite{Arkani-Hamed:2013jha} to the associahedron for bi-adjoint $\phi^3$ at tree level~\cite{Arkani-Hamed:2017mur} (with extensions to string scattering~\cite{Arkani-Hamed:2019mrd}). It is natural to look for hints of such structures underlying YM and gravity amplitudes; instead of directly working with tree amplitudes, one may decompose the problem and ask a somewhat strange question as a first step: are there combinatorial structures underlying BCJ numerators?

In this note, we take BCJ numerators from CCK duality~\cite{Cheung:2021zvb} as inputs and present preliminary evidence for such structures: in addition to the more familiar quasi-shuffle Hopf algebras~\cite{Hoffman}, we find hidden combinatorial permutohedra \cite{wiki:Permutohedron} for BCJ numerators. Any BCJ numerator can be written as the sum over all boundaries of a permutohedron (or terms from a quasi-shuffle product); for a co-dimension $d$ boundary (length-$d$ term), it contains a product of $d{+}1$ factors each with a spurious pole and a gauge-invariant numerator. We will focus on the case with two scalars and $n{-}2$ gluons, which corresponds to a $(n{-}3)$-dimensional permutohedron, and it has Fubini number ${\cal F}_{n{-}2}$ boundaries with co-dimensions $d=0, 1, \ldots, n{-}3$; each boundary is labeled by $d+1$ subsets, and for each factor labeled by such a set both the numerator (which is gauge invariant in the gluons) and the (spurious pole) denominator are given by Lorentz products of momenta and polarizations, as well as the reference momenta. Apart from being the most illustrative BCJ numerators of YMS cases, we will also see that they give nice BCJ numerators in the heavy-mass effective theory (HEFT)~\cite{Georgi:1990um,Luke:1992cs,Neubert:1993mb,manohar_wise_2000,Damgaard:2019lfh,Brandhuber:2021kpo, Brandhuber:2021bsf, Brandhuber:2022enp} as well as decoupling limit into pure YM amplitudes. BCJ numerators in HEFT have attracted lots of interest recently for their roles in the computation of gravitational amplitudes for black-hole scattering and gravitational waves~\cite{Brandhuber:2021eyq} ({\it c.f.} ~\cite{Kosower:2018adc,Bern:2019nnu,Damour:2019lcq,Bern:2021dqo,DiVecchia:2021bdo,Herrmann:2021tct,Bjerrum-Bohr:2021din,Bjerrum-Bohr:2021wwt,Jakobsen:2021lvp} for some recent works). We will take the heavy-mass limit of YMs amplitude, and (as we have checked up to $n=10$) highly nontrivial cancellations lead to a nice formula for BCJ numerators in HEFT which corresponds to ${\cal P}_{n{-}3}$ (one dimensional lower)!

Furthermore, our results imply new recursion relations and surprisingly, ``factorization'' properties of BCJ numerators on facets of permutohedra; all these can be extended to BCJ numerators of general YMS amplitudes, which in turn combine into a formula for the YM case as well. For the latter, we can then turn the logic around: since the BCJ numerators are manifestly gauge invariant in $n{-}1$ gluons, by showing that all spurious poles indeed cancel in the amplitude based on such ``factorizations", it follows from the uniqueness theorem of~\cite{Arkani-Hamed:2016rak} that they must give correct YM and gravity amplitude (after double copy) even without knowing the CCK duality.

Let us consider color-ordered YMS amplitude $A(1^\phi,2,\ldots,n{-}1,n^\phi)$ with scalars $1^\phi,n^\phi$. Its expansion onto the Kleiss-Kuij(KK) basis~\cite{Kleiss:1988ne} of bi-adjoint $\phi^3$ amplitudes has  BCJ master numerators as coefficients reads
\begin{equation}\label{eq: BCJnum}
 A(1^\phi,2,\ldots,n{-}1,n^\phi){=}\sum_{\beta\in S_{n-2}} K(1,\beta,n)A^{\phi^3}(1,\beta,n),
\end{equation}
where the sum is over $(n{-}2)!$ permutations of gluons and $A^{\phi^3}(1,\beta,n) \equiv m(1,2,\ldots,n | 1,\beta,n)$ denotes bi-adjoint $\phi^3$ amplitudes with the first ordering fixed to be $(1,2,\ldots,n)$.

Remarkably, the BCJ numerators from CCK duality $K(1,\beta,n)$ respect the Bose symmetry of all the $n{-}2$ gluons ~\cite{Cheung:2021zvb}: we only need a single numerator with the ordering chosen to be $\beta=(2,\ldots,n{-}1)$, and all others can be obtained by relabelling; they are also gauge invariant for the gluons, which becomes manifest since the dependence on polarizations is through Lorentz products of linearized field strengths $F_i^{\mu\nu}\equiv p_i^\mu \varepsilon_i^\nu-p_i^\nu\varepsilon_i^\mu$
\begin{equation}\label{Fprod}
[F_{\sigma}]^{\mu \nu}= [F_{\sigma_{1}}\cdot F_{\sigma_{2}}\cdots \cdot F_{\sigma_{|\sigma|}}]^{\mu \nu}
\end{equation}
for an ordered subset $\sigma$. The price to pay for these desirable properties is the presence of $2^{n{-}2}-1$ spurious poles, one for each nonempty subset $I\subset \{2, \ldots, n{-}1\}$:
\begin{equation}\label{Ddef}
  D_{I}:= p_{I}\cdot q_{I},\quad  \text{with} \quad p_I:=\sum_{i\in I}p_i,
\end{equation}
which depends on a reference momentum  $q_{I}$. These numerators can be simplified with some choices of $q_I$, and the final amplitude is independent of them.

\section{The permutohedron and algebra underlying BCJ numerators}
In this section, we show that all the terms in a BCJ master numerator obtained from CCK duality for YMS amplitudes are in one-to-one correspondence with all boundaries of permutohedron $\mathcal{P}_{n{-}2}$, or equivalently terms from a quasi-shuffle product.

\subsection{The (combinatorial) permutohedra and quasi-shuffle products}\label{sec:def}
Following~\cite{Cheung:2021zvb}, we organize $K(1,2,\ldots, n)$ according to the spurious pole structure, which is isomorphic to the boundary structure of the permutohedron $\mathcal{P}_{n-2}$.

The permutohedron $\mathcal{P}_{n-2}$ is an $(n{-}3)$-dimensional polytope~\cite{wiki:Permutohedron}, whose co-dimension $d$ boundary $\Gamma_d$ can be labeled by $d{+}1$ consecutive subsets
\begin{equation}\label{eq: lBD}
    \Gamma_d:=\{I_{0}, I_{1},\ldots, I_{d}\},
\end{equation}
where $I_d{\neq}\emptyset$ and $I_d\subset I_{d{-}1}\subset\ldots\subset I_0=\{2,3,\ldots,n{-}1\}$; the interior of ${\cal P}_{n{-}2}$ can be viewed as its co-dimension $0$ boundary, $\Gamma_0:=I_0$. ${\cal P}_{n{-}2}$ and its boundaries have appeared in the context of cubic tree graphs from the worldsheet~\cite{Gao:2017dek, He:2018pue}. Here each term of the BCJ numerator $K(1, 2, \ldots, n)$ with $d{+}1$ spurious poles corresponds to such a co-dimension $d$ boundary, thus the numerator can be expanded in terms of boundaries of $\mathcal{P}_{n-2}$
\begin{equation}\label{eq: KBDexp}
K(1,2,\ldots,n)=\sum_{d=0}^{n{-}3}\sum_{\Gamma_d\in\partial^d\mathcal{P}_{n{-}2}}K_{\Gamma_d}(1,2,\ldots,n),
\end{equation}
where we sum over all boundaries $\Gamma_d\in\partial^d\mathcal{P}_{n{-}2}$ with co-dimension $d=0, \ldots, n{-}3$, and the contribution from $\Gamma_d$, $K_{\Gamma_d}(1,2,\ldots,n)\equiv K_{\Gamma_d}$ reads
\begin{equation}\label{eq: KBD}
K_{\Gamma_d}=\prod_{k=0}^{d}\frac{p_{1
\Delta(I_{k},I_{k+1})
}\cdot F_{\tau_{k}}\cdot q_{I_{k}}}{D_{I_{k}}}.
\end{equation}
It has $d{+}1$ factors each with a denominator $D_{I_k}$ of \eqref{Ddef} and a numerator of the form $p_{1 \Delta(I_{k},I_{k+1})}\cdot F_{\tau_k} \cdot q_{I_k}$ for $k=0, \dots, d$. To specify the ordered subset $\tau_k$ of the Lorentz product as in \eqref{Fprod}, we introduce an alternative form of \eqref{eq: lBD} using ordered sets:
\begin{align}\label{eq: lBDt}
   \Gamma_d
   =&\{\tau_0\cup\tau_1\cup\ldots\cup\tau_d,\tau_1\cup\ldots\cup\tau_d,\ldots,\tau_{d{-}1}\cup\tau_d,\tau_d\}\\\nonumber
   \sim&\{\tau_0,\tau_1,\ldots,\tau_d\},
\end{align
where the first line is equivalent to \eqref{eq: lBD} but we use ordered sets $\tau_{k}=\mathrm{Id}(I_k/I_{k{+}1})$ with $I_{d{+}1}\equiv \emptyset$. $\mathrm{Id}(I)$ means sorting the subset $I$ in numerical ordering~\footnote{For $K(1,\beta,n)$, everything stays the same except that in \eqref{eq: lBDt} the definition of $\tau_i$'s become $\tau_{k}=\beta(I_k/I_{k{+}1})$, {\it i.e.} any subset $I$ is sorted according to $\beta$ ordering.}. Perhaps the most subtle point is that we also define
$\Delta(I_{k},I_{k+1}):=\left(\bar{I}_{k}\right)_{<\tau_{{k},1}}$, which refers to the elements in the set $\bar{I}_{k}=\{2,3,\ldots,n-1\}/I_k$ that are numerically smaller than the first element $\tau_{k,1}$ of the ordered set $\tau_{k}$.

For example, at $n=5$, we have a boundary $\Gamma_2=\{I_0=\{2,3,4\}, I_1=\{2,4\}, I_2=\{4\}\}$; equivalently, we have $\tau_0=\{3\}, \tau_1=\{2\}, \tau_2=\{4\}$, thus we have a term
\begin{equation}
K_{234,24,4}= \frac{p_1\cdot F_3\cdot q_{234}  p_1\cdot F_2\cdot q_{24} p_{123}\cdot F_4\cdot q_4 }{D_{234}D_{24} D_4  }
\end{equation}
where we have used $\Delta(I_2,\emptyset)=\left.\left(\bar{I_2}\right)\right|_{<4}=\{2,3\}$. A more nontrivial example for the latter is for $n{=}9$, $\Delta(\{4,5,7,8\},\{4,8\})=\left.\{2,3,6\}\right|_{<5}=\{2,3\}$.

On the other hand, the boundaries of $\mathcal{P}_ {n-2}$ have a nice quasi-shuffle product interpretation. The quasi-shuffle product $\star$ can be defined between two arbitrary generators $\left(\sigma_0,\sigma_1,\ldots,\sigma_r\right)$ and $\left(\rho_0,\rho_1,\ldots,\rho_s\right)$, where $\sigma_i$ and $\rho_j$ are sets with arbitrary lengths and can also be prompted to the quasi-shuffle Hopf algebra~\cite{Hoffman, Brandhuber:2021bsf, Brandhuber:2022enp, Chen:2022nei}. We summarize the definitions in appendix \ref{sec:appAG} and here we just use the following result of the quasi-shuffle product $\hat{K}(2,\ldots, n{-}1)\equiv(2)\star(3)\star\ldots\star(n-1)$
\begin{equation}\label{eq: part}
    \hat{K}(2,\ldots,n{-}1)
    {=}\sum_{d{=}0}^{n{-}3}\sum_{\tau\in \text{part}^{(d{+}1)}(2,\ldots,n{-}1)}\!\!\!\!\!\!\!\!\!\!\!\!\!\!(-1)^{n+d-1}(\tau_0,\tau_1,\ldots,\tau_{d}),
\end{equation}
where $\text{part}^{(d+1)}(2,\ldots,n{-}1)$ denotes all the ordered partitions of $\{2,3,\ldots,n{-}1\}$ into $d+1$ nonempty subsets $(\tau_0,\tau_1,\ldots,\tau_d)$ (each $\tau_i$ is sorted according to $\beta$).  Terms on the RHS of \eqref{eq: part} are in one-to-one correspondence with boundaries of $\mathcal{P}_{n-2}$ as in~\eqref{eq: lBDt}, thus we can rewrite \eqref{eq: KBDexp} in terms of the quasi-shuffle product
\begin{equation}
\label{eq:qsp}
K(1,2,\ldots,n)=\langle\hat{K}(2,\ldots n{-}1)\rangle,
\end{equation}
where we have defined a {\it linear map} $\langle\cdot\rangle$ from \eqref{eq: KBD} for any partition $(\tau_0,\tau_1,\ldots,\tau_{d})$
\begin{equation}\label{eq:linmap}
    \langle(\tau_0,\tau_1,\ldots,\tau_{d})\rangle{=}(-1)^{n+d-1}\prod_{k=0}^{d}\frac{p_{1\Delta(I_{k},I_{k+1})
}\cdot F_{\tau_{k}}\cdot q_{I_{k}}}{D_{I_{k}}}.
\end{equation}

\subsection{The counting and some examples}
By definition, the permutohedron $\mathcal{P}_m$ contains $m!$ vertices and $2^{m}-2$ co-dimension one facets. More generally, the number of co-dimension $d$ boundaries of this polytope is $(d+1)!S(m,d+1)$, where $S(m,d)$ is the second kind of Stirling number \cite{wiki:Stirling_numbers_of_the_second_kind}.  Algebraically, $S(m,d+1)$ also counts the number of ways to partition a set of $m$ labeled objects into $d+1$ nonempty unlabeled subsets $\{\tau_0,\tau_1,\ldots,\tau_d\}$ \cite{wiki:Stirling_numbers_of_the_second_kind}, so after considering the ordering between these sets, there are $(d+1)!S(m,d+1)$ terms in the summation for any $d$.
The total number is the Fubini number $\mathcal{F}_{m}$, where $\mathcal{F}_{m}=\sum_{d=1}^{m}d!S(m,d)$ \cite{Mezo}, thus the $n$-point BCJ numerator has $\mathcal{F}_{n-2}$ terms.
\begin{table}[H]
  \centering
    \begin{tabular}{|c|c|c|c|c|c|}
    \hline
    \diagbox[innerwidth=0.5cm]{n}{d}      & 0     & 1     & 2     & 3 &Total \\
    \hline
    3     & 1 &       &       &  &1\\
    \hline
    4     & 1 & 2 &       &  &3\\
    \hline
    5     & 1 & 6 & 6 & &13 \\
    \hline
    6     & 1 & 14 & 36 & 24&75 \\
    \hline
    \end{tabular}%
  \label{tab:BDc}%
  \caption{Counting co-dimension-$d$ boundaries of $\mathcal{P}_{n-2}$}
\end{table}%

Let us illustrate  \eqref{eq: KBDexp} and \eqref{eq: KBD} with some examples. The most trivial case is $n=3$, where the BCJ numerator corresponds to the zero-dimensional permutohedron $\mathcal{P}_{1}$ which is just a point. It contains one term with $\Gamma_0=\{I_0\}=\{2\}$, thus $K(1,2,3)=K_{2}(1,2,3)=\frac{p_1\cdot F_{2}\cdot q_2}{D_2}$, where we have used $\Delta(I_{0},I_{1})=\emptyset$, this is generally true since $\bar{I_0}=\emptyset$.

For $n=4$, the permutohedron $\mathcal{P}_{2}$ is a line segment, where the interior ($d=0$) is labelled by $I_0=\{23\}$, and the two vertices ($d=1$) are labeled by $\{23,2\}$ and $\{23,3\}$; we show these three terms in figure \ref{fig:P2}.

\begin{figure}[H]
    \centering
    \subfloat[]{ \label{fig:P2}
    \begin{minipage}[c]{1\linewidth}
    \centering
    \includegraphics[scale=1.5]{figs/P2K.pdf}
    \end{minipage}
    }

    \subfloat[]{ \label{fig:P3}
    \begin{minipage}[c]{1\linewidth}
    \centering
    \includegraphics[scale=1.3]{figs/P3D.pdf}
    \end{minipage}
    }
    \caption{Permutohedra ${\cal P}_2$ for $K(1,2,3,4)$(top) and ${\cal P}_3$ for  $K(1,2,3,4,5)$(bottom)}
\end{figure}

Equivalently, in \eqref{eq: part} the partition $\text{part}^{(1)}$ of $\{2,3\}$ has $(\{\tau_0=\{2, 3\})$ and $\text{part}^{(2)}$ has  $(\tau_0=\{2\},\tau_1=\{3\})$ and $(\tau_0=\{3\},\tau_1=\{2\})$: they are nothing but the interior and the two vertices, according to \eqref{eq: lBDt}.

Thus the BCJ numerator $K(1,2,3,4)$ has three terms, $K_{23}$, $K_{23,2}$ and $K_{23,3}$, which read
\begin{align} \label{eq: BCJnum4}
\frac{p_{1}{\cdot}F_{23}{\cdot} q_{23}}{D_{23}}{+}\frac{p_{1}{\cdot} F_{3}{\cdot} q_{23} p_{1}{\cdot} F_{2}{\cdot} q_{2}}{D_{23}D_{2}}{+}\frac{p_{1}{\cdot} F_{2}{\cdot} q_{23} p_{12}{\cdot} F_{3}{\cdot} q_{3} }{D_{23}D_{3}
\end{align}
Notice that the last term in the above equation is from the boundary $\{23,3\}$, so the second factor in the numerator is $p_{1\left(\overline{\{3\}}\right)_{<3}}\cdot F_{3}\cdot q_{3}=p_{12}\cdot F_{3}\cdot q_{3}$. Meanwhile \eqref{eq: BCJnum4} shows that the four-point numerator $K(1,2,3,4)$ has an overall pole $D_{I_0}=D_{23}$. This can be easily seen from \eqref{eq: lBD} since the first set of any co-dimension $d$ boundary $\Gamma_d$ is always labeled by $I_0=\{2,3,\ldots,n{-}1\}$.

Notice that \eqref{eq: KBD} means each term in the BCJ numerator contains spurious poles, and the co-dimension $d$ contribution will have $d{+}1$ spurious poles where $D_{I_0}$ is an overall pole for every boundary.
Except for the overall one, the simple poles can be written as $D_I$ where $I$ is a nonempty proper subset of $I_0=\{2,3,\ldots,n{-}1\}$ and two simple poles $D_I$ and $D_J$ are compatible if and only if $I\subset J$ or $J\subset I$. For example, at five-point, except for the overall $D_{234}$, the simple poles are $D_{2},D_{3},D_{4},D_{23},D_{24},D_{34}$ which correspond to the six co-dimension one boundaries of $\mathcal{P}_3$, and the compatible double poles are
\begin{equation}
    \{D_{23}D_{2},D_{24}D_{2},D_{23}D_{3},D_{34}D_{3},D_{24}D_{4},D_{34}D_{4}\},
\end{equation}
which correspond to six vertices of the hexagon $\mathcal{P}_3$. We show the boundary contribution formally in figure \ref{fig:P3}.
These $13$ terms form a two-dimensional polytope $\mathcal{P}_{3}$. These boundaries can also be realized in quasi-shuffle product $\hat{K}(2,3,4)$, which will be discussed in appendix \ref{sec: N5}. To be precise, we give some explicit examples of different co-dimension here
\begin{equation}
\begin{aligned}
    &K_{234}{=}\frac{p_1{\cdot}F_{234}\cdot q_{234}}{D_{234}},\\
    &K_{234,23}{=}\frac{p_1{\cdot}F_4{\cdot}q_{234} p_1{\cdot} F_{23}{\cdot}q_{23}}{D_{234}D_{23}},\;\,K_{234,2}{=}\frac{p_1{\cdot} F_{34}{\cdot} q_{234} p_1{\cdot}F_2{\cdot}q_2}{D_{234}D_2}\\
    &K_{234,23,2}{=}\frac{ p_1{\cdot}F_4{\cdot}q_{234} p_1{\cdot} F_3{\cdot} q_{23} p_1{\cdot} F_2{\cdot} q_2}{D_{234}D_{23}D_2 }.
\end{aligned}
\end{equation}
The complete result for the BCJ numerator $K(1,2,3,4,5)$ is shown in the appendix \ref{sec: N5}.

Moreover, we emphasize that all the spurious poles are canceled in the final amplitude, and the amplitude does not depend on the reference momenta. The proof will be put into the following paper \cite{toapp}.

For $n=6$, $\mathcal{P}_4$ is a three-dimensional truncated octahedron shown in figure \ref{fig: P4}. As we have counted, it contains $14$ co-dimension one boundaries (six squares and eight hexagons), $36$ edges and $4!=24$ vertices, thus $75$ boundaries in total. Some terms with co-dimension $d=0,1,2,3$ are
\begin{align}
    &K_{2345}=\frac{p_1{\cdot} F_{2345}{\cdot} q_{2345}}{D_{2345}}\\\nonumber
    &K_{2345,234}=\frac{p_1{\cdot} F_{5}{\cdot} q_{2345}p_1{\cdot} F_{234}{\cdot} q_{234}}{D_{2345}D_{234}}\\\nonumber
    &K_{2345,234,23}=\frac{p_1{\cdot} F_{5}{\cdot} q_{2345}p_1{\cdot} F_{4}{\cdot} q_{234}p_1{\cdot} F_{23}{\cdot} q_{23}}{D_{2345}D_{234}D_{23}}\\\nonumber
    &K_{2345,234,23,2}=\frac{p_1{\cdot} F_{5}{\cdot} q_{2345}p_1{\cdot} F_{4}{\cdot} q_{234}p_1{\cdot} F_{3}{\cdot} q_{23}p_1{\cdot} F_2{\cdot} q_2}{D_{2345}D_{234}D_{23}D_2}
\end{align}
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.4]{figs/P4label.pdf}
    \caption{The permutohedron $\mathcal{P}_4$ for $K(1,\ldots, 6)$}
    \label{fig: P4}
\end{figure}
\section{Numerators for general YMS and pure YM amplitudes}
\label{sec:YM}
More generally, the CCK duality has provided closed-formulae for BCJ numerators of $n$-point YMS amplitude with $r\geq 2$ scalars~\cite{Cheung:2021zvb}. It turns out that any such numerator corresponds to a permutohedron ${\cal P}_{n{-}r}$ (with dimension $n{-}r{-}1$): everything we have discussed above for $r=2$ case still applies if we replace $\{2, 3,\ldots, n{-}1\}$ by the set of $n{-}r$ gluons. For example, in the other extreme with $r=n$ we can formally define ${\cal P}_0:={\cal P}_{\emptyset}$ and the $n$-scalar numerator (of bi-adjoint $\phi^3$ amplitude) is $1$ or $0$. For $r=n{-}1$, we have the zero-dimensional ${\cal P}_1$ and the numerator with a single gluon $i$ reads
\begin{equation}
K^{1-{\rm gluon}}(1, \beta, n)=\frac{p_{1{\rightarrow }i} \cdot F_i \cdot q_i}{D_i}
\end{equation} where $1 {\rightarrow} i$ denotes all the scalars preceding $i$ in $(1\beta n)$. We shall not repeat this for general cases but leave detailed discussions to a separate paper~\cite{toapp}.

Moreover, since the pure YM amplitude can be expanded as a linear combination of these YMS ones~\cite{Lam:2016tlk, Fu:2017uzt, Du:2017kpo, Cheung:2017ems, Dong:2021qai}, we obtain its BCJ numerators for free; the resulting numerator naively contains $2\mathcal{F}_{n-2}$ terms as derived in \cite{Cheung:2021zvb}. However, we can still organize the terms according to pole structures and immediately combine them in pairs as $\mathcal{F}_{n-2}$ terms: the resulting numerator has the same form as the $2{-}$scalar case and corresponds to boundaries of permutohedron $\mathcal{P}_{n-2}$. By expanding $A^{\text{YM}}(1,2,\ldots,n)$ in exactly the same way as \eqref{eq: BCJnum}
each master BCJ numerator, {\it e.g.}  $K^{\text{YM}}(1,2,\ldots,n)$ is given by a sum over boundaries of $\mathcal{P}_{n-2}$ as in \eqref{eq: KBDexp}:
\begin{align}
    K^{\text{YM}}(1,2,\ldots,n)&=\sum_{d=0}^{n{-}3}\sum_{\Gamma_d\in\partial^d\mathcal{P}_{n{-}2}}K^{\text{YM}}_{\Gamma_d}(1,2,\ldots,n). \nonumber
\end{align}
where the contribution from each boundary
is identical to \eqref{eq: KBD} except for the $k=0$ factor which becomes
\begin{equation}\label{eq: N0YM}
    \frac{\varepsilon_n\cdot F_{1\tau_0}\cdot q_{I_{0}}+\varepsilon_1\cdot F_{\tau_0}\cdot\left(\varepsilon_n p_{1n}\cdot q_{I_{0}}-q_{I_{0}} p_{1}\cdot\varepsilon_n\right)}{D_{I_0}}
\end{equation}

Of course, similar to the YMS case, all spurious poles cancel in the final amplitude, which does not depend on $q_I$. Therefore we are free to choose them to simplify the expression \eqref{eq: N0YM}. One such choice is $q_{I_{0}}=\varepsilon_n$, and the $k=0$ factor \eqref{eq: N0YM} takes a simpler form
\begin{equation}\label{eq: KYMSim}
    \frac{\varepsilon_n{\cdot} F_{1\tau_0}{\cdot} \varepsilon_n}{\varepsilon_n{\cdot} p_{23\ldots n{-}1}}=-\frac{\varepsilon_n{\cdot} F_{1\tau_0}{\cdot} \varepsilon_n}{\varepsilon_n{\cdot} p_{1}}.
\end{equation}
It is easy to see that the BCJ numerators become manifestly gauge invariant in particles $1,2,\ldots,n{-}1$.

For example, the BCJ numerator $K^{\text{YM}}(1,2,3,4)$ reads
\begin{equation*}
    {-}\frac{\varepsilon_{4}{\cdot}F_{123}{\cdot} \varepsilon_4}{p_1{\cdot}\varepsilon_4}{-}\frac{\varepsilon_4{\cdot} F_{13}{\cdot} \varepsilon_4 p_{1}{\cdot} F_{2}{\cdot} q_{2}}{p_1{\cdot}\varepsilon_4 D_{2}}{-}\frac{\varepsilon_4{\cdot} F_{12}{\cdot} \varepsilon_4 p_{12}{\cdot} F_{3}{\cdot} q_{3} }{p_1{\cdot}\varepsilon_4 D_{3}}.
\end{equation*}

Furthermore, similar to the discussion in section \ref{sec:def}, BCJ numerators of YM amplitudes can also be interpreted in terms of quasi-shuffle products, and the only change is that in the linear map \eqref{eq:linmap} the $k=0$ factor is modified to  \eqref{eq: N0YM}.

Before ending the section, we mention the obvious double-copy from YM to GR
\begin{equation}
M^{\rm GR}_n=\sum_{\alpha, \beta} K^{\rm YM}(1, \alpha, n) m(1,\alpha ,n|1, \beta ,n) K^{\rm YM}(1, \beta, n)
\end{equation}
where we sum over a pair of permutations $\alpha,\beta$ of $\{2,3,\ldots, n{-}1\}$, with $m$ denoting bi-adjoint $\phi^3$ amplitudes; if we replace YM by YMS with $1, n$ being scalars, it gives the amplitude with $n{-}2$ gravitons and two scalars.

\section{Recursions and factorizations}
In this section, we propose recursion relations and factorization properties (on spurious poles $D_I$) for the BCJ numerators, which are implied by the combinatorial and algebraic structure. The argument can be equally applied to both two-scalar YMS and pure YM numerators.

\subsection{Recursion relations}
First, in quasi-shuffle product~\eqref{eq: part}, one can collect the terms with the same $\tau_d$ and then apply the linear map \eqref{eq:linmap} to obtain the following recursion relation,
\begin{equation}\label{eq: Krec}
    K(1,2,\ldots,n){=
    \sum_{I\subset\{2,\ldots,n{-}1\}}\frac{p_{1\Delta(I,\emptyset
    }{\cdot} F_I {\cdot} q_I}{D_I}\tilde{K}^I(1,\bar{I},n),
\end{equation}
where the summation is over all the nonempty subsets of $\{2,3,\ldots,n{-}1\}$. The definition of $\tilde{K}^I(1,\bar{I},n)$ is slightly different from \eqref{eq: KBD} in the denominator: it is given by the sum over boundaries of the permutohedron $\mathcal{P}_{\bar{I}}$ with vertices labeled by all permutations of set $\bar{I}$, and for each boundary $\Gamma_d=\{J_0=\bar{I},J_1,\ldots,J_{d}\}$ where $J_d\subset J_{d{-}1}\ldots\subset J_0$, we have a  contribution
\begin{equation}\label{eq: tKBD}
\tilde{K}^I_{\Gamma_d}=\prod_{k=0}^{d}\frac{p_{1\Delta(J_{k},J_{k+1})
}\cdot F_{\tau_k}\cdot q_{IJ_{k}}}{D_{IJ_{k}}},
\end{equation}
where $D_{IJ_k}=p_{IJ_k}\cdot q_{IJ_k}$, $\tau_k= \mathrm{Id}(J_k/J_{k{+}1})$ with $J_{d{+}1}\equiv \emptyset$ and the complement of the set $J_k$ appears in $\Delta(J_{k},J_{k+1})$ is defined as $\bar{I}/J_k$. For $\abs{I}{=}n{-}2$ ($\bar{I}=\emptyset$), we define $\tilde{K}^I(1,n)=1$. Formally, this numerator corresponds to the permutohedron $\mathcal{P}_0$.

For example, the recursion relation of $K(2,3,4)\equiv K(1,2,3,4,5)$ \footnote{Here we have omitted the labels of the scalar particles $1$ and $n$ in the numerators $K$ and $\tilde{K}$.} reads
\begin{small}
\begin{align}\label{eq: rec5}
    K&(2,3,4)=\frac{p_1{{\cdot}} F_{234}{\cdot} q_{234}}{D_{234}}\\\nonumber
    {+}&\frac{p_1{\cdot} F_{23}{\cdot} q_{23}}{D_{23}}\tilde{K}^{23}(4){+}\frac{p_{1}{\cdot} F_{24}{\cdot} q_{24}}{D_{24}}\tilde{K}^{24}(3){+}\frac{p_{12}{\cdot} F_{34}{\cdot} q_{34}}{D_{34}}\tilde{K}^{34}(2)\\\nonumber
    {+}&\frac{p_1{\cdot} F_2{\cdot} q_2}{D_2}\tilde{K}^2(3,4)
    {+}\frac{p_{12}{\cdot} F_3{\cdot} q_3}{D_3}\tilde{K}^3(2,4){+}\frac{p_{123}{\cdot} F_4{\cdot} q_4}{D_4}\tilde{K}^4(2,3).
\end{align}
\end{small}
Geometrically, the recursion relation~\eqref{eq: Krec} tells us how the co-dimension one boundaries of permutohedron are glued together. In the above five-point example, the term with $\abs{I}=3$ in the first line has only one pole and corresponds to the interior (co-dimension $0$ boundary) of $\mathcal{P}_{3}$, depicted in figure~\ref{fig:rec5}. For the three terms with $\abs{I}=2$ in the second line, each factor $\tilde{K}^I(1,\bar{I},n)$ corresponds to a zero-dimensional permutohedron; on the other hand, each term is mapped to a co-dimension one boundary of $\mathcal{P}_3$ without vertices. For the remaining three terms with $\abs{I}=1$ in the last line, each $\tilde{K}^I(1,\bar{I},n)$ corresponds to a one-dimensional permutohedron and it is mapped to a co-dimension one boundary with two vertices.
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.8]{figs/rec5.pdf}
    \caption{Recursion relation at $n=5$}
    \label{fig:rec5}
\end{figure}
\subsection{Factorization properties on spurious poles}
Next, we move to certain intriguing factorization properties of the BCJ numerator on spurious poles. Combinatorially, any co-dimension one boundary of the permutohedron $\mathcal{P}_{n{-}2}$ is the product of two lower-dimensional permutohedra $\mathcal{P}_{I}\times\mathcal{P}_{\bar{I}}$. Remarkably, we find that on any pole $D_I=0$, the residue of the BCJ numerator factorizes into the product of a $(\abs{I}{+}2)$-point numerator and a $(n{-}\abs{I})$-point numerator! Unlike the usual factorization on the physical poles of the amplitude, these factorizations on the spurious poles stem from the combinatorial picture without any known physical origin. Explicitly
\begin{equation}\label{eq: Kfac}
\begin{aligned}
\left.\mathrm{Res}\right|_{D_I{=}0}K(1,2,\ldots,n){=}&D_I K( 1,\mathrm{Id}(I),P) \\
    &\times\tilde{K}^I(1,\mathrm{Id}(\bar{I}),n),
\end{aligned}
\end{equation}
where $P\equiv \bar{I}n$ denotes an effective scalar. For the definition of $\Delta(I_{k}, I_{k{+}1})$ in $K( 1,\mathrm{Id}(I),P)$, the complement of the set $I_k$ is still defined as $\bar{I_k}=\{2,3,\ldots,n{-}1\}/I_k$ while for $\Delta(J_{k}, J_{k{+}1})$ in $\tilde{K}^I(1,\mathrm{Id}(\bar{I}),n)$ the complement of $J_k$ is defined as $\bar{I}/J_k$.
The factor $D_I K(1,\mathrm{Id}(I),P)$ in \eqref{eq: Kfac} means that the overall pole $D_I$ of $K(1,\mathrm{Id}(I),P)$ is excluded. The factorization properties~\eqref{eq: Kfac} can be proved directly by plugging in the definitions on both sides.

For instance, at six points as shown in figure \ref{fig: P4}, there are $14$ co-dimension one boundaries $D_I=0$ including eight poles with $\abs{I}=1$ or $3$ corresponding to hexagons and six poles with $\abs{I}=2$ corresponding to squares. On any of the hexagon boundary, {\it i.e.} when $D_I=0$ with $\abs{I}=1$ or $3$, the residue factorizes into $\mathcal{F}_3=13$ terms (times $\mathcal{F}_1=1$ term). Similarly when $D_I=0$ with $\abs{I}=2$, the residue factorizes differently, {\it e.g.} as $D_{23}K(1,2,3,456)\times\tilde{K}^{23}(1,4,5,6)$ when $D_{23}=0$ (the square is the product of two line segments $\mathcal{P}_{\{23\}}\times\mathcal{P}_{\{45\}}$).

Algebraically, the quasi-shuffle algebra can be prompted to a bialgebra by introducing the coproduct map \cite{Hoffman}, and one can show the factorization properties from the coproduct. Similar to \cite{Brandhuber:2021bsf, Brandhuber:2022enp}, we can also define the antipode map to make the bialgebra a quasi-shuffle Hopf algebra. Acting on the BCJ numerators, the antipode map does nothing but changes its overall sign. The detail is given in the appendix \ref{sec:appAG}.

We expect the factorization properties of BCJ numerators to be the key for showing the cancellation of spurious poles in the amplitude. Such properties also suggest certain positive geometries (rather than just combinatorics) underlying these BCJ numerators, and we leave further investigations to future works.

\section{Heavy-mass effective field theory} \label{sec:heavy1}
In this section, we study YMS amplitudes and their BCJ numerators in the heavy-mass effective field theory (HEFT), which are obtained by taking the heavy-mass limit for a pair of massive scalars with momenta \cite{Brandhuber:2021kpo, Brandhuber:2021bsf,Brandhuber:2022enp}
\begin{equation}
    p_1^\mu=mv^\mu,\qquad p_n^\mu=-m v^\mu - k^\mu,
\end{equation}
where $v^2=1$ and we are interested in the limit $m \to \infty$; in other words, we will study the the expansion in $1/m$ of the BCJ numerators which we denote as $K_\mathrm{H}(1,2,\ldots,n)$, as well as that of $\phi^3$ amplitudes, which combine to give the resulting HEFT amplitude $A^\mathrm{H}(1,2,\ldots,n)$ at the leading order in $1/m$. Here $k^\mu$ is at the same order as gluon momenta, which stay finite at ${\cal O}(m^0)$ as $m \to \infty$.

\subsection{Heavy limit of YMS amplitudes}
We will make a  particular choice of the reference momenta: $q_I=v$ for all $I$, which dramatically simplifies formulae for BCJ numerators and give rise to poles similar to HEFT numerators  in~\cite{Brandhuber:2021bsf}. In fact, for $n=4$ such a choice reduces the BCJ numerator to one term, since $v\cdot F_a \cdot v$ vanishes for a single particle $a$
\begin{equation}\label{eq: hamp4}
    K_\mathrm{H}(1,2,3,4)= \frac{p_{1}\cdot F_{23}\cdot v} {p_{23} \cdot v} = -\frac{2m^2}{k^2}\ v \cdot F_{23}\cdot v,
\end{equation}
where in the second equality we have used $v\cdot k=-k^2/(2m)$ implied by the on-shell condition $p_n^2=m^2$.
Notice that $K_{\mathrm{H}}(1,3,2,4)=K_{\mathrm{H}}(1,2,3,4)$, thus the amplitude  $A^\mathrm{H}(1,2,3,4)$ becomes
\begin{align} \label{eq:Hamp4}
&(\frac{1}{s_{12}}+ \frac{1}{s_{23}}) K_\mathrm{H}(1,2,3,4)  -\frac{1}{s_{23}}K_\mathrm{H}(1,3,2,4)  \\\nonumber
=&-\frac{m}{k^2} \frac{v \cdot F_{23}\cdot v}{ v\cdot p_2}.
\end{align}
Physically, the final HEFT amplitude has the leading order $\mathcal{O}(m)$ \cite{Brandhuber:2021kpo}. In the above example, we can see that the numerators are at $\mathcal{O}(m^2)$, and the sum of the leading contribution of $\phi^3$ amplitudes at $\mathcal{O}(m^0)$, say $1/s_{23}$ times the corresponding numerators vanishes. Therefore, the sum of the contribution of $\phi^3$ amplitudes at the next order, {\it i.e.} $1/s_{12}$ from $A^{\phi^3}(1,2,3,4)$ times the numerator produces the HEFT amplitude as the first non-vanishing order. This is also the case for $n=5$.
However, for higher $n$, the numerator contains some additional terms with higher power of $m$.
To obtain the leading order contribution of HEFT final amplitude, we expand the numerators and the $\phi^3$ amplitudes in $m^{-1}$. Note the overall pole $D_{23\ldots n-1}=v\cdot k$ for BCJ numerators is proportional to $m^{-1}$, we first collect the numerator according to its {\it superficial order} of $m^{-1}$, {\it i.e.} terms with $(i-1)$ $p_1$'s in the numerator,
\begin{equation} \label{eq: Kexpand}
    K_{\mathrm{H}}(1,2,\ldots,n)= \sum_{i=2}^{\lfloor n/2 \rfloor } K_\mathrm{H}^{(i)}(1,2,\ldots,n),
\end{equation}
where the upper bound of the summation is $\lfloor n/2 \rfloor$ since $p_1 \cdot F_{a}\cdot v=0$
implies that the numerator should contain as many $p_1 \cdot F_{ab} \cdot v$ as possible to have the highest power of $p_1$. In the above expansion, $K_\mathrm{H}^{(i)}\equiv K_\mathrm{H}^{(i)}(1,2,\ldots,n)$ refers to terms with the superficial order $\mathcal{O}(m^i)$. For example, at six points we have the following terms for $K_\mathrm{H}^{(2)}$ and $K_\mathrm{H}^{(3)}$ respectively:
\begin{equation*}
       \frac{p_1\cdot F_{23} \cdot v\ p_{23}\cdot F_5 \cdot v\ p_{23}\cdot F_4 \cdot v } {v\cdot k\ v \cdot p_{45}\ v\cdot p_4 }, \frac{p_1 \cdot F_{25}\cdot v\ p_1\cdot F_{34}\cdot v} {v\cdot k\ v\cdot p_{34} }.
\end{equation*}
In fact, as explained in appendix \ref{sec:heavy}, the actual order of $K_\mathrm{H}^{(i)}$ is $\mathcal{O}(m^{2})$ for $i=2$ and $\mathcal{O}(m^{i-1})$ otherwise.

In the HEFT amplitude, we sum over all cubic graphs relevant at leading order, and for each graph with its propagator structure, its numerator is given by the corresponding commutator of $K_\mathrm{H}^{(i)}$~\cite{Bern:2011ia}. Nicely we observe that certain commutators of $K_\mathrm{H}^{(i)}$ actually vanish, and the end result is that only $K_\mathrm{H}^{(2)}$ contributes to the amplitude at the leading order! We have checked such vanishing results up to $n=10$, but we do not have an all-$n$ proof at the moment.

In fact, such vanishing results are better than what we need here, {\it i.e.} only $K^{(2)}_{\rm H}$ contributes to gauge-theory amplitudes at leading order. We have checked up to $n=10$ that the stronger vanishing results actually ensure that only $K^{(2)}_{\rm H}$ contributes to gravity amplitude, which is at order $\mathcal{O}(m^2)$, as obtained by double copy in HEFT. We leave more details in the appendix \ref{sec:heavy} with a proof of the simplest case. As a result of this conjecture, the amplitude is given by
\begin{equation} \label{eq: Hamp}
\begin{aligned}
    &A^\mathrm{H}(1,2 \ldots,n)=\sum_{\Theta^1 } \frac{K_\mathrm{H}^{(2)}(1,\Theta^1,n)}{d_{\Theta^1}}, \\
\end{aligned}
\end{equation}
where the summation is over nested commutators of depth $n{-}4$ (``co-depth" $1$) of the ordered set $(2,3,\ldots,n-1)$. For instance, at 5-point we sum over $\Theta^1=([2,3],4),(2,[3,4])$; $d_{\Theta^1}$ denotes the propagator denominator corresponding to the cubic tree associated with $\Theta^1$ (two sub-trees on the scalar line $(1 n)$):
\begin{equation*}
\begin{aligned}\includegraphics[width=0.23\linewidth]{figs/general1.pdf}\end{aligned} \leftrightarrow d_{\Theta^{1}}, {\it e.g.} \begin{aligned}\includegraphics[width=0.23\linewidth]{figs/theta511.pdf}\end{aligned} \leftrightarrow  d_{[2,3],4}=s_{123}s_{23}.
\end{equation*}
Moreover, it is easy to show (see appendix \ref{sec:heavy} for details) that the effective BCJ numerator $K_\mathrm{H}^{(2)}(1,2,\ldots,n)$ contains $\mathcal{F}_{n{-}3}$ terms, and its pole structure corresponds to the permutohedron $\mathcal{P}_{\{34\ldots n{-}1\}}$, which means that
\begin{equation}\label{eq: BDexph}
    K_\mathrm{H}^{(2)}(1,2,\ldots,n)=\sum_{d=0}^{n-2}\sum_{\Gamma_d\in\partial^d\mathcal{P}_{n{-}3}}K^{ (2)}_{H,\Gamma_d}(1,2,\ldots,n).
\end{equation}
For the boundary $\Gamma_d{=}\{I_0,I_1,\ldots,I_d\}{\in}\partial^d\mathcal{P}_{\{34\ldots n{-}1\}}$ where $I_d{\subset}I_{d-1}{\subset}\ldots{\subset}I_0{=}\{3,4,\ldots,n{-}1\}$ and $I_d{\neq}\emptyset$, the contribution is
\begin{align}
    K^{(2)}_{\mathrm{H},\Gamma_d}&= \frac{m v\cdot F_{\tau_{0}}\cdot v}{p_{23\ldots n{-}1}\cdot v} \prod_{k=1}^{d}\frac{p_{\Delta(I_{k},I_{k{+}1})}\cdot F_{\tau_d}\cdot v}{v\cdot p_{I_k}}\\\nonumber
    &=-\frac{2m^2 v\cdot F_{\tau_{0}}\cdot v}{k^2} \prod_{k=1}^{d}\frac{p_{\Delta(I_{k},I_{k{+}1})}\cdot F_{\tau_d}\cdot v}{v\cdot p_{I_k}},
\end{align}
where in the calculation of $\Delta(I_{k},I_{k{+}1})$, the complement set of $I_k$ is still taken to be $\{2,3,\ldots,n{-}1\}/I_k$. For $n=4$, there is no commutator in $\Theta^1$ and the result is \eqref{eq:Hamp4}.
For $n=5$, the amplitude becomes
\begin{equation}
\begin{aligned}
        A^\mathrm{H}(1,2,3,4,5){=} &\frac{1}{s_{12}} \frac{K_\mathrm{H}^{(2)}(1,2,[3,4],5)}{s_{34}}  \\
        & {+}\frac{1}{s_{123}} \frac{K_\mathrm{H}^{(2)}(1,[2,3],4,5)}{s_{23}},
\end{aligned}
\end{equation}
where $K_\mathrm{H}^{(2)}(1,2,3,4,5)$ is given by{\small
\begin{equation*}
 -\frac{2m^2}{k^2}(v \cdot F_{234} \cdot v+\frac{ v\cdot F_{24} \cdot v\ p_{2}\cdot F_3\cdot v }{v\cdot p_3 }
     + \frac{ v\cdot F_{23} \cdot v \ p_{23}\cdot F_4\cdot v }{v\cdot p_4}).
\end{equation*}}
Let us give a final example for $n=6$ amplitude
\begin{equation}
\begin{aligned}
    & \frac{1}{s_{12}} \left(\frac{K_\mathrm{H}^{(2)}(1,2,[[3,4],5],6)}{s_{34}s_{345}}{+} \frac{K_\mathrm{H}^{(2)}(1,2,[3,[4,5]],6)}{s_{45}s_{345}} \right) \\
    {+}& \frac{1}{s_{123}} \frac{K_\mathrm{H}^{(2)}(1,[2,3],[4,5],6)}{s_{23}s_{45}}  \\
    {+}&\frac{1}{s_{1234}} \left(\frac{K_\mathrm{H}^{(2)}(1,[[2,3],4],5,6)}{s_{23}s_{234}}{+} \frac{K_\mathrm{H}^{(2)}(1,[2,[3,4]],5,6)}{s_{34}s_{234}} \right).
\end{aligned}
\end{equation}
It is interesting to notice the numerators we present here only differ from those in \cite{Brandhuber:2021bsf} denoted by $N(1,2,\ldots,n)$ by an overall prefacto
\begin{equation}\label{eq: HBCJnumerator}
    K_\mathrm{H}^{(2)}(1,2,\ldots,n)=(-1)^{n}(n{-}2)\frac{2m}{k^2} v\cdot p_2 N(1,2,\ldots,n).
\end{equation}
It is highly nontrivial, however, that these two sets of effective BCJ numerators give the same HEFT amplitude. In~\cite{Brandhuber:2021bsf}, the expression involves the sum of cubic graphs corresponding to nested commutators of depth $n-3$ of the ordered set $(2,3,\ldots,n{-}1)$, thus the propagator denominator contains an overall factor $s_{23\ldots n-1}$, which in our case is replaced by different $s_{1 \sigma}$ for different terms. In addition, the numerator of~\cite{Brandhuber:2021bsf} for each cubic graph is given by a nested commutator of $N(1,2,\ldots,n)$, thus the number of terms in it is twice as ours. Nevertheless, we have analytically checked up to $n=10$ that the amplitude \eqref{eq: Hamp} agrees with \cite{Brandhuber:2021bsf}. Moreover, we have checked that although they look very different, the HEFT gravity amplitude via double copy also agrees with that in~\cite{Brandhuber:2021bsf}, and we expect both agreements to hold for all $n$.

\subsection{Decoupling into pure YM}
Given the explicit result of the $n$-point heavy mass BCJ numerators, the $(n-1)$-point pure YM BCJ numerators, as well as the amplitudes, can be easily obtained via the decoupling limit:
$mv\to \varepsilon_n$, $p_{23\ldots n-1}^2\to 0$
to obtain the BCJ numerator $K^{\prime \text{YM}}(2,3,\ldots,n)$ \cite{Brandhuber:2021kpo,Brandhuber:2021bsf}. Under this kinematics, the overall factor $k^2$ vanishes, which we ignore in the decoupling limit. For instance, the three-point BCJ numerator is given by $K^{\prime \text{YM}}(2,3,4)=-2\varepsilon_4 \cdot F_{23}\cdot \varepsilon_4$. Therefore, the three-point amplitude is
\begin{align}
    A^{\text{YM}}(2,3,4)&=  -\frac{\varepsilon_4 \cdot F_{23}\cdot \varepsilon_4}{ \varepsilon_4 \cdot p_2}\\\nonumber
    &=\varepsilon _4\cdot \varepsilon _2 p_2\cdot \varepsilon _3{-}\varepsilon _2\cdot \varepsilon _3 p_2\cdot \varepsilon _4{-}\varepsilon _4\cdot \varepsilon _3 p_3\cdot \varepsilon _2.
\end{align}
For the 4-point YM amplitude, the numerator $K^{\prime \text{YM}}(2,3,4,5)$ reads
\begin{equation*}
    {-}2\left(\varepsilon_{5}{\cdot}F_{234}{\cdot} \varepsilon_5{+}\frac{\varepsilon_5{\cdot} F_{24}{\cdot} \varepsilon_5 p_{2}{\cdot} F_{3}{\cdot} \varepsilon_5}{p_3{\cdot}\varepsilon_5 }{+}\frac{\varepsilon_5{\cdot} F_{23}{\cdot} \varepsilon_5 p_{23}{\cdot} F_4{\cdot} \varepsilon_5 }{ p_4{\cdot}\varepsilon_5}\right).
\end{equation*}
Note that $K^{\prime \text{YM}}(2,3,\ldots,n)$ also manifests the gauge invariance of particles $2,3,\ldots,n-1$. Moreover, it is related to the BCJ numerator given in sec.\ref{sec:YM} via
\begin{equation*}
    K^{\prime \text{YM}}(2,3,\ldots,n)=2\left. \varepsilon_n \cdot p_2 K^{\text{YM}}(2,3,\ldots,n)\right|_{q_I\to \varepsilon_n}.
\end{equation*}
These numerators, accompanied by different $\phi^3$ amplitudes, produce the same YM amplitude.

\section{Conclusions and outlook}

In this note, we established a correspondence between BCJ numerators from covariant color-kinematics duality and the combinatorial permutohedra, which are closely related to the quasi-shuffle Hopf algebra. This apply to all YMS amplitudes, but the most interesting case is that with two scalars, whose numerators share the same combinatorial structure as the pure YM ones: each term is mapped to a boundary of $\mathcal{P}_{n-2}$; the contribution from each boundary is almost identical in these two cases, except that we need to modify one factor to take into account the remaining two gluons. We also found nice recursion relations and factorization properties implied by this picture. Finally, based on highly nontrivial cancellations which are needed for both YMS and gravity amplitudes (via double copy) in HEFT, we conjectured a compact formula for their effective numerators; they become closely related to permutohedra $\mathcal{P}_{n-3}$, which, while producing the same amplitude, differ by an overall factor from the numerators in ~\cite{Brandhuber:2021bsf, Brandhuber:2022enp}.

There are numerous open questions for further investigations. First, as we will present in \cite{toapp}, it is interesting to see how lower-dimensional permutohedra for general YMS numerators combine into $\mathcal{P}_{n-2}$ which corresponds to the pure YM ones; we also find interesting combinatorial structures underlying BCJ numerators of amplitudes in NLSM {\it etc.}. Moreover, the somewhat miraculous cancellations that simplify these numerators in HEFT still remain to be proven, which would also be important to establish the correct double copy in HEFT. Since the final amplitudes are independent of reference momenta, all the spurious poles must cancel, which still calls for a direct understanding (without relying on the CCK duality); such an understanding could connect this combinatorial picture (especially the factorizations) to the uniqueness theorem for YM amplitude~\cite{Arkani-Hamed:2016rak, Rodina:2016jyz} and YMS ones via the universal expansion~\cite{Dong:2021qai}. Last but not least, it is tempting to ask: could we combine the permutohedra for BCJ numerators with the associahedra for bi-adjoint $\phi^3$ amplitudes, and obtain a unified geometric understanding of gluon and graviton scattering?

\begin{acknowledgments}
We thank Linghui Hou, Guanda Lin, Tianheng Wang for discussions and collaborations on related projects. This research is supported in part by the Key Research Program of CAS, Grant No. XDPB15 and National
Natural Science Foundation of China under Grant No.
11935013,11947301,12047502,12047503
\end{acknowledgments}

\newpage

\widetext

\section{Introduction}
Despite very different natures, gauge theories and gravity have deep connections; one of the oldest and the most prominent example is the double copy structure \cite{Kawai:1985xq, Bern:2008qj, Bern:2010ue}. Originally it was discovered from Kawai-Lewellen-Tye (KLT) relations \cite{Kawai:1985xq} in string theory, and a modern realization of double copy has relied on the duality between color and kinematics for gauge theory amplitudes, where the Bern-Carrasco-Johansson (BCJ) kinematic numerators satisfy the same Jacobi relations as the color factors. The duality and double copy have led to tremendous progress in the study of amplitudes both in gauge theory and gravity (see~\cite{Bern:2019prr, Bern:2022wqg, Adamo:2022dcm} and references therein). More recently, the authors of~\cite{Cheung:2021zvb} have revealed a so-called covariant color-kinematics (CCK) duality for a large class of theories including Yang-Mills theory (YM) and its coupling to bi-adjoint $\phi^3$ (YMS). As a consequence, the duality implies new, closed-form expression for BCJ numerators of all tree-level amplitudes in YMS and YM theory.
Previous works on BCJ numerators and kinematic algebras include~\cite{Mafra:2011kj,Bargheer:2012gv,CHY3,He:2015wgf,Fu:2017uzt,Teng:2017tbo, Du:2017kpo, He:2017spx, Edison:2020ehu,He:2021lro, Monteiro:2011pc,Monteiro:2013rya,Cheung:2016prv,Chen:2019ywi,Edison:2022jln} and references therein.

On the other hand, recent years have seen progress on revealing new geometric/combinatorial structures underlying scattering amplitudes {\it e.g.} from the (all-loop) amplituhedron of supersymmetric Yang-Mills ~\cite{Arkani-Hamed:2013jha} to the associahedron for bi-adjoint $\phi^3$ at tree level~\cite{Arkani-Hamed:2017mur} (with extensions to string scattering~\cite{Arkani-Hamed:2019mrd}). It is natural to look for hints of such structures underlying YM and gravity amplitudes; instead of directly working with tree amplitudes, one may decompose the problem and ask a somewhat strange question as a first step: are there combinatorial structures underlying BCJ numerators?

In this note, we take BCJ numerators from CCK duality~\cite{Cheung:2021zvb} as inputs and present preliminary evidence for such structures: in addition to the more familiar quasi-shuffle Hopf algebras~\cite{Hoffman}, we find hidden combinatorial permutohedra \cite{wiki:Permutohedron} for BCJ numerators. Any BCJ numerator can be written as the sum over all boundaries of a permutohedron (or terms from a quasi-shuffle product); for a co-dimension $d$ boundary (length-$d$ term), it contains a product of $d{+}1$ factors each with a spurious pole and a gauge-invariant numerator. We will focus on the case with two scalars and $n{-}2$ gluons, which corresponds to a $(n{-}3)$-dimensional permutohedron, and it has Fubini number ${\cal F}_{n{-}2}$ boundaries with co-dimensions $d=0, 1, \ldots, n{-}3$; each boundary is labeled by $d+1$ subsets, and for each factor labeled by such a set both the numerator (which is gauge invariant in the gluons) and the (spurious pole) denominator are given by Lorentz products of momenta and polarizations, as well as the reference momenta. Apart from being the most illustrative BCJ numerators of YMS cases, we will also see that they give nice BCJ numerators in the heavy-mass effective theory (HEFT)~\cite{Georgi:1990um,Luke:1992cs,Neubert:1993mb,manohar_wise_2000,Damgaard:2019lfh,Brandhuber:2021kpo, Brandhuber:2021bsf, Brandhuber:2022enp} as well as decoupling limit into pure YM amplitudes. BCJ numerators in HEFT have attracted lots of interest recently for their roles in the computation of gravitational amplitudes for black-hole scattering and gravitational waves~\cite{Brandhuber:2021eyq} ({\it c.f.} ~\cite{Kosower:2018adc,Bern:2019nnu,Damour:2019lcq,Bern:2021dqo,DiVecchia:2021bdo,Herrmann:2021tct,Bjerrum-Bohr:2021din,Bjerrum-Bohr:2021wwt,Jakobsen:2021lvp} for some recent works). We will take the heavy-mass limit of YMs amplitude, and (as we have checked up to $n=10$) highly nontrivial cancellations lead to a nice formula for BCJ numerators in HEFT which corresponds to ${\cal P}_{n{-}3}$ (one dimensional lower)!

Furthermore, our results imply new recursion relations and surprisingly, ``factorization'' properties of BCJ numerators on facets of permutohedra; all these can be extended to BCJ numerators of general YMS amplitudes, which in turn combine into a formula for the YM case as well. For the latter, we can then turn the logic around: since the BCJ numerators are manifestly gauge invariant in $n{-}1$ gluons, by showing that all spurious poles indeed cancel in the amplitude based on such ``factorizations", it follows from the uniqueness theorem of~\cite{Arkani-Hamed:2016rak} that they must give correct YM and gravity amplitude (after double copy) even without knowing the CCK duality.

Let us consider color-ordered YMS amplitude $A(1^\phi,2,\ldots,n{-}1,n^\phi)$ with scalars $1^\phi,n^\phi$. Its expansion onto the Kleiss-Kuij(KK) basis~\cite{Kleiss:1988ne} of bi-adjoint $\phi^3$ amplitudes has  BCJ master numerators as coefficients reads
\begin{equation}\label{eq: BCJnum}
 A(1^\phi,2,\ldots,n{-}1,n^\phi){=}\sum_{\beta\in S_{n-2}} K(1,\beta,n)A^{\phi^3}(1,\beta,n),
\end{equation}
where the sum is over $(n{-}2)!$ permutations of gluons and $A^{\phi^3}(1,\beta,n) \equiv m(1,2,\ldots,n | 1,\beta,n)$ denotes bi-adjoint $\phi^3$ amplitudes with the first ordering fixed to be $(1,2,\ldots,n)$.

Remarkably, the BCJ numerators from CCK duality $K(1,\beta,n)$ respect the Bose symmetry of all the $n{-}2$ gluons ~\cite{Cheung:2021zvb}: we only need a single numerator with the ordering chosen to be $\beta=(2,\ldots,n{-}1)$, and all others can be obtained by relabelling; they are also gauge invariant for the gluons, which becomes manifest since the dependence on polarizations is through Lorentz products of linearized field strengths $F_i^{\mu\nu}\equiv p_i^\mu \varepsilon_i^\nu-p_i^\nu\varepsilon_i^\mu$
\begin{equation}\label{Fprod}
[F_{\sigma}]^{\mu \nu}= [F_{\sigma_{1}}\cdot F_{\sigma_{2}}\cdots \cdot F_{\sigma_{|\sigma|}}]^{\mu \nu}
\end{equation}
for an ordered subset $\sigma$. The price to pay for these desirable properties is the presence of $2^{n{-}2}-1$ spurious poles, one for each nonempty subset $I\subset \{2, \ldots, n{-}1\}$:
\begin{equation}\label{Ddef}
  D_{I}:= p_{I}\cdot q_{I},\quad  \text{with} \quad p_I:=\sum_{i\in I}p_i,
\end{equation}
which depends on a reference momentum  $q_{I}$. These numerators can be simplified with some choices of $q_I$, and the final amplitude is independent of them.

\section{The permutohedron and algebra underlying BCJ numerators}
In this section, we show that all the terms in a BCJ master numerator obtained from CCK duality for YMS amplitudes are in one-to-one correspondence with all boundaries of permutohedron $\mathcal{P}_{n{-}2}$, or equivalently terms from a quasi-shuffle product.

\subsection{The (combinatorial) permutohedra and quasi-shuffle products}\label{sec:def}
Following~\cite{Cheung:2021zvb}, we organize $K(1,2,\ldots, n)$ according to the spurious pole structure, which is isomorphic to the boundary structure of the permutohedron $\mathcal{P}_{n-2}$.

The permutohedron $\mathcal{P}_{n-2}$ is an $(n{-}3)$-dimensional polytope~\cite{wiki:Permutohedron}, whose co-dimension $d$ boundary $\Gamma_d$ can be labeled by $d{+}1$ consecutive subsets
\begin{equation}\label{eq: lBD}
    \Gamma_d:=\{I_{0}, I_{1},\ldots, I_{d}\},
\end{equation}
where $I_d{\neq}\emptyset$ and $I_d\subset I_{d{-}1}\subset\ldots\subset I_0=\{2,3,\ldots,n{-}1\}$; the interior of ${\cal P}_{n{-}2}$ can be viewed as its co-dimension $0$ boundary, $\Gamma_0:=I_0$. ${\cal P}_{n{-}2}$ and its boundaries have appeared in the context of cubic tree graphs from the worldsheet~\cite{Gao:2017dek, He:2018pue}. Here each term of the BCJ numerator $K(1, 2, \ldots, n)$ with $d{+}1$ spurious poles corresponds to such a co-dimension $d$ boundary, thus the numerator can be expanded in terms of boundaries of $\mathcal{P}_{n-2}$
\begin{equation}\label{eq: KBDexp}
K(1,2,\ldots,n)=\sum_{d=0}^{n{-}3}\sum_{\Gamma_d\in\partial^d\mathcal{P}_{n{-}2}}K_{\Gamma_d}(1,2,\ldots,n),
\end{equation}
where we sum over all boundaries $\Gamma_d\in\partial^d\mathcal{P}_{n{-}2}$ with co-dimension $d=0, \ldots, n{-}3$, and the contribution from $\Gamma_d$, $K_{\Gamma_d}(1,2,\ldots,n)\equiv K_{\Gamma_d}$ reads
\begin{equation}\label{eq: KBD}
K_{\Gamma_d}=\prod_{k=0}^{d}\frac{p_{1
\Delta(I_{k},I_{k+1})
}\cdot F_{\tau_{k}}\cdot q_{I_{k}}}{D_{I_{k}}}.
\end{equation}
It has $d{+}1$ factors each with a denominator $D_{I_k}$ of \eqref{Ddef} and a numerator of the form $p_{1 \Delta(I_{k},I_{k+1})}\cdot F_{\tau_k} \cdot q_{I_k}$ for $k=0, \dots, d$. To specify the ordered subset $\tau_k$ of the Lorentz product as in \eqref{Fprod}, we introduce an alternative form of \eqref{eq: lBD} using ordered sets:
\begin{align}\label{eq: lBDt}
   \Gamma_d
   =&\{\tau_0\cup\tau_1\cup\ldots\cup\tau_d,\tau_1\cup\ldots\cup\tau_d,\ldots,\tau_{d{-}1}\cup\tau_d,\tau_d\}\\\nonumber
   \sim&\{\tau_0,\tau_1,\ldots,\tau_d\},
\end{align
where the first line is equivalent to \eqref{eq: lBD} but we use ordered sets $\tau_{k}=\mathrm{Id}(I_k/I_{k{+}1})$ with $I_{d{+}1}\equiv \emptyset$. $\mathrm{Id}(I)$ means sorting the subset $I$ in numerical ordering~\footnote{For $K(1,\beta,n)$, everything stays the same except that in \eqref{eq: lBDt} the definition of $\tau_i$'s become $\tau_{k}=\beta(I_k/I_{k{+}1})$, {\it i.e.} any subset $I$ is sorted according to $\beta$ ordering.}. Perhaps the most subtle point is that we also define
$\Delta(I_{k},I_{k+1}):=\left(\bar{I}_{k}\right)_{<\tau_{{k},1}}$, which refers to the elements in the set $\bar{I}_{k}=\{2,3,\ldots,n-1\}/I_k$ that are numerically smaller than the first element $\tau_{k,1}$ of the ordered set $\tau_{k}$.

For example, at $n=5$, we have a boundary $\Gamma_2=\{I_0=\{2,3,4\}, I_1=\{2,4\}, I_2=\{4\}\}$; equivalently, we have $\tau_0=\{3\}, \tau_1=\{2\}, \tau_2=\{4\}$, thus we have a term
\begin{equation}
K_{234,24,4}= \frac{p_1\cdot F_3\cdot q_{234}  p_1\cdot F_2\cdot q_{24} p_{123}\cdot F_4\cdot q_4 }{D_{234}D_{24} D_4  }
\end{equation}
where we have used $\Delta(I_2,\emptyset)=\left.\left(\bar{I_2}\right)\right|_{<4}=\{2,3\}$. A more nontrivial example for the latter is for $n{=}9$, $\Delta(\{4,5,7,8\},\{4,8\})=\left.\{2,3,6\}\right|_{<5}=\{2,3\}$.

On the other hand, the boundaries of $\mathcal{P}_ {n-2}$ have a nice quasi-shuffle product interpretation. The quasi-shuffle product $\star$ can be defined between two arbitrary generators $\left(\sigma_0,\sigma_1,\ldots,\sigma_r\right)$ and $\left(\rho_0,\rho_1,\ldots,\rho_s\right)$, where $\sigma_i$ and $\rho_j$ are sets with arbitrary lengths and can also be prompted to the quasi-shuffle Hopf algebra~\cite{Hoffman, Brandhuber:2021bsf, Brandhuber:2022enp, Chen:2022nei}. We summarize the definitions in appendix \ref{sec:appAG} and here we just use the following result of the quasi-shuffle product $\hat{K}(2,\ldots, n{-}1)\equiv(2)\star(3)\star\ldots\star(n-1)$
\begin{equation}\label{eq: part}
    \hat{K}(2,\ldots,n{-}1)
    {=}\sum_{d{=}0}^{n{-}3}\sum_{\tau\in \text{part}^{(d{+}1)}(2,\ldots,n{-}1)}\!\!\!\!\!\!\!\!\!\!\!\!\!\!(-1)^{n+d-1}(\tau_0,\tau_1,\ldots,\tau_{d}),
\end{equation}
where $\text{part}^{(d+1)}(2,\ldots,n{-}1)$ denotes all the ordered partitions of $\{2,3,\ldots,n{-}1\}$ into $d+1$ nonempty subsets $(\tau_0,\tau_1,\ldots,\tau_d)$ (each $\tau_i$ is sorted according to $\beta$).  Terms on the RHS of \eqref{eq: part} are in one-to-one correspondence with boundaries of $\mathcal{P}_{n-2}$ as in~\eqref{eq: lBDt}, thus we can rewrite \eqref{eq: KBDexp} in terms of the quasi-shuffle product
\begin{equation}
\label{eq:qsp}
K(1,2,\ldots,n)=\langle\hat{K}(2,\ldots n{-}1)\rangle,
\end{equation}
where we have defined a {\it linear map} $\langle\cdot\rangle$ from \eqref{eq: KBD} for any partition $(\tau_0,\tau_1,\ldots,\tau_{d})$
\begin{equation}\label{eq:linmap}
    \langle(\tau_0,\tau_1,\ldots,\tau_{d})\rangle{=}(-1)^{n+d-1}\prod_{k=0}^{d}\frac{p_{1\Delta(I_{k},I_{k+1})
}\cdot F_{\tau_{k}}\cdot q_{I_{k}}}{D_{I_{k}}}.
\end{equation}

\subsection{The counting and some examples}
By definition, the permutohedron $\mathcal{P}_m$ contains $m!$ vertices and $2^{m}-2$ co-dimension one facets. More generally, the number of co-dimension $d$ boundaries of this polytope is $(d+1)!S(m,d+1)$, where $S(m,d)$ is the second kind of Stirling number \cite{wiki:Stirling_numbers_of_the_second_kind}.  Algebraically, $S(m,d+1)$ also counts the number of ways to partition a set of $m$ labeled objects into $d+1$ nonempty unlabeled subsets $\{\tau_0,\tau_1,\ldots,\tau_d\}$ \cite{wiki:Stirling_numbers_of_the_second_kind}, so after considering the ordering between these sets, there are $(d+1)!S(m,d+1)$ terms in the summation for any $d$.
The total number is the Fubini number $\mathcal{F}_{m}$, where $\mathcal{F}_{m}=\sum_{d=1}^{m}d!S(m,d)$ \cite{Mezo}, thus the $n$-point BCJ numerator has $\mathcal{F}_{n-2}$ terms.
\begin{table}[H]
  \centering
    \begin{tabular}{|c|c|c|c|c|c|}
    \hline
    \diagbox[innerwidth=0.5cm]{n}{d}      & 0     & 1     & 2     & 3 &Total \\
    \hline
    3     & 1 &       &       &  &1\\
    \hline
    4     & 1 & 2 &       &  &3\\
    \hline
    5     & 1 & 6 & 6 & &13 \\
    \hline
    6     & 1 & 14 & 36 & 24&75 \\
    \hline
    \end{tabular}%
  \label{tab:BDc}%
  \caption{Counting co-dimension-$d$ boundaries of $\mathcal{P}_{n-2}$}
\end{table}%

Let us illustrate  \eqref{eq: KBDexp} and \eqref{eq: KBD} with some examples. The most trivial case is $n=3$, where the BCJ numerator corresponds to the zero-dimensional permutohedron $\mathcal{P}_{1}$ which is just a point. It contains one term with $\Gamma_0=\{I_0\}=\{2\}$, thus $K(1,2,3)=K_{2}(1,2,3)=\frac{p_1\cdot F_{2}\cdot q_2}{D_2}$, where we have used $\Delta(I_{0},I_{1})=\emptyset$, this is generally true since $\bar{I_0}=\emptyset$.

For $n=4$, the permutohedron $\mathcal{P}_{2}$ is a line segment, where the interior ($d=0$) is labelled by $I_0=\{23\}$, and the two vertices ($d=1$) are labeled by $\{23,2\}$ and $\{23,3\}$; we show these three terms in figure \ref{fig:P2}.

\begin{figure}[H]
    \centering
    \subfloat[]{ \label{fig:P2}
    \begin{minipage}[c]{1\linewidth}
    \centering
    \includegraphics[scale=1.5]{figs/P2K.pdf}
    \end{minipage}
    }

    \subfloat[]{ \label{fig:P3}
    \begin{minipage}[c]{1\linewidth}
    \centering
    \includegraphics[scale=1.3]{figs/P3D.pdf}
    \end{minipage}
    }
    \caption{Permutohedra ${\cal P}_2$ for $K(1,2,3,4)$(top) and ${\cal P}_3$ for  $K(1,2,3,4,5)$(bottom)}
\end{figure}

Equivalently, in \eqref{eq: part} the partition $\text{part}^{(1)}$ of $\{2,3\}$ has $(\{\tau_0=\{2, 3\})$ and $\text{part}^{(2)}$ has  $(\tau_0=\{2\},\tau_1=\{3\})$ and $(\tau_0=\{3\},\tau_1=\{2\})$: they are nothing but the interior and the two vertices, according to \eqref{eq: lBDt}.

Thus the BCJ numerator $K(1,2,3,4)$ has three terms, $K_{23}$, $K_{23,2}$ and $K_{23,3}$, which read
\begin{align} \label{eq: BCJnum4}
\frac{p_{1}{\cdot}F_{23}{\cdot} q_{23}}{D_{23}}{+}\frac{p_{1}{\cdot} F_{3}{\cdot} q_{23} p_{1}{\cdot} F_{2}{\cdot} q_{2}}{D_{23}D_{2}}{+}\frac{p_{1}{\cdot} F_{2}{\cdot} q_{23} p_{12}{\cdot} F_{3}{\cdot} q_{3} }{D_{23}D_{3}
\end{align}
Notice that the last term in the above equation is from the boundary $\{23,3\}$, so the second factor in the numerator is $p_{1\left(\overline{\{3\}}\right)_{<3}}\cdot F_{3}\cdot q_{3}=p_{12}\cdot F_{3}\cdot q_{3}$. Meanwhile \eqref{eq: BCJnum4} shows that the four-point numerator $K(1,2,3,4)$ has an overall pole $D_{I_0}=D_{23}$. This can be easily seen from \eqref{eq: lBD} since the first set of any co-dimension $d$ boundary $\Gamma_d$ is always labeled by $I_0=\{2,3,\ldots,n{-}1\}$.

Notice that \eqref{eq: KBD} means each term in the BCJ numerator contains spurious poles, and the co-dimension $d$ contribution will have $d{+}1$ spurious poles where $D_{I_0}$ is an overall pole for every boundary.
Except for the overall one, the simple poles can be written as $D_I$ where $I$ is a nonempty proper subset of $I_0=\{2,3,\ldots,n{-}1\}$ and two simple poles $D_I$ and $D_J$ are compatible if and only if $I\subset J$ or $J\subset I$. For example, at five-point, except for the overall $D_{234}$, the simple poles are $D_{2},D_{3},D_{4},D_{23},D_{24},D_{34}$ which correspond to the six co-dimension one boundaries of $\mathcal{P}_3$, and the compatible double poles are
\begin{equation}
    \{D_{23}D_{2},D_{24}D_{2},D_{23}D_{3},D_{34}D_{3},D_{24}D_{4},D_{34}D_{4}\},
\end{equation}
which correspond to six vertices of the hexagon $\mathcal{P}_3$. We show the boundary contribution formally in figure \ref{fig:P3}.
These $13$ terms form a two-dimensional polytope $\mathcal{P}_{3}$. These boundaries can also be realized in quasi-shuffle product $\hat{K}(2,3,4)$, which will be discussed in appendix \ref{sec: N5}. To be precise, we give some explicit examples of different co-dimension here
\begin{equation}
\begin{aligned}
    &K_{234}{=}\frac{p_1{\cdot}F_{234}\cdot q_{234}}{D_{234}},\\
    &K_{234,23}{=}\frac{p_1{\cdot}F_4{\cdot}q_{234} p_1{\cdot} F_{23}{\cdot}q_{23}}{D_{234}D_{23}},\;\,K_{234,2}{=}\frac{p_1{\cdot} F_{34}{\cdot} q_{234} p_1{\cdot}F_2{\cdot}q_2}{D_{234}D_2}\\
    &K_{234,23,2}{=}\frac{ p_1{\cdot}F_4{\cdot}q_{234} p_1{\cdot} F_3{\cdot} q_{23} p_1{\cdot} F_2{\cdot} q_2}{D_{234}D_{23}D_2 }.
\end{aligned}
\end{equation}
The complete result for the BCJ numerator $K(1,2,3,4,5)$ is shown in the appendix \ref{sec: N5}.

Moreover, we emphasize that all the spurious poles are canceled in the final amplitude, and the amplitude does not depend on the reference momenta. The proof will be put into the following paper \cite{toapp}.

For $n=6$, $\mathcal{P}_4$ is a three-dimensional truncated octahedron shown in figure \ref{fig: P4}. As we have counted, it contains $14$ co-dimension one boundaries (six squares and eight hexagons), $36$ edges and $4!=24$ vertices, thus $75$ boundaries in total. Some terms with co-dimension $d=0,1,2,3$ are
\begin{align}
    &K_{2345}=\frac{p_1{\cdot} F_{2345}{\cdot} q_{2345}}{D_{2345}}\\\nonumber
    &K_{2345,234}=\frac{p_1{\cdot} F_{5}{\cdot} q_{2345}p_1{\cdot} F_{234}{\cdot} q_{234}}{D_{2345}D_{234}}\\\nonumber
    &K_{2345,234,23}=\frac{p_1{\cdot} F_{5}{\cdot} q_{2345}p_1{\cdot} F_{4}{\cdot} q_{234}p_1{\cdot} F_{23}{\cdot} q_{23}}{D_{2345}D_{234}D_{23}}\\\nonumber
    &K_{2345,234,23,2}=\frac{p_1{\cdot} F_{5}{\cdot} q_{2345}p_1{\cdot} F_{4}{\cdot} q_{234}p_1{\cdot} F_{3}{\cdot} q_{23}p_1{\cdot} F_2{\cdot} q_2}{D_{2345}D_{234}D_{23}D_2}
\end{align}
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.4]{figs/P4label.pdf}
    \caption{The permutohedron $\mathcal{P}_4$ for $K(1,\ldots, 6)$}
    \label{fig: P4}
\end{figure}
\section{Numerators for general YMS and pure YM amplitudes}
\label{sec:YM}
More generally, the CCK duality has provided closed-formulae for BCJ numerators of $n$-point YMS amplitude with $r\geq 2$ scalars~\cite{Cheung:2021zvb}. It turns out that any such numerator corresponds to a permutohedron ${\cal P}_{n{-}r}$ (with dimension $n{-}r{-}1$): everything we have discussed above for $r=2$ case still applies if we replace $\{2, 3,\ldots, n{-}1\}$ by the set of $n{-}r$ gluons. For example, in the other extreme with $r=n$ we can formally define ${\cal P}_0:={\cal P}_{\emptyset}$ and the $n$-scalar numerator (of bi-adjoint $\phi^3$ amplitude) is $1$ or $0$. For $r=n{-}1$, we have the zero-dimensional ${\cal P}_1$ and the numerator with a single gluon $i$ reads
\begin{equation}
K^{1-{\rm gluon}}(1, \beta, n)=\frac{p_{1{\rightarrow }i} \cdot F_i \cdot q_i}{D_i}
\end{equation} where $1 {\rightarrow} i$ denotes all the scalars preceding $i$ in $(1\beta n)$. We shall not repeat this for general cases but leave detailed discussions to a separate paper~\cite{toapp}.

Moreover, since the pure YM amplitude can be expanded as a linear combination of these YMS ones~\cite{Lam:2016tlk, Fu:2017uzt, Du:2017kpo, Cheung:2017ems, Dong:2021qai}, we obtain its BCJ numerators for free; the resulting numerator naively contains $2\mathcal{F}_{n-2}$ terms as derived in \cite{Cheung:2021zvb}. However, we can still organize the terms according to pole structures and immediately combine them in pairs as $\mathcal{F}_{n-2}$ terms: the resulting numerator has the same form as the $2{-}$scalar case and corresponds to boundaries of permutohedron $\mathcal{P}_{n-2}$. By expanding $A^{\text{YM}}(1,2,\ldots,n)$ in exactly the same way as \eqref{eq: BCJnum}
each master BCJ numerator, {\it e.g.}  $K^{\text{YM}}(1,2,\ldots,n)$ is given by a sum over boundaries of $\mathcal{P}_{n-2}$ as in \eqref{eq: KBDexp}:
\begin{align}
    K^{\text{YM}}(1,2,\ldots,n)&=\sum_{d=0}^{n{-}3}\sum_{\Gamma_d\in\partial^d\mathcal{P}_{n{-}2}}K^{\text{YM}}_{\Gamma_d}(1,2,\ldots,n). \nonumber
\end{align}
where the contribution from each boundary
is identical to \eqref{eq: KBD} except for the $k=0$ factor which becomes
\begin{equation}\label{eq: N0YM}
    \frac{\varepsilon_n\cdot F_{1\tau_0}\cdot q_{I_{0}}+\varepsilon_1\cdot F_{\tau_0}\cdot\left(\varepsilon_n p_{1n}\cdot q_{I_{0}}-q_{I_{0}} p_{1}\cdot\varepsilon_n\right)}{D_{I_0}}
\end{equation}

Of course, similar to the YMS case, all spurious poles cancel in the final amplitude, which does not depend on $q_I$. Therefore we are free to choose them to simplify the expression \eqref{eq: N0YM}. One such choice is $q_{I_{0}}=\varepsilon_n$, and the $k=0$ factor \eqref{eq: N0YM} takes a simpler form
\begin{equation}\label{eq: KYMSim}
    \frac{\varepsilon_n{\cdot} F_{1\tau_0}{\cdot} \varepsilon_n}{\varepsilon_n{\cdot} p_{23\ldots n{-}1}}=-\frac{\varepsilon_n{\cdot} F_{1\tau_0}{\cdot} \varepsilon_n}{\varepsilon_n{\cdot} p_{1}}.
\end{equation}
It is easy to see that the BCJ numerators become manifestly gauge invariant in particles $1,2,\ldots,n{-}1$.

For example, the BCJ numerator $K^{\text{YM}}(1,2,3,4)$ reads
\begin{equation*}
    {-}\frac{\varepsilon_{4}{\cdot}F_{123}{\cdot} \varepsilon_4}{p_1{\cdot}\varepsilon_4}{-}\frac{\varepsilon_4{\cdot} F_{13}{\cdot} \varepsilon_4 p_{1}{\cdot} F_{2}{\cdot} q_{2}}{p_1{\cdot}\varepsilon_4 D_{2}}{-}\frac{\varepsilon_4{\cdot} F_{12}{\cdot} \varepsilon_4 p_{12}{\cdot} F_{3}{\cdot} q_{3} }{p_1{\cdot}\varepsilon_4 D_{3}}.
\end{equation*}

Furthermore, similar to the discussion in section \ref{sec:def}, BCJ numerators of YM amplitudes can also be interpreted in terms of quasi-shuffle products, and the only change is that in the linear map \eqref{eq:linmap} the $k=0$ factor is modified to  \eqref{eq: N0YM}.

Before ending the section, we mention the obvious double-copy from YM to GR
\begin{equation}
M^{\rm GR}_n=\sum_{\alpha, \beta} K^{\rm YM}(1, \alpha, n) m(1,\alpha ,n|1, \beta ,n) K^{\rm YM}(1, \beta, n)
\end{equation}
where we sum over a pair of permutations $\alpha,\beta$ of $\{2,3,\ldots, n{-}1\}$, with $m$ denoting bi-adjoint $\phi^3$ amplitudes; if we replace YM by YMS with $1, n$ being scalars, it gives the amplitude with $n{-}2$ gravitons and two scalars.

\section{Recursions and factorizations}
In this section, we propose recursion relations and factorization properties (on spurious poles $D_I$) for the BCJ numerators, which are implied by the combinatorial and algebraic structure. The argument can be equally applied to both two-scalar YMS and pure YM numerators.

\subsection{Recursion relations}
First, in quasi-shuffle product~\eqref{eq: part}, one can collect the terms with the same $\tau_d$ and then apply the linear map \eqref{eq:linmap} to obtain the following recursion relation,
\begin{equation}\label{eq: Krec}
    K(1,2,\ldots,n){=
    \sum_{I\subset\{2,\ldots,n{-}1\}}\frac{p_{1\Delta(I,\emptyset
    }{\cdot} F_I {\cdot} q_I}{D_I}\tilde{K}^I(1,\bar{I},n),
\end{equation}
where the summation is over all the nonempty subsets of $\{2,3,\ldots,n{-}1\}$. The definition of $\tilde{K}^I(1,\bar{I},n)$ is slightly different from \eqref{eq: KBD} in the denominator: it is given by the sum over boundaries of the permutohedron $\mathcal{P}_{\bar{I}}$ with vertices labeled by all permutations of set $\bar{I}$, and for each boundary $\Gamma_d=\{J_0=\bar{I},J_1,\ldots,J_{d}\}$ where $J_d\subset J_{d{-}1}\ldots\subset J_0$, we have a  contribution
\begin{equation}\label{eq: tKBD}
\tilde{K}^I_{\Gamma_d}=\prod_{k=0}^{d}\frac{p_{1\Delta(J_{k},J_{k+1})
}\cdot F_{\tau_k}\cdot q_{IJ_{k}}}{D_{IJ_{k}}},
\end{equation}
where $D_{IJ_k}=p_{IJ_k}\cdot q_{IJ_k}$, $\tau_k= \mathrm{Id}(J_k/J_{k{+}1})$ with $J_{d{+}1}\equiv \emptyset$ and the complement of the set $J_k$ appears in $\Delta(J_{k},J_{k+1})$ is defined as $\bar{I}/J_k$. For $\abs{I}{=}n{-}2$ ($\bar{I}=\emptyset$), we define $\tilde{K}^I(1,n)=1$. Formally, this numerator corresponds to the permutohedron $\mathcal{P}_0$.

For example, the recursion relation of $K(2,3,4)\equiv K(1,2,3,4,5)$ \footnote{Here we have omitted the labels of the scalar particles $1$ and $n$ in the numerators $K$ and $\tilde{K}$.} reads
\begin{small}
\begin{align}\label{eq: rec5}
    K&(2,3,4)=\frac{p_1{{\cdot}} F_{234}{\cdot} q_{234}}{D_{234}}\\\nonumber
    {+}&\frac{p_1{\cdot} F_{23}{\cdot} q_{23}}{D_{23}}\tilde{K}^{23}(4){+}\frac{p_{1}{\cdot} F_{24}{\cdot} q_{24}}{D_{24}}\tilde{K}^{24}(3){+}\frac{p_{12}{\cdot} F_{34}{\cdot} q_{34}}{D_{34}}\tilde{K}^{34}(2)\\\nonumber
    {+}&\frac{p_1{\cdot} F_2{\cdot} q_2}{D_2}\tilde{K}^2(3,4)
    {+}\frac{p_{12}{\cdot} F_3{\cdot} q_3}{D_3}\tilde{K}^3(2,4){+}\frac{p_{123}{\cdot} F_4{\cdot} q_4}{D_4}\tilde{K}^4(2,3).
\end{align}
\end{small}
Geometrically, the recursion relation~\eqref{eq: Krec} tells us how the co-dimension one boundaries of permutohedron are glued together. In the above five-point example, the term with $\abs{I}=3$ in the first line has only one pole and corresponds to the interior (co-dimension $0$ boundary) of $\mathcal{P}_{3}$, depicted in figure~\ref{fig:rec5}. For the three terms with $\abs{I}=2$ in the second line, each factor $\tilde{K}^I(1,\bar{I},n)$ corresponds to a zero-dimensional permutohedron; on the other hand, each term is mapped to a co-dimension one boundary of $\mathcal{P}_3$ without vertices. For the remaining three terms with $\abs{I}=1$ in the last line, each $\tilde{K}^I(1,\bar{I},n)$ corresponds to a one-dimensional permutohedron and it is mapped to a co-dimension one boundary with two vertices.
\begin{figure}[H]
    \centering
    \includegraphics[scale=0.8]{figs/rec5.pdf}
    \caption{Recursion relation at $n=5$}
    \label{fig:rec5}
\end{figure}
\subsection{Factorization properties on spurious poles}
Next, we move to certain intriguing factorization properties of the BCJ numerator on spurious poles. Combinatorially, any co-dimension one boundary of the permutohedron $\mathcal{P}_{n{-}2}$ is the product of two lower-dimensional permutohedra $\mathcal{P}_{I}\times\mathcal{P}_{\bar{I}}$. Remarkably, we find that on any pole $D_I=0$, the residue of the BCJ numerator factorizes into the product of a $(\abs{I}{+}2)$-point numerator and a $(n{-}\abs{I})$-point numerator! Unlike the usual factorization on the physical poles of the amplitude, these factorizations on the spurious poles stem from the combinatorial picture without any known physical origin. Explicitly
\begin{equation}\label{eq: Kfac}
\begin{aligned}
\left.\mathrm{Res}\right|_{D_I{=}0}K(1,2,\ldots,n){=}&D_I K( 1,\mathrm{Id}(I),P) \\
    &\times\tilde{K}^I(1,\mathrm{Id}(\bar{I}),n),
\end{aligned}
\end{equation}
where $P\equiv \bar{I}n$ denotes an effective scalar. For the definition of $\Delta(I_{k}, I_{k{+}1})$ in $K( 1,\mathrm{Id}(I),P)$, the complement of the set $I_k$ is still defined as $\bar{I_k}=\{2,3,\ldots,n{-}1\}/I_k$ while for $\Delta(J_{k}, J_{k{+}1})$ in $\tilde{K}^I(1,\mathrm{Id}(\bar{I}),n)$ the complement of $J_k$ is defined as $\bar{I}/J_k$.
The factor $D_I K(1,\mathrm{Id}(I),P)$ in \eqref{eq: Kfac} means that the overall pole $D_I$ of $K(1,\mathrm{Id}(I),P)$ is excluded. The factorization properties~\eqref{eq: Kfac} can be proved directly by plugging in the definitions on both sides.

For instance, at six points as shown in figure \ref{fig: P4}, there are $14$ co-dimension one boundaries $D_I=0$ including eight poles with $\abs{I}=1$ or $3$ corresponding to hexagons and six poles with $\abs{I}=2$ corresponding to squares. On any of the hexagon boundary, {\it i.e.} when $D_I=0$ with $\abs{I}=1$ or $3$, the residue factorizes into $\mathcal{F}_3=13$ terms (times $\mathcal{F}_1=1$ term). Similarly when $D_I=0$ with $\abs{I}=2$, the residue factorizes differently, {\it e.g.} as $D_{23}K(1,2,3,456)\times\tilde{K}^{23}(1,4,5,6)$ when $D_{23}=0$ (the square is the product of two line segments $\mathcal{P}_{\{23\}}\times\mathcal{P}_{\{45\}}$).

Algebraically, the quasi-shuffle algebra can be prompted to a bialgebra by introducing the coproduct map \cite{Hoffman}, and one can show the factorization properties from the coproduct. Similar to \cite{Brandhuber:2021bsf, Brandhuber:2022enp}, we can also define the antipode map to make the bialgebra a quasi-shuffle Hopf algebra. Acting on the BCJ numerators, the antipode map does nothing but changes its overall sign. The detail is given in the appendix \ref{sec:appAG}.

We expect the factorization properties of BCJ numerators to be the key for showing the cancellation of spurious poles in the amplitude. Such properties also suggest certain positive geometries (rather than just combinatorics) underlying these BCJ numerators, and we leave further investigations to future works.

\section{Heavy-mass effective field theory} \label{sec:heavy1}
In this section, we study YMS amplitudes and their BCJ numerators in the heavy-mass effective field theory (HEFT), which are obtained by taking the heavy-mass limit for a pair of massive scalars with momenta \cite{Brandhuber:2021kpo, Brandhuber:2021bsf,Brandhuber:2022enp}
\begin{equation}
    p_1^\mu=mv^\mu,\qquad p_n^\mu=-m v^\mu - k^\mu,
\end{equation}
where $v^2=1$ and we are interested in the limit $m \to \infty$; in other words, we will study the the expansion in $1/m$ of the BCJ numerators which we denote as $K_\mathrm{H}(1,2,\ldots,n)$, as well as that of $\phi^3$ amplitudes, which combine to give the resulting HEFT amplitude $A^\mathrm{H}(1,2,\ldots,n)$ at the leading order in $1/m$. Here $k^\mu$ is at the same order as gluon momenta, which stay finite at ${\cal O}(m^0)$ as $m \to \infty$.

\subsection{Heavy limit of YMS amplitudes}
We will make a  particular choice of the reference momenta: $q_I=v$ for all $I$, which dramatically simplifies formulae for BCJ numerators and give rise to poles similar to HEFT numerators  in~\cite{Brandhuber:2021bsf}. In fact, for $n=4$ such a choice reduces the BCJ numerator to one term, since $v\cdot F_a \cdot v$ vanishes for a single particle $a$
\begin{equation}\label{eq: hamp4}
    K_\mathrm{H}(1,2,3,4)= \frac{p_{1}\cdot F_{23}\cdot v} {p_{23} \cdot v} = -\frac{2m^2}{k^2}\ v \cdot F_{23}\cdot v,
\end{equation}
where in the second equality we have used $v\cdot k=-k^2/(2m)$ implied by the on-shell condition $p_n^2=m^2$.
Notice that $K_{\mathrm{H}}(1,3,2,4)=K_{\mathrm{H}}(1,2,3,4)$, thus the amplitude  $A^\mathrm{H}(1,2,3,4)$ becomes
\begin{align} \label{eq:Hamp4}
&(\frac{1}{s_{12}}+ \frac{1}{s_{23}}) K_\mathrm{H}(1,2,3,4)  -\frac{1}{s_{23}}K_\mathrm{H}(1,3,2,4)  \\\nonumber
=&-\frac{m}{k^2} \frac{v \cdot F_{23}\cdot v}{ v\cdot p_2}.
\end{align}
Physically, the final HEFT amplitude has the leading order $\mathcal{O}(m)$ \cite{Brandhuber:2021kpo}. In the above example, we can see that the numerators are at $\mathcal{O}(m^2)$, and the sum of the leading contribution of $\phi^3$ amplitudes at $\mathcal{O}(m^0)$, say $1/s_{23}$ times the corresponding numerators vanishes. Therefore, the sum of the contribution of $\phi^3$ amplitudes at the next order, {\it i.e.} $1/s_{12}$ from $A^{\phi^3}(1,2,3,4)$ times the numerator produces the HEFT amplitude as the first non-vanishing order. This is also the case for $n=5$.
However, for higher $n$, the numerator contains some additional terms with higher power of $m$.
To obtain the leading order contribution of HEFT final amplitude, we expand the numerators and the $\phi^3$ amplitudes in $m^{-1}$. Note the overall pole $D_{23\ldots n-1}=v\cdot k$ for BCJ numerators is proportional to $m^{-1}$, we first collect the numerator according to its {\it superficial order} of $m^{-1}$, {\it i.e.} terms with $(i-1)$ $p_1$'s in the numerator,
\begin{equation} \label{eq: Kexpand}
    K_{\mathrm{H}}(1,2,\ldots,n)= \sum_{i=2}^{\lfloor n/2 \rfloor } K_\mathrm{H}^{(i)}(1,2,\ldots,n),
\end{equation}
where the upper bound of the summation is $\lfloor n/2 \rfloor$ since $p_1 \cdot F_{a}\cdot v=0$
implies that the numerator should contain as many $p_1 \cdot F_{ab} \cdot v$ as possible to have the highest power of $p_1$. In the above expansion, $K_\mathrm{H}^{(i)}\equiv K_\mathrm{H}^{(i)}(1,2,\ldots,n)$ refers to terms with the superficial order $\mathcal{O}(m^i)$. For example, at six points we have the following terms for $K_\mathrm{H}^{(2)}$ and $K_\mathrm{H}^{(3)}$ respectively:
\begin{equation*}
       \frac{p_1\cdot F_{23} \cdot v\ p_{23}\cdot F_5 \cdot v\ p_{23}\cdot F_4 \cdot v } {v\cdot k\ v \cdot p_{45}\ v\cdot p_4 }, \frac{p_1 \cdot F_{25}\cdot v\ p_1\cdot F_{34}\cdot v} {v\cdot k\ v\cdot p_{34} }.
\end{equation*}
In fact, as explained in appendix \ref{sec:heavy}, the actual order of $K_\mathrm{H}^{(i)}$ is $\mathcal{O}(m^{2})$ for $i=2$ and $\mathcal{O}(m^{i-1})$ otherwise.

In the HEFT amplitude, we sum over all cubic graphs relevant at leading order, and for each graph with its propagator structure, its numerator is given by the corresponding commutator of $K_\mathrm{H}^{(i)}$~\cite{Bern:2011ia}. Nicely we observe that certain commutators of $K_\mathrm{H}^{(i)}$ actually vanish, and the end result is that only $K_\mathrm{H}^{(2)}$ contributes to the amplitude at the leading order! We have checked such vanishing results up to $n=10$, but we do not have an all-$n$ proof at the moment.

In fact, such vanishing results are better than what we need here, {\it i.e.} only $K^{(2)}_{\rm H}$ contributes to gauge-theory amplitudes at leading order. We have checked up to $n=10$ that the stronger vanishing results actually ensure that only $K^{(2)}_{\rm H}$ contributes to gravity amplitude, which is at order $\mathcal{O}(m^2)$, as obtained by double copy in HEFT. We leave more details in the appendix \ref{sec:heavy} with a proof of the simplest case. As a result of this conjecture, the amplitude is given by
\begin{equation} \label{eq: Hamp}
\begin{aligned}
    &A^\mathrm{H}(1,2 \ldots,n)=\sum_{\Theta^1 } \frac{K_\mathrm{H}^{(2)}(1,\Theta^1,n)}{d_{\Theta^1}}, \\
\end{aligned}
\end{equation}
where the summation is over nested commutators of depth $n{-}4$ (``co-depth" $1$) of the ordered set $(2,3,\ldots,n-1)$. For instance, at 5-point we sum over $\Theta^1=([2,3],4),(2,[3,4])$; $d_{\Theta^1}$ denotes the propagator denominator corresponding to the cubic tree associated with $\Theta^1$ (two sub-trees on the scalar line $(1 n)$):
\begin{equation*}
\begin{aligned}\includegraphics[width=0.23\linewidth]{figs/general1.pdf}\end{aligned} \leftrightarrow d_{\Theta^{1}}, {\it e.g.} \begin{aligned}\includegraphics[width=0.23\linewidth]{figs/theta511.pdf}\end{aligned} \leftrightarrow  d_{[2,3],4}=s_{123}s_{23}.
\end{equation*}
Moreover, it is easy to show (see appendix \ref{sec:heavy} for details) that the effective BCJ numerator $K_\mathrm{H}^{(2)}(1,2,\ldots,n)$ contains $\mathcal{F}_{n{-}3}$ terms, and its pole structure corresponds to the permutohedron $\mathcal{P}_{\{34\ldots n{-}1\}}$, which means that
\begin{equation}\label{eq: BDexph}
    K_\mathrm{H}^{(2)}(1,2,\ldots,n)=\sum_{d=0}^{n-2}\sum_{\Gamma_d\in\partial^d\mathcal{P}_{n{-}3}}K^{ (2)}_{H,\Gamma_d}(1,2,\ldots,n).
\end{equation}
For the boundary $\Gamma_d{=}\{I_0,I_1,\ldots,I_d\}{\in}\partial^d\mathcal{P}_{\{34\ldots n{-}1\}}$ where $I_d{\subset}I_{d-1}{\subset}\ldots{\subset}I_0{=}\{3,4,\ldots,n{-}1\}$ and $I_d{\neq}\emptyset$, the contribution is
\begin{align}
    K^{(2)}_{\mathrm{H},\Gamma_d}&= \frac{m v\cdot F_{\tau_{0}}\cdot v}{p_{23\ldots n{-}1}\cdot v} \prod_{k=1}^{d}\frac{p_{\Delta(I_{k},I_{k{+}1})}\cdot F_{\tau_d}\cdot v}{v\cdot p_{I_k}}\\\nonumber
    &=-\frac{2m^2 v\cdot F_{\tau_{0}}\cdot v}{k^2} \prod_{k=1}^{d}\frac{p_{\Delta(I_{k},I_{k{+}1})}\cdot F_{\tau_d}\cdot v}{v\cdot p_{I_k}},
\end{align}
where in the calculation of $\Delta(I_{k},I_{k{+}1})$, the complement set of $I_k$ is still taken to be $\{2,3,\ldots,n{-}1\}/I_k$. For $n=4$, there is no commutator in $\Theta^1$ and the result is \eqref{eq:Hamp4}.
For $n=5$, the amplitude becomes
\begin{equation}
\begin{aligned}
        A^\mathrm{H}(1,2,3,4,5){=} &\frac{1}{s_{12}} \frac{K_\mathrm{H}^{(2)}(1,2,[3,4],5)}{s_{34}}  \\
        & {+}\frac{1}{s_{123}} \frac{K_\mathrm{H}^{(2)}(1,[2,3],4,5)}{s_{23}},
\end{aligned}
\end{equation}
where $K_\mathrm{H}^{(2)}(1,2,3,4,5)$ is given by{\small
\begin{equation*}
 -\frac{2m^2}{k^2}(v \cdot F_{234} \cdot v+\frac{ v\cdot F_{24} \cdot v\ p_{2}\cdot F_3\cdot v }{v\cdot p_3 }
     + \frac{ v\cdot F_{23} \cdot v \ p_{23}\cdot F_4\cdot v }{v\cdot p_4}).
\end{equation*}}
Let us give a final example for $n=6$ amplitude
\begin{equation}
\begin{aligned}
    & \frac{1}{s_{12}} \left(\frac{K_\mathrm{H}^{(2)}(1,2,[[3,4],5],6)}{s_{34}s_{345}}{+} \frac{K_\mathrm{H}^{(2)}(1,2,[3,[4,5]],6)}{s_{45}s_{345}} \right) \\
    {+}& \frac{1}{s_{123}} \frac{K_\mathrm{H}^{(2)}(1,[2,3],[4,5],6)}{s_{23}s_{45}}  \\
    {+}&\frac{1}{s_{1234}} \left(\frac{K_\mathrm{H}^{(2)}(1,[[2,3],4],5,6)}{s_{23}s_{234}}{+} \frac{K_\mathrm{H}^{(2)}(1,[2,[3,4]],5,6)}{s_{34}s_{234}} \right).
\end{aligned}
\end{equation}
It is interesting to notice the numerators we present here only differ from those in \cite{Brandhuber:2021bsf} denoted by $N(1,2,\ldots,n)$ by an overall prefacto
\begin{equation}\label{eq: HBCJnumerator}
    K_\mathrm{H}^{(2)}(1,2,\ldots,n)=(-1)^{n}(n{-}2)\frac{2m}{k^2} v\cdot p_2 N(1,2,\ldots,n).
\end{equation}
It is highly nontrivial, however, that these two sets of effective BCJ numerators give the same HEFT amplitude. In~\cite{Brandhuber:2021bsf}, the expression involves the sum of cubic graphs corresponding to nested commutators of depth $n-3$ of the ordered set $(2,3,\ldots,n{-}1)$, thus the propagator denominator contains an overall factor $s_{23\ldots n-1}$, which in our case is replaced by different $s_{1 \sigma}$ for different terms. In addition, the numerator of~\cite{Brandhuber:2021bsf} for each cubic graph is given by a nested commutator of $N(1,2,\ldots,n)$, thus the number of terms in it is twice as ours. Nevertheless, we have analytically checked up to $n=10$ that the amplitude \eqref{eq: Hamp} agrees with \cite{Brandhuber:2021bsf}. Moreover, we have checked that although they look very different, the HEFT gravity amplitude via double copy also agrees with that in~\cite{Brandhuber:2021bsf}, and we expect both agreements to hold for all $n$.

\subsection{Decoupling into pure YM}
Given the explicit result of the $n$-point heavy mass BCJ numerators, the $(n-1)$-point pure YM BCJ numerators, as well as the amplitudes, can be easily obtained via the decoupling limit:
$mv\to \varepsilon_n$, $p_{23\ldots n-1}^2\to 0$
to obtain the BCJ numerator $K^{\prime \text{YM}}(2,3,\ldots,n)$ \cite{Brandhuber:2021kpo,Brandhuber:2021bsf}. Under this kinematics, the overall factor $k^2$ vanishes, which we ignore in the decoupling limit. For instance, the three-point BCJ numerator is given by $K^{\prime \text{YM}}(2,3,4)=-2\varepsilon_4 \cdot F_{23}\cdot \varepsilon_4$. Therefore, the three-point amplitude is
\begin{align}
    A^{\text{YM}}(2,3,4)&=  -\frac{\varepsilon_4 \cdot F_{23}\cdot \varepsilon_4}{ \varepsilon_4 \cdot p_2}\\\nonumber
    &=\varepsilon _4\cdot \varepsilon _2 p_2\cdot \varepsilon _3{-}\varepsilon _2\cdot \varepsilon _3 p_2\cdot \varepsilon _4{-}\varepsilon _4\cdot \varepsilon _3 p_3\cdot \varepsilon _2.
\end{align}
For the 4-point YM amplitude, the numerator $K^{\prime \text{YM}}(2,3,4,5)$ reads
\begin{equation*}
    {-}2\left(\varepsilon_{5}{\cdot}F_{234}{\cdot} \varepsilon_5{+}\frac{\varepsilon_5{\cdot} F_{24}{\cdot} \varepsilon_5 p_{2}{\cdot} F_{3}{\cdot} \varepsilon_5}{p_3{\cdot}\varepsilon_5 }{+}\frac{\varepsilon_5{\cdot} F_{23}{\cdot} \varepsilon_5 p_{23}{\cdot} F_4{\cdot} \varepsilon_5 }{ p_4{\cdot}\varepsilon_5}\right).
\end{equation*}
Note that $K^{\prime \text{YM}}(2,3,\ldots,n)$ also manifests the gauge invariance of particles $2,3,\ldots,n-1$. Moreover, it is related to the BCJ numerator given in sec.\ref{sec:YM} via
\begin{equation*}
    K^{\prime \text{YM}}(2,3,\ldots,n)=2\left. \varepsilon_n \cdot p_2 K^{\text{YM}}(2,3,\ldots,n)\right|_{q_I\to \varepsilon_n}.
\end{equation*}
These numerators, accompanied by different $\phi^3$ amplitudes, produce the same YM amplitude.

\section{Conclusions and outlook}

In this note, we established a correspondence between BCJ numerators from covariant color-kinematics duality and the combinatorial permutohedra, which are closely related to the quasi-shuffle Hopf algebra. This apply to all YMS amplitudes, but the most interesting case is that with two scalars, whose numerators share the same combinatorial structure as the pure YM ones: each term is mapped to a boundary of $\mathcal{P}_{n-2}$; the contribution from each boundary is almost identical in these two cases, except that we need to modify one factor to take into account the remaining two gluons. We also found nice recursion relations and factorization properties implied by this picture. Finally, based on highly nontrivial cancellations which are needed for both YMS and gravity amplitudes (via double copy) in HEFT, we conjectured a compact formula for their effective numerators; they become closely related to permutohedra $\mathcal{P}_{n-3}$, which, while producing the same amplitude, differ by an overall factor from the numerators in ~\cite{Brandhuber:2021bsf, Brandhuber:2022enp}.

There are numerous open questions for further investigations. First, as we will present in \cite{toapp}, it is interesting to see how lower-dimensional permutohedra for general YMS numerators combine into $\mathcal{P}_{n-2}$ which corresponds to the pure YM ones; we also find interesting combinatorial structures underlying BCJ numerators of amplitudes in NLSM {\it etc.}. Moreover, the somewhat miraculous cancellations that simplify these numerators in HEFT still remain to be proven, which would also be important to establish the correct double copy in HEFT. Since the final amplitudes are independent of reference momenta, all the spurious poles must cancel, which still calls for a direct understanding (without relying on the CCK duality); such an understanding could connect this combinatorial picture (especially the factorizations) to the uniqueness theorem for YM amplitude~\cite{Arkani-Hamed:2016rak, Rodina:2016jyz} and YMS ones via the universal expansion~\cite{Dong:2021qai}. Last but not least, it is tempting to ask: could we combine the permutohedra for BCJ numerators with the associahedra for bi-adjoint $\phi^3$ amplitudes, and obtain a unified geometric understanding of gluon and graviton scattering?

\begin{acknowledgments}
We thank Linghui Hou, Guanda Lin, Tianheng Wang for discussions and collaborations on related projects. This research is supported in part by the Key Research Program of CAS, Grant No. XDPB15 and National
Natural Science Foundation of China under Grant No.
11935013,11947301,12047502,12047503
\end{acknowledgments}

\newpage

\widetext