Oscillating Composite Asymmetric Dark Matter
Masahiro Ibe Shin Kobayashi Ryo Nagai Wakutaka Nakano
hep-ph astro-ph.HE




α
β̱
ε̧
δ̣

ϕ
γ
θ
κ̨
łλ
μ
ν
ψ
∂
ρ̊
σ
τ
ῠ
φ̌
ω
ξ
η
ζ
Δ
Γ
Θ̋
Ł
Φ
Ψ
Σ

ø

    diag Spin SOØ O SU U Sp SL trM_PLInt. J. Mod. Phys. Mod. Phys. Lett. Nucl. Phys. Phys. Lett. Phys. Rev. Phys. Rev. Lett. Prog. Theor. Phys. Z. Phys. df BHinf inf evap eq smM_ PL GeV
    [e-mail: ]ibe@icrr.u-tokyo.ac.jpKavli IPMU (WPI), UTIAS, University of Tokyo, Kashiwa, Chiba 277-8583, JapanICRR, University of Tokyo, Kashiwa, Chiba 277-8582, Japan[e-mail: ]shinkoba@icrr.u-tokyo.ac.jpICRR, University of Tokyo, Kashiwa, Chiba 277-8582, Japan[e-mail: ]rnagai@icrr.u-tokyo.ac.jpICRR, University of Tokyo, Kashiwa, Chiba 277-8582, Japan[e-mail: ]m156077@icrr.u-tokyo.ac.jpICRR, University of Tokyo, Kashiwa, Chiba 277-8582, Japan
    The asymmetric dark matter (ADM) scenario can solve the coincidence problem between
    the baryon and the dark matter (DM) abundance when the DM mass is of O(1) GeV.
    In the ADM scenarios, composite dark matter is particularly motivated,
    as it can naturally provide the DM mass in the O(1) GeV  range 
    and a large annihilation cross section simultaneously.
    In this paper, we discuss the indirect detection constraints on the composite ADM model.
    The portal operators connecting the B-L asymmetries in the dark and the Standard Model(SM) sectors
    are assumed to be generated in association with the seesaw mechanism.
    In this model, composite dark matter inevitably obtains a tiny Majorana mass 
    which induces a pair-annihilation of ADM at late times.
    We show that the model can be efficiently tested by the searches for the γ-ray 
    from the dwarf spheroidal galaxies and the interstellar electron/positron flux.
    Oscillating Composite Asymmetric Dark Matter
        Wakutaka Nakano
        December 30, 2023
    ============================================
    
    IPMU 19-0103
    
    
    
    
    § INTRODUCTION
    
    Asymmetric dark matter (ADM) scenario sheds light on 
    the coincidence problem between the observed baryon and dark matter (DM) abundances in the universe
     <cit.> 
    (see also <cit.> for reviews).
    If the DM abundance is provided by a mechanism which is unrelated to the baryogenesis, it is quite puzzling why those 
    abundances are close with each other despite the fact that the baryon abundance is dominated by 
    the contribution from the matter-antimatter asymmetry. 
    In the ADM scenario the coincidence problem can be explained when the DM mass is of O(1) GeV,
    where the matter-antimatter asymmetry is thermally distributed between the dark and the Standard Model (visible) sectors.
    
    Among various ADM scenarios, composite baryonic DM in QCD-like dynamics is particularly motivated
    since it can naturally provide a large annihilation cross section and  the DM mass in the GeV range 
    simultaneously <cit.>.
    Recently, a minimal composite ADM model and its ultraviolet (UV) completion <cit.> have been proposed where the asymmetry 
    generated by the thermal leptogenesis <cit.> 
    (see also <cit.> for review) 
    is thermally distributed between the two sectors through 
    a portal operator associated with the seesaw 
    mechanism <cit.>.
    The dark sector of the model consists of QCD-like dynamics and QED-like interaction, 
    which are called as dark QCD and dark QED, respectively.
    The lightest baryons of dark QCD play the role of ADM.
    The dark QED photon (dark photon) obtains a mass of 10100 MeV,
    which plays a crucial role to transfer the excessive entropy of the dark sector
    into the visible sector before neutrino decoupling <cit.>.
    
    In this paper, we discuss the indirect detection of the composite ADM model in <cit.>.
    The portal operator in this model is generated in association with the seesaw mechanism.
    In this model, the dark-neutron, one of the lightest dark baryons, inevitably obtains a tiny Majorana mass. 
    Such a tiny Majorana mass induces the oscillation between DM particle and the antiparticle,
    which induces a pair-annihilation of ADM at late times <cit.>.
    A pair of DM particle and the antiparticle annihilates into multiple dark pions, 
    and the (neutral) dark pion subsequently decays into a pair of the dark photons.
    The dark photon eventually decays into an electron-positron pair. 
    Thus, the late time annihilation of ADM results in multiple 
    soft electrons/positrons.
    In addition, soft photons are also emitted as final state radiation.
    As we will see, the model can be efficiently tested 
    by the searches for the γ-ray from the dwarf spheroidal galaxies (dSphs) by the Fermi-LAT.
    We also discuss the constraints from 
    the observations of the interstellar electron/positron flux by the Voyager-1.
    
    The organization of the paper is as follows.
    In section <ref>, we review the composite ADM model in <cit.> and show 
    how the tiny Majorana mass of the dark neutron appears associated with the seesaw mechanism.
    In section <ref>, we derive the expected γ-ray flux from the dSphs and discuss the constraints 
    on the model by comparing the flux with the Fermi-LAT results. 
    We also estimate the interstellar electron/positron flux in cosmic ray from the late time annihilation and compare it with the Voyger-1 result.
    The final section is devoted to the conclusions. 
    
    
    
    
    § DM ANTI-DM OSCILLATION IN THE COMPOSITE ADM MODEL
    
    
     §.§ A Model of Composite ADM
    
    In this subsection, we briefly review the composite ADM model in <cit.>.
    The model is based on N_g-generation dark quarks with SU(3)_D × U(1)_D gauge symmetry.
    SU(3)_D provides the dark QCD dynamics and U(1)_D the dark QED interaction.
    The dark quarks are the fundamental representations of SU(3)_D.
    They are charged under the dark QED and the B-L in analogy to
    the up-type and the down-type quarks in the visible sector (see Tab. <ref>).
    They have tiny masses, 
    
        L_mass = m_U' U' U' + m_D'D' D' +  ,
    
    with m_U' and m_D' being the mass parameters.
    Hereafter, we put primes on the parameters and the fields in the dark sector
    when there  are counterparts in the visible sector.
    
    The dark QCD exhibits confinement below the dynamical scale of SU(3)_D, 
    Λ_ QCD', which leads to the emergence of the dark baryons and the dark mesons.
    Throughout this paper, we assume that only one generation of the dark quarks have masses 
    smaller than Λ_ QCD'.
    [For N_g >1, we assume the heavier dark quarks decay into the lighter 
    ones by emitting the dark Higgs boson which has the dark QED charge of 1.
    It should be noted that the dark quark masses are not generated by the vacuum expectation 
    value of the dark Higgs boson <cit.>,
    and hence, the dark Higgs couplings generically violate the flavor symmetry in the dark sector.]
    The lightest dark baryons, i.e. the dark nucleons, 
    
        p' ∝ U'U'D' ,   p̅' ∝U̅'U̅'D̅' ,   
        n' ∝ U'D'D' ,   n̅' ∝U̅'D̅'D̅' ,
    
    are stable in the decoupling limit from the visible sector due to their B-L charges.
    Once the B-L asymmetry is shared between the visible and the dark sector, 
    the dark nucleon abundance is dominated by the asymmetric component due to their large annihilation cross section.
    Therefore, the dark nucleon with a mass in the GeV range is a good candidate for ADM.
    
    
    
    
    
    
    When the B-L asymmetry is thermally distributed between the visible and the dark sectors,
    the ratio of the B-L asymmetry stored in each sector is given by A_ DM/A_ SM = 44 N_g/237
    for the B-L charges given in Tab. <ref> <cit.>.
    [In the presence of additional B-L charged fields in the dark sector, such as dark leptons,
    the ratio can be modified.
    Besides, the neutrality condition of U(1)_D and the contributions from the dark Higgs sector
    also change the ratio by some tens percent for a given N_g.]
    Thus, the observed ratio of the DM and the baryon abundance can be reproduced when the dark nucleon mass is
    
        m_N' ≃Ω_ DM/Ω_ BA_ B/A_ SMA_ SM/A_ DM× m_N
        ≃8.5  GeV/N_g .
    
    Here, we have used the ratio of the baryon asymmetry to the B-L asymmetry in the visible sector, A_ B/A_ SM = 30/97 <cit.>.
    The dark nucleon mass in this range can be naturally realized when Λ_ QCD' is in the GeV range.
    
    
    
    
    
    The lightest dark mesons,
    
        π'^0∝ U' U̅' - D' D̅'  ,   π'^+ ∝ U'D̅'
         ,   π'^- ∝ D'U̅' ,
    
    annihilate or decay into the dark photons.
    As a result, they do not contribute to the effective number of neutrino degrees of freedom
    nor to the dark matter abundance significantly even if they are stable.
    In the following analysis, we assume that the dark charged pions are stable for simplicity.
    [If U(1)_D is broken by the vacuum expectation value of  a dark Higgs with the dark QED charge of 2, 
    a ℤ_2 symmetry remains unbroken which makes the dark charged pion stable.
    If U(1)_D is broken by the dark Higgs with the charge 1, the neutral and the charged pions can mix each other, and hence, the charged pions decay.]
    The decay of the dark neutral pion into a pair of dark photons, on the other hand, is inevitable due to the chiral anomaly.
    As we will see, the decay of the neutral pion plays a central role for the indirect detection of ADM.
    
    The dark photon obtains its mass by the dark Higgs mechanism, 
    and it decays into the visible fermions thought the kinetic mixing with the visible QED photon, 
    
        L_γ^' = ϵ/2 F_μν F^'μν + m_γ^'^2/2 A_μ^' A^'μ .
    
    Here, F_μν and F^'_μν denote the field strengths of 
    the visible and the dark QED with A^'_μ being the dark photon gauge field. 
    In the following, we assume the kinetic mixing parameters of ϵ = 10^-10–10^-8
    and the dark photon mass in 10100 MeV range which 
    satisfies all the constraints <cit.> (see also <cit.>).
    [See <cit.> for discussion on the origin of the tiny kinetic mixing parameters.]
    In this parameter range, the dark photon decays when the cosmic temperature is above O(1) MeV.
    
    Finally, let us comment on the ratio between the abundances of the dark protons and the dark neutrons.
    In the present model, there is no dark leptons nor dark weak gauge bosons.
    Besides, it is expected that 
    the mass difference 
    between the dark neutron and the dark proton is 
    smaller than the mass of the dark pion
    when the dark quark masses are smaller than the dynamical
    scale of SU(3)_D.
    Thus, the dark neutron is stable in the limit of the 
    vanishing B-L portal interactions (see below).
    The ratio between the dark proton abundance 
    and the dark neutron abundance is given by <cit.>,
    
        n_n'/n_p' ∼ e^-(m_n'-m_p')/T_F .
    
    Here, n_n,p' and m_n,p' are the number densities 
    and the masses of the dark neutron and the dark proton, respectively.
    T_F denotes the freeze-out temperature of the dark pion annihilation,
    T_F ≃ m_π'/ O(10).
    Thus, for m_n'-m_p' ≪ m_π', the dark neutron abundance is comparable 
    to that of the dark proton.
    In the following, we take n_n' = n_p'.
    
    
    
    
     §.§ The  portal operator
    
    The B-L asymmetry generated by thermal leptogenesis is thermally distributed between the visible
    and the dark sectors.
    For this purpose, there need to be portal interactions which connect the B-L symmetry in the two sectors.
    In the model in <cit.> (see also <cit.>), the following operators are assumed 
    as the portal operators,
    
        L_portal∼1/M_*^3 (U' D' D') (L H)
        + 1/M_*^3 (U'^† D'^†D') (L H) +  ,
    
    where L and H are the lepton and the Higgs doublets in the visible sector, and M_* is a dimensional parameter.
    [The portal operators require the gauge invariant operators which are charged under the B-L symmetry. 
    This is the reason why we need both the up-type and the down-type quarks in the dark sector.]
    Here, we omit the 1 coefficients.
    The effects of the above operators decouple at the cosmic temperature below 
    T_*∼ M_* (M_* / )^1/5.
    Here, = 2.4 × 10^18 GeV denotes the reduced Planck mass.
    For successful ADM with thermal leptogenesis, the decoupling 
    temperature, T_*, is required to be lower than the temperature, T_B-L, at which leptogenesis completes.
    In the following, we consider the so-called strong washout regime of thermal leptogenesis, 
    where the leptogenesis completes at the temperature about T_B-L≃ M_R/z_B-L with
    z_B-L≃ 10 <cit.>.
    
    In <cit.>, the UV model has been proposed 
    in which the portal operators in Eq. (<ref>) 
    are generated by integrating out the right-handed neutrinos, N̅, and the dark colored Higgs boson, H_C'.
    The gauge charges of H_C' are identical to those of D', while H_C'
    has the B-L charge -2/3.
    The right-handed neutrinos couple to both sectors via,
    
        L = M_R/2N̅N̅ + y_N L H N̅ +1/2M_C^2 |H_C'|^2-Y_N H_C'D'N̅ 
        -Y_C H_C' U' D' - Y_C̅ H_C'^†U' D'
        +
    
    Here, M_C denotes the dark colored Higgs mass, M_R the mass of the right-handed neutrinos, 
    and y_N and Y's are the Yukawa coupling constants.
    The flavor and the gauge indices are suppressed.
    It should be noted that the mass terms of the right-handed neutrino 
    break B-L symmetry explicitly.
    The first two terms are relevant for the seesaw mechanism.
    
    
    By integrating out N̅ and H_C' from Eq. (<ref>), the portal operators in Eq. (<ref>) are obtained 
    where M_* corresponds to 
    
        1/M_*^3 = y_NY_NY_C̅/2M_C^2 M_R ,   1/M_*^3 = y_NY_NY_C^*/2M_C^2 M_R ,
    
    for each term of Eq. (<ref>), respectively.
    From the condition of T_* < T_B-L, the mass of the dark colored Higgs should satisfy,
    [Hereafter, we take M_R >0 and neglect the complex phases of Y_N, Y_C and Y_C̅.
    We also assume Y_C = Y_C̅ for simplicity.]
        M_R/z_B-L≲ M_C≲10/z_B-L^5/4(m̂_ν/0.1 eV)^1/4√(Y_N Y_C) M_R .
    
    The first inequality comes from a consistency condition of the decoupling limit of 
    the dark colored Higgs
    at the temperature T_B-L.
    In the right hand side, we have reparameterized the neutrino Yukawa coupling by using a tiny neutrino mass parameter, 
    m̂_ν,
    
        |y_N^2| ∼ 10^-5( m̂_ν/0.1  eV) ( M_R/10^9 GeV)  .
    
    
    
    Incidentally, 
    the dark nucleon can decay into the dark pion and the anti-neutrino in the visible sector through the B-L portal operator in Eq. (<ref>) <cit.>.
    The lifetime is roughly given by,
    
        τ_N' ≃ 10^33  sec(2  GeV/Λ_ QCD')^4
              (0.1  eV/m̂_ν)
              (M_R/10^9  GeV)
            (M̃_C/3× 10^9  GeV)^4
                (10  GeV/m_ DM) ,
    
    where M̃_C = M_C/√(Y_NY_C).
    Thus, the lifetime of the dark nucleons is 
    much longer than the age of the universe 
    for M_R ∼ M_C ∼ 10^9 GeV.
    
    
    
    
     §.§ The Majorana mass of the dark neutron
    
    The portal operators in Eq. (<ref>) are generated in association with the seesaw mechanism.
    As a notable feature of the UV completion model in Eq. (<ref>), it 
    also leads to the Majorana mass term of the dark neutron.
    This can be observed by integrating out H_C' and N̅ one by one.
    In the case of M_C>M_R, we first integrate out H_C' from Eq. (<ref>), which reads
    	
        L =    N̅^† iσ^μ∂_μN̅ + 
        		M_R/2N̅N̅ + y_N L H N̅
           -Y_N^2/2M_C^2D̅'N̅(D̅'N̅)^† -Y_N/2M_C^2D̅'N̅[ (Y_CU'D')^† + Y_C̅U̅'D̅']
        		+ 
           +(quartic in dark quark fields).
    
    Here, we show the kinetic term of N̅ explicitly which were implicit in Eq. (<ref>).
    This formula is of the form
    	
        L=(AN̅N̅ + BN̅+)-CN̅N̅^†,
    
    where
    [Here, χ^†σ^μ↔∂_μη 
    =χ^†σ^μ∂_μη
    -∂_μχ^†σ^μη
    for the Weyl fermions, χ  and η.
    ]
        A=M_R/2,   B=y_NLH-Y_N/2M_C^2D̅'[ Y_C̅U̅'D̅' + (Y_CU'D')^†],   C=-i/2σ^μ↔∂_μ+Y_N^2/2M_C^2D̅'D̅'^†.
    
    To make N̅ integrated out, it is convenient to complete the square of Eq. (<ref>) with respect to N̅.
    For this purpose, we shift N̅ by N̅→N̅+ψ̅, with which we can eliminate the linear term 
    in Eq. (<ref>).
    The condition ψ̅ must satisfy is 2Aψ̅+B-Cψ̅^†=0, which reads
    	
        ψ̅=-2A^*B+CB^†/4A^2-C^2≃ -1/2M_R^2(1+C^2/M_R^2)(M_R^*B+CB^† )
    
    After the shift, we integrate out N̅ to obtain
    	
        L   =(Aψ̅ψ̅ + Bψ̅+)-Cψ̅ψ̅^†
        
        					   =1/2Bψ̅+
        
        					   =-(1+C^2/M_R^2)(1/2M_RBB+1/2M_R^2CB^† B)+
    
    From Eq. (<ref>), we find that BB term includes the Mojorana mass term of the dark neutron
    	
        1/2M_R BB ⊃Y_N^2Y_C̅^2/8M_R M_C^4(U̅'D̅'D̅')^2∼Y_N^2Y_C̅^2Λ_ QCD'^6/8M_R M_C^4n̅'n̅'.
    
    In this way, Eq. (<ref>) leads to the Majorana mass,
    
        m_M = Y_N^2Y_C̅^2Λ_ QCD'^6/4M_R M_C^4= Λ_ QCD'^6/4M_R M̃_C^4 ,
    
    in addition to the B-L portal  operators in Eq. (<ref>).
    
    Once the dark neutron obtains the Majorana mass, the dark neutron and the anti-dark neutron
    oscillate with a time scale of t_ osc = m_M^-1 <cit.>.
    The probability to find an anti-dark neutron at a time t
    is given by,
    
        P(n'↔n̅') = sin^2(m_Mt)  .
    
    Here, we assume that the initial state at t = 0 is a pure dark neutron state.
    As we will see in the next section, the oscillation 
    induces a pair-annihilation of ADM which ends up with multiple 
    soft electrons/positrons/photons.
    
    
    
     §.§ Washout Interactions and On-Shell Portal
    
    Before closing this section, let us discuss the B-L washout interactions which 
    are also induced from Eq. (<ref>).
    In fact, the term CB^† B  in 
    Eq. (<ref>)
    includes
    	
        1/2M_R^2 CB^† B⊃y_NY_N^*/2M_C^2 M_R^2[Y_C̅^*(U'^†D'^†D'^†)
        																			(iσ^μ∂_μ)(L H)
        																				
        																				+Y_C(U' D' D'^†)(iσ^μ∂_μ)(L H)] .
    
    	In these interaction terms, and those in Eq. (<ref>), L couples to the dark sector operators which have the opposite B-L charges with each other.
    	Thus, if these operators are also in equilibrium
    	at T_B-L, the B-L asymmetry generated by leptogenesis is washed out.
    	To avoid such problems, it is required that
    	
        10/z_B-L^7/4(m̂_ν/0.1 eV)^1/4 M_R	≲M̃_C .
    
    	By comparing Eqs. (<ref>) and (<ref>), we find that the allowed parameter region for the ADM scenario is highly restricted due to the washout interaction when the portal operators are generated from the UV model in Eq. (<ref>).
    	
    	This constraint can be easily relaxed by introducing additional B-L portals.
    	For example, we may introduce a pair of gauge singlet fermions, (X,X̅) 
    	with new scalar fields, H_p, and H_Cp', whose gauge and B-L charges are 
        the same with those of the Higgs doublet of the SM and the dark colored Higgs, respectively.
    	In this case, there can be additional operators,
    	
        L= 
        	    M_X X X̅ +1/2M_H^2 |H_p|^2+1/2M_Cp^2 |H_Cp'|^2 
        	    +y_X L H_p X̅ -Y_X H_Cp'D'X̅ .
    
    	Here, M_X, M_H and M_Cp are the mass parameters of (X,X̅), H_p and H_Cp', respectively, and y_X and Y_X are Yukawa coupling constants.
    	[N̅ and X̅ can be distinguished by an approximate discrete symmetry under which (X,X̅), H_p and H_Cp' are charged. With the discrete symmetry, we can avoid unnecessarily mixing between N̅ and X̅. 
    	]
    	As the mass of X̅ is the Dirac type, the interaction terms in Eq. (<ref>) do not violate the B-L symmetry.
    	Thus, these interactions do not washout the asymmetry generated by leptogenesis 
    	but thermally distribute the asymmetry between the visible and the dark sector for M_X,H,Cp < T_B-L.
    	
    	
    	
    	In the following analysis, we divide the parameter region into two. 
    	
      *  Off-shell B-L portal scenario:
    	    
        10/z_B-L^7/4(m̂_ν/0.1 eV)^1/4 M_R	≲M̃_C≲10/z_B-L^5/4(m̂_ν/0.1 eV)^1/4M_R .
    
      *  On-Shell B-L portal scenario:
    	 	    
        10/z_B-L^5/4(m̂_ν/0.1 eV)^1/4 M_R≲M̃_C  .
    
    In the on-shell portal scenario, we assume that there are lighter particles than T_B-L which mediate the B-L asymmetry between two sectors as in Eq. (<ref>).
    [In the on-shell scenario, we may take Y_N = 0, and hence, the Majorana dark neutron mass is not inevitable.]
    It should be emphasized that the B-L asymmetries in the two sectors are thermally distributed in both the scenarios.
    [In the absence of the on-shell portal, 
    the region with M_C < T_B-L results 
    in a dark sector asymmetry which depends on the branching ratio
    of N̅ for small Y_N's <cit.>.
    If Y_N's are large for M_C < T_B-L, on the other hand, the B-L asymmetry 
    is washed out very strongly and results in too small asymmetry. ]
    
    
    § GAMMA-RAY AND ELECTRON/POSITRON FLUXES 
    
    
    
    As we have seen in the previous section, the dark neutron obtains a Majorana mass 
    when the portal operator is generated in association with the seesaw mechanism.  
    Due to the Majorana mass of the dark neutron, the dark neutron can oscillate into the anti-dark neutron. 
    The typical time scale of the oscillation, t_osc= m_M^-1, is estimated as
    
        t_osc   ≃3.3× 10^21  (Λ'_QCD/2 )^-6(M̃_C/3× 10^9 )^4(M_R/10^9 ) .
    
    We now see that some fraction of n' can convert into n̅' at  late time,
    and then n'/p' and n̅' annihilate into the dark pions. The neutral dark pions decay into the dark photons, and the dark photons finally decay into e^+e^- pairs. γ can be also emitted by the final state radiation (FSR) process as depicted in figure <ref>. 
    In this section, we discuss the constraints 
    on the late-time annihilation from the 
    observations of the γ-ray from the dSphs and the interstellar e^++e^- flux.
    
    
    
     §.§ Gamma-ray flux from the Dwarf Spheroidal Galaxies
    
    The γ-ray signal is one of the most promising channels
    to search for dark matter annihilation (e.g., <cit.> for review). 
    In particular, dSphs in our galaxy
    are the ideal targets to search for the γ-ray signal,
    since they have high dynamical mass-to-light ratios, (M/L ∼ 10-1000), 
    while they lack contaminating astrophysical γ-ray sources <cit.>.
    In this subsection, we estimate the γ-ray 
    fluxes from the dSphs and compare them with the upper limits 
    on the fluxes put by the Fermi-LAT.
    
    First, we calculate the γ-ray spectrum at production by the n'n̅' annihilation processes:
    
        n'n̅'̅→
        mπ'^0
        +
        lπ'^+ 
        + 
        lπ'^- 
        ,   
        (m,l = 0,1,2,⋯) .
    
    The cascade spectrum can be calculated by using the technique developed in <cit.>.
    
    We start to calculate the γ-ray spectrum at the rest frame of γ'. 
    For m_γ'≫ m_e, the spectrum is given by the Altarelli-Parisi approximation 
    formula <cit.>,
    [In the appendix <ref>, we compare the direct calculation of
    the FSR with the Altarelli-Parisi approximation formula, and confirm the validity of 
    the approximation in the parameter region we are interested in.]
        dÑ_γ/dx_0
        =
        α_EM/π1+(1-x_0)^2/x_0[-1+ln(4(1-x_0)/ϵ^2_0)] ,
    
    where ϵ_0=2m_e/m_γ' and
    x_0=2E_0/m_γ' with E_0 being the energy of γ at the rest frame of γ'. α_EM denotes the fine structure constant of SM QED.
    
    
    
    
    
    The next step is to translate the spectrum in the rest frame of γ' to that in the rest frame of π'^0. 
    For the case where m_π'≫ m_γ', the spectrum is calculated as
    
        dÑ_γ/dx_1
        =
        2∫^1_x_1dx_0/x_0dÑ_γ/dx_0
        f(2x_1/x_0-1)
        +
        𝒪(m^2_γ'/m^2_π') ,
    
    where x_1=2E_1/m_π' with E_1 being the energy of γ at the rest frame of π'. 
    The function f represents the effect of the anisotropy of the γ' decay. 
    According to <cit.>, we take
    
        f(cosθ)
        =
        3/8(1+cos^2θ) ,
     
    with θ being the angle between the γ emission line and the boost axis of γ'.
    Note that the angle θ is kinematically constrained as
    
        cosθ
        =
        2x_1/x_0-1
        +
        𝒪(m^2_γ'/m^2_π') .
    
    This is the reason why we put f(2x_1/x_0-1) in Eq. (<ref>).
    
    
    We next translate the spectrum Eq. (<ref>) to that in the center of mass (CM) frame for the ADM annihilation.
    In order to do that, we need to know how much π'^0 is boosted.
    If the total number of the dark pions is two (m+2l=2), 
    we can exactly know the energy/boost of the dark pions since they should be emitted back to back in the CM frame.
    In this case, the γ spectrum is calculated as
    
        dÑ^(m,l)_γ/dx_2
        =
        2∫^1_x_2dx_1/x_1dÑ_γ/dx_1
        +
        𝒪(m^2_π'/m^2_DM) ,
             m+2l=2 ,
    
    where x_2=E_2/m_DM with E_2 being the energy of γ at the CM frame.
     
    On the other hand, in the case of m+2l≥3, 
    it becomes highly non-trivial to know how much the π'^0 can be boosted 
    even when we assume that the matrix element of the annihilation is constant as a function 
    of the final state momenta.
    This is because, in this case, the energy spectrum of the dark pion is given as
    
        dN_π'/dξ
        =
        1/R_ndR_n/dξ ,
    
    where ξ=E_π'/m_DM and R_n is the n=m+2l body phase space integration <cit.>.
    E_π' denotes the energy of the dark pion in the CM frame.
    In general, it is difficult to perform the phase space integration for n≥ 3. 
    However, as discussed in <cit.>, under the assumption that m_π'^0=m_π'^+≡ m_π'≪ m_DM, we can perform the phase space integrations analytically as
    
        dN_π'/dξ
        =
        (n-1)(n-2)(1-ξ)^n-3ξ
        +
        𝒪(m^2_π'/m^2_DM),
    
    for n=m+2l≥3. Using the results, we finally obtain
    
        dÑ^(m,l)_γ/dx_2
        =
        2(n-1)(n-2)
        ∫^1_x_2dξ
        (1-ξ)^n-3∫^1_x_2/ξdx_1/x_1dÑ_γ/dx_1
        +
        𝒪(m^2_π'/m^2_DM)
         ,
    
    for n=m+2l ≥ 3 where we assume m_π'^0=m_π'^+≡ m_π'.
    
    Finally, we sum over the possible intermediate states
    and take into account the number of the final states.
    It turns out that the total γ spectrum from the n'n̅' annihilation is expressed as
    
        dN^(n'n̅')_γ/dx_2
        =
        ∑_m,l2m
        (
        ^(n'n̅')(m,l)
        dÑ^(m,l)_γ/dx_2)
         ,
    
    where ^(n'n̅')(m,l) denotes the branching ratio for the n'n̅'̅→ mπ'^0+lπ'^+ + lπ'^- annihilation process. The factor 2m corresponds to the number of e^+e^- pairs in the annihilation process.
    
    In the same way, we can estimate the γ spectrum from the p'n̅' annihilation processes:
    
        p'n̅'̅→
        mπ'^0
        +
        lπ'^+ 
        + 
        (l-1)π'^-
         ,   
        (m = 0,1,2,⋯, l=1,2,⋯) .
    
    The γ spectrum is calculated as
    
        dN^(p'n̅')_γ/dx_2
        =
        ∑_m,l2m
        (
        ^(p'n̅')(m,l)
        dÑ^(m,l)_γ/dx_2) ,
    
    with replacing n=m+2l by n=m+2l-1 in the calculation of dÑ^(m,l)_γ/dx_2.
    
    In the following analysis, we simply assume that the branching ratio of the dark nucleon annihilation can be estimated as that of nucleon-antinucleon annihilation. According to <cit.>, 
    we approximate the branching ratios by the fireball model,
    [In this approximation, the Parity violating mode, (m,l)=(2,0), is allowed, 
    although it is not significant numerically.]
        ^(n'n̅')(m,l)
           =
        2α^2l/(1+α)^n+(1-α)^n 
        _n C _2l 
        P_n ,    
        n=m+2l ,
        ^(p'n̅')(m,l)
           =
        2α^2l-1/(1+α)^n+(1-α)^n 
        _n C _2l-1 
        P_n ,    
        n=m+2l-1 ,
    
    where 
    
        P_n
        =
        1/√(2π)σexp(-(n-⟨n⟩)^2/2σ^2) ,
    
    with a=1/4, ⟨n⟩=5.05, σ^2=a⟨n⟩ and
    
        α
         =
           {[ √(2)     n=2 ,; 1.5     n≠2  . ].
    
    
    
    We are now ready to estimate the γ-ray spectrum emitted from the ADM annihilation. 
    Figure <ref> shows the value of the γ-ray spectrum. 
    Here, we take m_DM=10 GeV, m_π'=1 and m_γ'=40.
    The black solid and the dashed lines correspond 
    to the spectra predicted from the n'n̅' and p'n̅' annihilation, respectively.
    In the analysis, we ignore the contributions from the annihilation with large (m,l) since the branching ratios of them are much suppressed. 
    We stop taking the sum over (m,l) if the size of contribution 
    is less than 1% of the total amount.
    
    The figure shows that the ADM annihilation predict the continuous γ-ray spectrum peaked at the energy 
    of O(m_DM/10).
    This is expected as the typical number of the dark pions for an annihilation is five, and the neutral dark pion
    decays into two pairs of e^+e^-.
    
    It should be reminded that the γ-ray emission from the ADM annihilation can happen at the present universe since the ADM
    oscillation effectively happens at the late time scale.
    The ADM signals can  therefore be tested by γ-ray telescope experiments from nearby sources, 
    while evading the constraints from the observations of the cosmic microwave observations (see e.g. <cit.>).
    
    
    
    
    The γ-ray flux from the dSphs for an energy bin from E_min to  E_max is calculated as
    
        Φ
        =
        ∫^E_max_E_mindEE
        ∫_ΔΩdΩ/4π∫_l.o.s. dl
        (
        n_n'n_n̅'⟨σ v⟩_n'n̅'dN^(n'n̅')_γ/dE
        +
        n_p'n_n̅'⟨σ v⟩_p'n̅'dN^(p'n̅')_γ/dE) ,
    
    where we perform the integrations over
    a solid angle, ΔΩ, and the line-of-sight (l.o.s.).
    Here n_i and ⟨σ v⟩_ij denote the number density of a particle i at the dSphs
    and the kinematically averaged cross section for ij annihilation, respectively. N^(n'n̅')_γ and N^(p'n̅')_γ are the photon spectra from n'n̅' and p'n̅' annihilation which can be calculated from Eqs. (<ref>) and (<ref>), respectively.
    
    
    It should be noted that 
    the total amount of the γ-ray flux can be large enough to be tested by the γ-ray searches on the dSphs
    although the flux is suppressed by the factor,
    
        n_n̅'/n_n'   ≃(t_0/t_osc)^2
           ≃
        1.6× 10^-8(Λ'_QCD/2 )^12(M̃_C/3× 10^9 )^-8(M_R/10^9 )^-2 .
    
    where t_0≃ 4.3× 10^17 sec is the age of the universe. 
    This is because the thermally-averaged cross section can be large due to the strong interaction. 
    In the following analysis, we take the annihilation cross sections to be
    
        ⟨σ v⟩_n'n̅'
        =
        ⟨σ v⟩_p'n̅'
        =
        4π/m^2_DM ,
    
    to give rough estimation.
    Such a large annihilation cross section multiplied by the relative velocity is supported 
    by the cross section measurements of the non-relativistic nucleon and anti-nucleon 
    annihilation <cit.> 
    (see also <cit.>).
    [The cross section multiplied by the relative velocity in Eq. (<ref>) 
    is much smaller than the unitarity limit.
    
    In the Appendix <ref>, we discuss
    the Sommerfeld enhancement effects by the exchange of the dark pions.
    There, we find that the enhancement effects are not 
    significant in the present setup.
    ]
    
     In Figure <ref>, we show the predicted γ-ray flux from the Draco dSph. 
     The black solid, dashed, and dotted lines correspond to the γ-ray flux when we take M̃_C=10^9, 1.5× 10^9, and 2× 10^9, respectively. Here, we assume m_n'=m_p'=m_DM=10 and fix m_π'=1 and m_γ'=40. 
     To obtain the predicted γ-ray spectrum, we use the J-factors estimated in <cit.>
    which takes into account the effects of the non-sphericity of the dSphs.
    [As for the J-factor of the Ursa Minor classical dSphs, we use the value given in
    <cit.> as it is not analyzed in <cit.>.
    ]
    The green line corresponds to the upper bound (95% C.L.) on the γ-ray flux based on the  
    6 years of Pass 8 data by the Fermi-LAT collaboration <cit.>.
    The figure shows that the γ-ray flux from the 
    late-time annihilation becomes comparable to the 
    upper limit on the observed flux for M̃_C = 10^9 GeV and M_R = 10^10 GeV, which corresponds to the oscillation time scale of 
    t_osc=10^21sec.
    We discuss the constraints on the model parameters 
    by the Fermi-LAT in subsection <ref>.
    
    
    
    
     §.§ Interstellar Electron/Positron Flux
    
    The Fermi-LAT observation does not constrain 
    the late-time annihilation for m_DM≲ 3 GeV, 
    since the Fermi-LAT is sensitive to the γ-ray with 
    energy higher than 500 MeV.
    For such a rather light ADM, 
    the most stringent constraint is put by 
    the observation of the interstellar e^++e^- flux by the Voyager-1 <cit.> (see also <cit.>).
    In this subsection, we estimate the e^++e^- flux from the late-time annihilation in the Milky Way.
    
    The energy spectrum of e^++e^- at production 
    by the late-time ADM 
    annihilation is obtained by replacing dÑ_γ/dx_0 in Eq. (<ref>) with
    the e^+/e^- spectrum in the dark photon rest frame,
    
        dÑ_e/dx_0 = δ(x_0 - 1) .
    
    Here, x_0 = 2E_0/m_γ' with E_0 being 
    the energy of either e^- or e^+.
    By repeating the same analysis in the previous section,
    we can convert this spectrum to the one in the rest frame 
    of the ADM annihilation.
    In Figure <ref>, we show the e^+/e^- spectrum at production for 	m_e≪ m_γ'≪ m_π'≪ m_DM. 
    
    
    
    For a given e^+/e^- spectra at production, the interstellar e^++e^- flux at around the location of the Earth is 
    given by <cit.>,
    [A typical propagation time of the 
    cosmic ray to travel of 1 kpc is 
    much shorter than the 
    age of the universe.]
        dΦ/dE_e = 1/4π b(E)(ρ_DM/m_DM)^2
        (t_0/t_osc)^2
        ∑_i=n'n̅',p'n̅'⟨σ v ⟩_i∫_E^E_ DMdE_s 
        I(E,E_s) dN_e i/dE(E_s) .
    
    Here, ρ_DM denotes a local dark matter 
    density at around the location of the Earth, 
    I(E,E_s) is a Green function which encodes the propagation of 
    e^± from a source with a given energy E_s to 
    any energy E, and b(E) is the e^± energy loss function.
    [The Green function is dimensionless while  b(E) has a unit of GeV/sec
    which is typically 
    b(E)≃ 10^-(16-15) GeV/sec 
    for E=10 MeV to 1 GeV <cit.>.
    ]
    
    In Figure <ref>, we show the 
    interstellar e^++e^- flux at around the location of the Earth from the late-time ADM annihilation. 
    Here, the annihilation cross section and the oscillation time scale is set to be (t_0/t_osc)^2×⟨σ v⟩ = 1 pb.
    The Green function, I(E,E_s), and the energy loss rate, b(E), are those provided by <cit.>.
    In the figure, the solid lines assume the MED propagation model,
    while the upper and the lower dotted lines assume 
    the MAX and the MIN propagation models, respectively (see <cit.>).
    The dark matter profile is assumed to be the NFW profile <cit.>,
    [We numerically checked that the
    spectra are not significantly changed even 
    for a cored Burkert profile <cit.>,
    though they are slightly suppressed.]
    with the local dark matter density at around the Earth
    to be ρ_DM = 0.3 GeV/cm^3.
    
    In the figure, we also show the interstellar e^++e^- spectrum observed by the Voyager-1 <cit.>, where the data is taken from <cit.>. 
    The figure shows that the e^++e^- flux 
    from the late-time ADM annihilation is much smaller
    than the observed flux for (t_0/t_osc)^2×⟨σ v⟩ = 1 pb.
    We will summarize the constraints from the Voyger-1 
    in the next subsection.
    
    
    
     §.§ Constraints on Parameter Space
    
    
    
    As we have seen in the previous subsections, 
    we can probe the time scale of the matter-antimatter oscillation by the γ-ray observation up to t_osc=10^21 sec for m_DM≃ 10 GeV. 
    This oscillation time scale corresponds to the effective annihilation cross section,
    [The effective cross section into the γ-ray is further suppressed by Eq. (<ref>).]
        (t_0/t_ osc)^2 ⟨σ v⟩∼ 10  pb(10  GeV/m_ DM)^2
            (10^21sec/t_ osc)^2 .
    
    A lighter ADM can be also tested by the observation of the interstellar e^++e^- flux.
    
    
    In Figure <ref>, we show the constraints on the oscillation time scale from 
    the observations by the Fermi-LAT and the Voyager-1.
    Here, we assume m_e≪ m_γ'≪ m_π'≪ m_DM while we fix m_γ'=40 MeV.
    [The constraints do not depend on m_γ' significantly, as long as m_e≪ m_γ'≪ m_π'≪ m_DM.]
    The green region corresponds to the 95% C.L. excluded region from the Fermi-LAT 
    observations (see also <cit.>),
    where we take into account the γ-ray fluxes
    from the 8-classical dSphs.
    The yellow shaded region corresponds to the 95% C.L. excluded region from the Voyager-1 observation
    for the MED propagation model with the NFW dark halo profile.
    
    
    We see that, for m_DM≃ 5 10, the more stringent 
    constraints are put by the Fermi-LAT observation, where the 
    oscillation time scale shorter than t_osc∼ 10^21 is excluded.
    For a lighter mass region,
    the Voyager-1 observation excludes the oscillation time scale shorter than 
    t_osc∼ 10^21-22 sec.
    
    
    
    
    In Figure <ref>,
    we translate the constraints on the oscillation time scale to those on the 
    parameters of the present model.
    In the figure, we consider m_DM=2,5,10 GeV.
    We also take  Λ'_QCD=2 ×(m_DM/10 GeV)
    to mimic QCD for each choice of the dark matter mass.
    We also assume  m_e ≪ m_γ'≪ m_π'≪ m_DM.
    The green and yellow shaded regions correspond the 95% C.L. excluded regions by 
    the Fermi-LAT and the Voyager-1, respectively. 
    The lower gray region is excluded 
    where the B-L asymmetry is washed out (see Eq. (<ref>)).
    Above the solid line, we require an on-shell B-L portal sector (see Eq. (<ref>)).
    
    
    We now see that the composite ADM scenario 
    with M_C = 𝒪(10^9) GeV can be tested by the γ-ray searches from the dSphs by the Fermi-LAT for m_ DM≃ 10 GeV.
    Even for a lighter ADM scenario, we see that the region with M_C = 10^8 GeV
    has been excluded by the Voyager-1 observation.
    The resultant constraint is important in view of the fact that the parameter region with M_C∼ 10^9–10^10 GeV
    is highly motivated in the UV completion 
    model based on SU(4)  (⊃ SU(3)_D× U(1)_D) gauge theory <cit.>.
    In this UV completion, the tiny kinetic mixing of ϵ = 10^-10–10^-9 which evades all the phenomenological constraints on the dark photon <cit.>
    is achieved when the SU(4) breaking scale 
    is at around 10^9–10^10 GeV.
    The SU(4) breaking scale also leads 
    to the colored dark Higgs mass in a similar range.
    The γ-ray searches are already sensitive to such a well-motivated parameter region 
    for m_DM≃ 10 GeV.
    
    
    Several comments are in order.
    In our discussion, we consider only the γ-ray emitted by the FSR. 
    This should be justified as the γ-rays made by the Synchrotron radiation and the inverse Compton scattering 
    from the sub-GeV e^+/e^- are very soft and below the Fermi-LAT sensitivity <cit.>.
    It should be also noted that the γ-ray signal from the galactic center does not lead to more stringent constraints,
    despite the signal strength is higher than that from the dSphs.
    This is because the γ-ray background is much higher 
    for the galactic center, and hence, it is difficult to 
    distinguish the continuous signal spectrum from the background spectrum.
    
    Future γ-ray searches 
    such as 
    e-ASTROGAM <cit.>,
    SMILE <cit.>,
    GRAINE <cit.>,
    and 
    GRAMS <cit.>
    projects will be important to test the model further.
    It should be emphasized that those experiments 
    are sensitive to the MeV γ-rays, and hence,
    they are also able to test the models with m_DM= a few GeV to which the Fermi-LAT loses 
    sensitivity.
    In Figure <ref>, we show the prospected lower limit on t_osc at 95%CL 
    by the  γ-ray search 
    from the dSphs by e-ASTROGAM in one year of effective exposure.
    In our analysis, we used the effective area and the prospected sensitivities for a γ-ray flux
    from a point-like source at a high latitude (in Galactic coordinates) in <cit.>.
    The testable parameter region can be wider when the J-factors of the 
    ultra-faint dSphs are determined more precisely by 
    future spectroscopic observations such as the Prime Focus Spectrograph <cit.>.
    For example, if the J-factor of Triangulum II converges to the central value in <cit.>,
    i.e. log_10 J ≃ 20, the prospected lower limit on t_osc
    becomes higher for about a factor of 2^1/2.
    
    
    
    
    
    
    § CONCLUSIONS
    
    The composite ADM model is particularly motivated as it provides the DM mass of O(1) GeV 
    and a large annihilation cross section simultaneously.
    In this paper, we discussed the indirect detection of the composite ADM where the portal operators 
    of the B-L asymmetry is generated in association with the seesaw mechanism.
    In this model, the dark-neutron obtains a tiny Majorana mass,
    and hence, ADM can pair-annihilate at later times.
    
    As we have discussed, the late time annihilation of ADM results in multiple soft electrons/positrons and soft photons emitted as the FSR.
    As a result, some parameter region of the composite ADM which 
    is motivated by thermal leptogenesis and dark UV completion models
    has been excluded by the Fermi-LAT and the Voyager-1 observations.
    The obtained constraint is tighter than that from  the anti-neutrino flux  made by the 
    decay of  ADM via the B-L portal operator <cit.> (see Eq. (<ref>)).
    Future experiments which are 
    sensitive to sub-GeV γ-rays such as 
    e-ASTROGAM <cit.>,
    SMILE <cit.>,
    GRAINE <cit.>,
    and 
    GRAMS <cit.>
    projects will be important to test the oscillating ADM model
    further.
    
    
    
    
    § ACKNOWLEDGEMENTS
    
    
    MI thanks K. Kohri for useful comments on the possible constraints on the model from the cosmic electron/positron rays.
    This work is supported by JSPS KAKENHI Grant Numbers, 19K14701 (R. N.),  No. 15H05889, 
    No. 16H03991, NO. 17H02878 
    and No. 18H05542 (M. I.) and by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. 
    This work is also supported by the Advanced Leading Graduate Course for Photon Science (S. K.).
    
    
    § FINAL STATE RADIATION IN THE DARK PHOTON DECAY
    
    
    
    This appendix is devoted to the photon energy spectrum of the final state radiation in the 
    dark photon decay, γ' → e^+e^-γ.
    One of the diagrams is shown in the figure  <ref>.
    
    The invariant amplitude for this process is 
    
        ℳ=    - 4πϵα_EMu̅(p_1) [ ε^*(p_3) p_1+p_3+m_e/(p_1+p_3)^2-m_e^2ε(p_0) + ε(p_0) -p_2-p_3+m_e/(p_2+p_3)^2-m_e^2ε^*(p_3) ] v(p_2)   ,
    
    where ϵ represents the strength of kinetic mixing, α_EM the fine structure constant
    of QED, ε the polarization vector, m_e the electron mass, u and v spinors and p momentum vector.
    Here the subscripts (0, 1, 2, 3) denote the (γ', e^-, e^+, γ).
    
    Summing over the spins of the final state e^-,  e^+ and averaging over the helicity of initial state γ', we obtain
    
        1/3∑_spinℳ^2=   8 (4πϵα_EM)^2/31/(m_13^2-m_e^2)^2 (m_23^2-m_e^2)^2[ m_13^2 m_23^2 { 2m_12^4 + 2m_12^2 (m_13^2+m_23^2) +m_13^4 + m_23^4 }. 
            - m_e^2 (m_13^2+m_23^2) { 2m_12^4+4 m_12^2 (m_13^2+m_23^2) +3 (m_13^2+m_23^2)^2 }
            + m_e^4 { 2m_12^4 + 10 m_12^2 (m_13^2+m_23^2) + 11 (m_13^2+m_23^2)^2 }
           . - 4 m_e^6 { 2m_12^2 + 3 (m_13^2+m_23^2 ) } + 2 m_e^8 ]  ,
    
    by using the Mandelstam invariants, m_ij^2 =(p_i-p_j)^2, with the subscripts defined above.
    There is a relation between the invariants,
    m_γ'^2 + 2 m_e^2 = m_12^2+m_13^2+m_23^2, with m_γ' being
    the dark photon mass.
    This expression is symmetric under the exchange between m_13^2 and m_23^2 as expected.
    
    Now, let us calculate the decay rate with the final state radiation.
    In the following calculation, we use the center of mass frame in which three out-going particles lie in a same plane. 
    Thus, we can transform the three-body phase space integral into integration over the energy of two particles and three angles. 
    By taking into account of the energy-momentum conservation, the three-body phase space has 9-4=5 d.o.f.
    After fixing the energy of e^-, three d.o.f. remain.
    Two of them are angles (α, β) that specify the direction of p⃗_3.
    The last one is an angle δ which determines the plane of decay around p⃗_3.
    Thus, Γ_γ'→ e^+e^-γ can be written as
    
        Γ_γ'→ e^+e^-γ   =∫1/16m_γ'1/3∑_spinℳ^2dE_3 dE_1 dα d(cosβ) dδ/(2π)^5 , 
           =m_γ'/32(2π)^3∫1/3∑_spinℳ^2 dx dy   , 
        
        																	   =m_γ'/32(2π)^3∫ dx ∑_n=0^2ϵ_0^2n[ f_n(x,y_max(x,ϵ_0)) - f_n(x,y_min(x,ϵ_0)) ]  .
    
    Here we define x=2E_3/m_γ',  y=2E_1/m_γ' and ϵ_0 = 2 m_e/m_γ'.
    Each f_n(x,y) is defined as the integration of the invariant scattering amplitude over E_1, i.e., y.
    The analytical formula for each f_n(x,y) is as follows:
    		
        f_0(x,y)=8/3 (4πϵα_EM)^2 [2(1-y)-1+(1-x)^2/xln(1-y/1-x-y)]  , 
        
        			   f_1(x,y)=4/3 (4πϵα_EM)^2 [x+2y-2/(1-y)(1-x-y)+2ln(1-y/1-x-y)]  , 
        
        			   f_2(x,y)=2/3 (4πϵα_EM)^2 [x+2y-2/(1-y)(1-x-y)+2/xln(1-y/1-x-y)]  .
    
    Here y_min and y_max are the lower and the upper bounds of the integration region of y corresponding to the Dalitz region.
    The explicit forms of y_min and y_max are
    		
        y_min=max[1-x/2-x/2√(1-ϵ_0^2/1-x),ϵ_0]  , 
        
        			   y_max=min[1-x/2+x/2√(1-ϵ_0^2/1-x),1]  .
    
    	
    
    From above, we obtain the energy spectrum of the final state radiation photon.
    The energy spectrum is expressed as <cit.>
        1/N_γN_γx    = 1/Γ_γ'→ e^+e^-d Γ_γ'→ e^+e^-γ/dx .
    
    	Here, Γ_γ'→ e^+e^- = 1/3ϵ^2 α_EM m_γ' is the decay rate of the process γ'→ e^+e^-.
    We compare the result with twice the Altarelli-Parisi approximation formula <cit.>
        1/Γ_γ'→ e^+e^-d Γ_γ'→ e^+e^-γ/dx= α_EM/π1+(1-x)^2/xln(4(1-x)/ϵ_0^2)  ,
    
    in the figure <ref>.
    We take m_γ' = 40 MeV.
    We see that two formulae are in good agreement in a wide range of the photon momentum.
    	
    
    
    
    § SOMMERFELD ENHANCEMENT
    
    The dark pion exchange 
    between the dark nucleons 
    generates attractive/repulsive forces between 
    them depending on their spins and the isospins.
    [Since the dark quark masses are assumed to be
    much smaller than the dark dynamical scale, 
    the dark sector possesses the isospin symmetry 
    as in the case of the QCD in the SM sector.]
    For example, one dark pion exchange results in a
    static potential,
    
        V(𝐫) = g_A^' 2/16π f_π'^2 (τ_1·τ_2)
            (σ_1·∂)
            (σ_2·∂)
            1/re^-m_π'r ,
    
    which goes like 1/r^3 in the region of r≪ m_π'^-1.
    This potential is obtained from the axial-current interaction,
    
        L =g_A'/f_π' ∂^a
        π' N̅'γ_μγ_5(τ^a/2)N' ,
    
    where f_π' is the decay constant of
    the dark pion and g'_A is the form factor of the dark nucleon axial current.
    [We take the normalization such that 
    f_π≃ 93 MeV and g_A≃ 1.26
    in the case of the SM.
    ]
    The spin and the isospin indices are implicit,
    where σ and τ denote the Pauli matrices
    applying to the spin and the isospin 
    of each nucleon, respectively.
    The way of the isospin transition can be read off by noting τ_1 ij·τ_2,kℓ =
    2 (δ_iℓδ_jk-δ_ijδ_kℓ/2).
    
    As discussed in <cit.>, the attractive potential forces mediated by the pseudo-scalar field causes 
    the Sommerfeld enhancement of the dark matter annihilation <cit.>.
    In this appendix, we discuss the Sommerfeld enhancement caused by the dark pion exchange.
    In our analysis, we rely on the formalism of the Sommerfeld enhancement in <cit.>, in which the lower cut-off on the relative velocity is taken into account in a self-consistent way.
    
    Following <cit.>,
    we approximate the potential by a spherical one,
    
        V(𝐫) ≃ -g_A^' 2/16π f_π'^21/r^3e^-m_π' ,
    
    and estimate
    the enhancemnt of the s-wave annihilation.
    [Strictly speaking, 
    we need to solve a coupled equation 
    between the states with angular momenta,
    since the potential force in Eq. (<ref>) changes the nucleon angular momentum by Δℓ =±2.
    ] 
    Under this approximation, the Sommerfeld enhancement factor can be obtained by solving the effective Schrödinger equation,
    
        [-∇^2/2m_RED+ V(𝐫) + u δ^(3)(𝐫)]ψ(𝐫) = p^2/2m_REDψ(𝐫) .
    
    Here, m_RED = m_DM/2 is the reduced mass and p denotes relative momentum of the incident dark matter.
    The boundary condition of the wave function ψ(𝐫) is taken to be
    an incident plane wave with an outgoing spherical wave, i.e. ψ(𝐫)→ e^ipz+f e^ipr/r at r→∞.
    The complex parameter u encodes the annihilation cross section at a short distance without the Sommerfeld enhancement factor, 
    i.e. u = - i σ v_0/2.
    [The dark-nucleon self-scattering due to short-range forces can be also encoded in the real part of u. In our analysis, 
    we assume the self-scattering by short-range forces are subdominant and take u ≃ 0. ]
     
    Since the potential goes to infinity faster than r^-2 at the origin, it must be regularized at short distances.
    In our analysis, we introduce a short distance cutoff r_0 
    satisfying  V(r_0) = m_DM and regulate the scalar potential by replacing V(r)→ V_reg(r) = V(r+r_0) <cit.>.
    [Our conclusions do not depend on the choice of the regularization significantly.]
    With the regulated potential, the Sommerfeld enhancement factor 
    is given by <cit.>,
    
        S_ENF(v) = σ v/σ v_0≃S(v)/|1 - i m_RED^2/4πσ v_0(T(v)+i S(v))v|^2 .
    
    Here, T(v) and S(v) are given by,
    
        T(v) = 1/p(.dg_p/dr|_r=0
            - 
           .dg_p_0/dr|_r=0) , 
           S(v) = 1/p.d g_p/dr|_r=0 ,
    
    with the function g_p(r) being a solution of 
    
        [
            -d^2/dr^2
            +2 m_RED V_reg(r) - p^2
            ]    g_p(r) = 0 ,
           g_p(0) = 1 ,
        lim_r→∞   g_p(r)∝ e^ipr .
    
    The short distance cross section σ v_0 is 
    fixed at a high momentum p_0.
    
    In Eq. (<ref>), the factor S(v) corresponds to the naive Sommerfeld enhancement factor.
    The denominator, on the other hand, provides an IR cutoff in the limit of v→ 0 with which the unitarity violation by the naive Sommerfeld enhancement factor is regulated self-consistently.
    The regularization effect is particularly important 
    when the short-distance cross section is large as in the case of the ADM scenario.
    In Figure <ref>,
    we compare the naive enhancement factor shown in <cit.> 
    and the one in Eq. (<ref>) 
    by assuming σ v_0 = 4π/m_DM^2.
    [Due to a slightly different choice of V_reg(r), 
    the positions of the resonances appearing in S(v) are shifted from those in <cit.>.]
    
    The figure shows that the enhancement factors 
    at the resonances are significantly suppressed 
    when the short-distance annihilation cross section is large.
    
    
    
    
    
    
    
    
    Now, let us apply Eq. (<ref>) to the 
    dark nucleon annihilation.
    In Figure <ref>, we show the Sommerfeld enhancement factor as a function of m_DM for g_A = 1, f_π' = 1 GeV, and m_π'=1 GeV.
    
    
    The figure shows that the regularization effects
    are important at around the resonance, m_DM≃ 21 GeV.
    The figure also shows that the Sommerfeld enhancement factor for the mass region of the ADM, 
    m_DM≲ 10 GeV, is less significant.
    
    As we fix the short-range cross section of the ADM, σ v_0 ≃ 4π/m_DM^2,
    to mimic the measured nucleon annihilation cross section at v = 10^-1 <cit.>,
    the effective Sommerfeld enhancement factor corresponds to S_ENF(v)/S_ENF(10^-1).
    The figure shows that the effective enhancement factor is close to unity for m_DM≲ 10 GeV.
    
    In Figure <ref>, we also show the 
    Sommerfeld enhancement factor for more realistic relations between the parameters,
    
        f_π' = 0.1× m_ DM , 
               
            m_π' = 0.1× m_ DM ,
    
    which mimic QCD.
    
    
    The figure shows that no resonance appears 
    when the parameters satisfy these relations.
    As a result, we find that the effective enhancement factor, S_ENF(v)/S_ENF(10^-1), is of 1.
    [The Sommerfeld enhancement
    for coupled channels between different 
    angular momenta requires more careful analysis.
    However, as the centrifugal barriers of the higher angular momenta make the attractive potential wells shallower 
    and smaller in spatial size, 
    the resonances are expected to appear at a 
    higher dark nucleon mass than those 
    for ℓ = 0.
    Thus, the coupled equations do not lead to resonances in the mass range m_DM≲ 10 GeV.]
    We also numerically confirmed that the results do not depend on the dark pion mass as long as it is much lighter than the dark nucleon.
    Therefore, we conclude that the Sommerfeld enhancement is not significant in the present setup.
