Then we find that the spectrum of the observed final state, for instance a , is obtained by integrating over x$_{b}$ . Kinematics constrains x$_{a}$ to be in the region where the approximation is valid for small mass limit ϵ$_{a,b,c}$ ≪ 1, and the triangle function is defined as λ ( x, y, z ) = x 2 + y 2 + z 2 − 2 xy − 2 xz − 2 yz . The general process is S → ψ$_{1}$ψ$_{2}$φ$_{3}$ with amplitude given in Eq. (A1). In terms of x and ϵ , we can write the amplitude as follows |M| 2 =2( CC $_{I}$) $^{2}$Y 2 ( 1 + ϵ 2 3 − ( ϵ$_{2}$ + ϵ$_{1}$ ) 2 − x$_{3}$ ) , (B4) which is symmetric for the subscripts 1 and 2. The distribution function of the fermion ψ$_{1}$ is obtained by paring ψ$_{1}$ and φ$_{3}$ , i.e. , setting a = ψ$_{1}$, b = φ$_{3}$, c = ψ$_{2}$, where N = 1 / Γ is the normalization factor. In the massless limit, it is approximated to be a simple function x 2 $_{1}$. So the fermion spectrum is very hard. The distribution of ψ$_{2}$ is got simply by replacing the subscript 1 by 2 in the above formula. Similarly, we get the distribution function of scalar φ$_{3}$ by choosing a = φ$_{3}$, b = ψ$_{1}$ , and c = ψ$_{2}$ dN dx$_{3}$ = 2( CC $_{I}$) $^{2}$Y 4 256 π $^{3}$m$_{D}$ N ∫ R + 1 ( x 3 ) R − 1 ( x 3 ) dx$_{1}$ ( 1 + ϵ 2 3 − ( ϵ$_{2}$ + ϵ$_{1}$ ) 2 − x$_{3}$ ) = 2( CC $_{I}$) $^{2}$Y 4 256 π $^{3}$m$_{D}$ N ( 1 + ϵ 2 3 − ( ϵ$_{2}$ + ϵ$_{1}$ ) 2 − x$_{3}$ ) ( R + 1 ( x$_{3}$ ) − R − 1 ( x$_{3}$ ) ) . (B6) dN dx$_{1}$ = 2( CC $_{I}$) $^{2}$Y 4 256 π $^{3}$m$_{D}$ × N ∫ R + 3 ( x 1 ) R − 3 ( x 1 ) dx$_{3}$ ( 1 + ϵ 2 3 − ( ϵ$_{2}$ + ϵ$_{1}$ ) 2 − x$_{3}$ ) = 2( CC $_{I}$) $^{2}$Y 4 256 π $^{3}$m$_{D}$ N [ 1 + ϵ 2 3 − ( ϵ$_{2}$ + ϵ$_{1}$ ) 2 − 1 2 ( R + 3 ( x$_{1}$ ) + R − 3 ( x$_{1}$ ) ) ] × ( R + 3 ( x$_{1}$ ) − R − 3 ( x$_{1}$ ) ) , (B5) R ± b = 1 2 (1 − x$_{a}$ + ϵ 2 $_{a}$) − $^{1}$[(2 − x$_{a}$ )(1 + ϵ 2 a + ϵ 2 b − ϵ 2 c − x$_{a}$ ) ± √ x 2 a − 4 ϵ 2 $_{a}$λ 1 2 ( 1 + ϵ 2 a − x$_{a}$, ϵ 2 $_{b}$, ϵ 2 c ) ] ≈ 1 2 ((2 − x$_{a}$ ) ± x$_{a}$ ) , (B3) 2 ϵ$_{a}$ ⩽ x$_{a}$ ⩽ 1 + ( ϵ 2 a − ϵ 2 b − ϵ 2 c − 2 ϵ$_{b}$ϵ$_{c}$ ) = R max a , (B2) d$_{3}$ ( P S, Y → abc ) = Y 2 128 π $_{3}$dx$_{a}$ dx$_{b}$ . (B1) space is reduced to the two-body phase space by pairing a with b and further integrating over c where ϵ 2 a ≡ m 2 $_{a}$/Y $^{2}$. If all the final state masses are ignorable, the region is between 0 and 1. For fixed x$_{a}$ , x$_{b}$ is constrained within the region [ R − b , R + b ] 19