Cosmic-Ray Signatures of Dark Matter from a Flavor Dependent Gauge Symmetry Model with Neutrino Mass Mechanism

We propose an extension to the Standard Model accommodating two families of Dirac neutral fermions and Majorana fermions under additional ${U(1)_{e-\mu} \times Z_3\times Z_2}$ symmetries where ${U(1)_{e-\mu}}$ is a flavor dependent gauge symmetry related to the first and second family of the lepton sector, which features a two-loop induced neutrino mass model. The two families are favored by minimally reproducing the current neutrino oscillation data and two mass difference squares and canceling the gauge anomalies at the same time. As a result, we have a prediction for neutrino masses. The lightest Dirac neutral fermion is a dark matter candidate with tree-level interaction restricted to electron, muon and neutrinos, which makes it difficult to detect in direct dark matter search as well as indirect search focusing on the ${\tau}$-channel, such as through ${\gamma}$-rays. It may however be probed by search for dark matter signatures in electron and positron cosmic rays, and allows interpretation of a structure appearing in the CALET electron+positron spectrum around 350-400 GeV as its signature, with a boost factor $\approx$40 Breit-Wigner enhancement of the annihilation cross section.


I. INTRODUCTION
The cosmological standard model includes dark matter (DM) as an essential component, commonly considered to be a neutral particle not part of the standard model of particle physics (SM).
Assuming thermal production in the early Universe, a weakly interacting massive particle (WIMP) in the GeV-TeV mass is a strong candidate, since the Weak Interaction of the SM yields just the right annihilation cross section to predict the observed relic density of DM, a relation known as the WIMP miracle. This default candidate is the main target of experimental DM search, and since the weak interaction couples universally to all leptons and quarks, its parameter space is successively scanned and ruled out by direct detection experiments based on WIMP-nucleon interactions [1][2][3][4][5] and indirect searches looking for the products of annihilation into hadronic channels, such as anti-protons [6,7] and γ-rays [8].
Avoiding hadronic interaction of DM requires the introduction of a new force and corresponding charge, which is only carried by the DM and leptons. In the initial version of this Leptophilic Dark Matter [9], all lepton generations carry the same charge, resulting in equal branching ratios in the annihilation of DM. In this case, the strongest constraints on the DM annihilation cross section come from observation of dwarf galaxies in γ-rays based on the DM + DM → τ + + τ − channel, which due to its higher γ-ray multiplicity yields limits about half a magnitude more strict than those on DM + DM → e + + e − channel and DM + DM → µ + + µ − channel [10,11]. These limits are subject to about one order of magnitude variation from uncertainty on the halo shape and resulting J-factors [12], which however is independent of the annihilation channel. Most recently very strict limits on hadronic and the DM + DM → τ + + τ − channel based on the morphology of γ-ray flux from the galactic center have been brought forward [13], giving explicitly no such constraint on DM + DM → e + + e − channel and DM + DM → µ + + µ − channel.
On the other hand, search for DM annihilation in positron and electron cosmic rays with detectors such as AMS-02 [14][15][16][17][18][19], CALET [20,21] , DAMPE [22] and the Cosmic Ray Subsystem on the Voyager probes is most sensitive to the electron channel, since its signature is a sharp drop in the spectrum at the mass of the DM particle which can be recognized above a smooth astrophysical background [23][24][25]. For GeV-TeV range DM, the target region is the local DM halo within ∼kpc range due to the energy loss and resulting limited propagation distance of electron cosmic rays. This complementarity can reduce the possible impact of astrophysical uncertainties in the case of DM with universal coupling to leptons. For DM with selective coupling to the different lepton flavors, either search with γ-rays or charged cosmic rays may have preferential sensitivity.
Apart from DM, the other strong indication of physics beyond the SM is the neutrino mass, and many theoretical models extending the SM aim at solving both issues simultaneously, examples being radiative seesaw models at one-loop [26], two-loop [27,28], and three-loop [29][30][31]. Several models extending the SM by an additional U(1) gauge symmetry have been proposed, which favor annihilation or decay to tau and/or muon as a possible DM-only explanation of the positron excess [32][33][34], while also featuring a mechanism for giving the neutrinos mass.
In this context we investigate if a thermally produced DM candidate based on a flavor-specific U (1) e−µ gauge symmetry coupling only to electron and muon is also feasible, corroborated by simultaneous explanation of the neutrino sector. This kind of DM would be a favorable target to search in electron-positron cosmic rays while being less detectable by γ-ray search. After establishing the particle physics model defining the properties of the DM, we discuss its cosmic-ray signatures and implications from available CALET and AMS-02 data. While introduction of a new flavor-specific gauge interaction lacks the elegance of the classical WIMP, studying such a model seems worthwhile as it allows to keep a thermal production mechanism and a WIMP-like DM candidate. This should be seen against the trend of DM candidate theory becoming more and more diversified to avoid constraints on the WIMP and WIMP-like particles [35].
Our extension of the SM is based on a radiatively induced neutrino mass (scotogenic model), which originally provides us with an appropriate explanation of the hierarchy among the Yukawa sector of the SM. The ratio between the top Yukawa quark coupling(∼ 1) and the electron Yukawa coupling(∼ 10 −6 ) is of the order 10 6 , which respectively are the heaviest and lightest masses in the fermion sector of the SM. However, the ratio between the electron Yukawa coupling and the typical neutrino Yukawa coupling(∼ 10 −13 ) is of the order 10 7 . If we assume the neutrino mass to be of Dirac type and to be induced at tree level, which is the same as for the other matter sectors in the SM, this would suggest that there is a huge gap between the neutrino coupling and the other three Yukawa couplings. The scotogenic model generates neutrino mass at loop level, with newly introduced fields running inside the loop. It is found that with a 0.01 loop suppression factor and two Yukawa couplings at one-loop level in the neutrino mass formula, the order of Yukawa coupling at one-loop level is minimally 10 −6 , which is comparable to the electron Yukawa coupling. We fix the mass scale of one new field to be on the order of one TeV, which allows for the new scale to be tested by current experiments. Another advantage of this model is its predicted correlation between the DM candidate properties and the neutrino mass, since the DM field is running inside the neutrino loop. Therefore, the neutrino interacts with SM-like Higgs only though the DM field in the generation of the neutrino mass. This provides a natural explanation for the tininess of the neutrino mass, and phenomenology apart from direct and indirect DM search. Further phenomenology arises from the not so small strength of the Yukawa couplings and their nontrivial structure to induce the neutrino mixings as well as mass eigenvalues, which might cause lepton flavor violations (LFVs) that are severely restricted by current experiments such as MEG [36,37].
To realize a sufficiently high cross section yielding the observed relic density in thermal production of the DM candidate, the annihilation process should be s-channel dominated, which however is helicity-suppressed for a Majorana particle. Therefore we introduce a Z 3 discrete symmetry under which the DM is charged, giving it a Dirac nature and ensuring its stability. Also, we impose a Z 2 discrete symmetry to forbid tree level neutrino mass, where this symmetry is softly broken in the Higgs potential and its broken term contributes to generating the tiny neutrino mass. The neutrino mass is induced at two-loop level, where we introduce two types of neutral fermions; Dirac type and Majorana type. In the neutrino sector, the U (1) e−µ symmetry also plays an important role in predicting the neutrino mass. Because the nonzero charges (-1 or 1) have to be assigned to only two families, the minimal number of new fermions are two families, which is also the minimal number to explain the active neutrino oscillation data and their mass eigenvalues. Furthermore, the two families are required to allow gauge anomaly cancellation in a minimal manner. Thus, we predict one massless neutrino that causes the other two massive neutrinos to be uniquely determined by the experimental results, which are the squared solar mass difference and squared atmospheric mass difference, as we will discuss for both cases of normal hierarchy (NH) and inverted hierarchy (IH) in detail. This paper is organized as follows: In Section II, we explain our particle physics scenario and formulate the lepton sector and the Higgs potential, the masses and mixings for the two new fermions and the active neutrinos, and the mass of the new gauge boson and its interactions, also discussing LFVs. In Section III, we will discuss our DM candidate, in which we briefly explain why it is not subject to current bounds from direct detection searches, and explain calculation of the DM relic density. We also show that Breit-Wigner enhancement may lead to a significant boost factor (B) on the annihilation cross section, which may increase the signatures to the level detectable by current indirect DM search. In Section IV, the electron and positron cosmic-ray signature of the DM candidate X is explained, and after introduction of propagation and astrophysical background models, the e − + e + flux measured by CALET [21] and the e + -only flux measured by AMS-02 [19] are interpreted including the DM signature. It is shown that step-like spectral structures in the CALET spectrum could be identified with the signature of the DM candidate, identifying the best-fit regions in m X vs. B space. Finally we summarize and conclude our results in Section V.

II. PARTICLE PHYSICS MODEL
In this section, we review our model [38]. At first, we explain our motivation for introducing new symmetries and fields. Then, we construct the Lagrangian and Higgs potential, and continue with formulating the neutral fermions, LFVs, and the additional gauge boson sector.

A. Particle Contents and Lagrangian
We introduce three families of vector-like fermions (N e , N µ , N τ ), and two families of Majorana fermions (ν Re , ν Rµ ) in the fermion sector, so that we can construct a two-loop induced neutrino mass model. These fermions are minimally required to reproduce the neutrino oscillation data and cancel the anomaly for ν R . We extend the scalar sector by introducing an isospin doublet inert boson η, an isospin singlet inert boson S, and a singlet boson ϕ that gives nonzero Vacuum Expectation Values (VEVs) to spontaneously break the U (1) e−µ symmetry as shown later, where the SM-like scalar boson is symbolized by H. Here we denote their VEVs as ϕ ≡ v ϕ / √ 2 and H ≡ v H / √ 2, respectively. In addition, we impose three additional symmetries; gauged symmetry U (1) e−µ and discrete Abelian symmetries Z 3 and Z 2 . The first symmetry defines the newly introduced interaction with only the two first generations of leptons, giving the model the intended property of avoiding gauge interactions with the τ -lepton, while the second one provides stability of potential DM candidates N, η, S, and assures the Dirac feature of N . We associate the lightest Dirac particle N with DM, since the heavier ones can decay into the lighter ones via five-dimensional terms even though the decay is forbidden within the renormalizable theory. The field contents and their assignments are summarized in Table I for fermions and Table II where Z 2 is softly broken and all the fields are singlet under SU (3) C .  Yukawa Interactions: Under these fields and symmetries, the renormalizable Lagrangian for quark and lepton sector is given by whereη is defined by iσ 2 η * , σ 2 being the second Pauli matrix, andN C R ν R S * is also allowed by our symmetries but it does not contribute to any phenomenologies. Thus we neglect this term. Z 2 symmetry forbids the Dirac termL LH ν R at tree level, where Z 2 is softly broken at the Higgs potential below.
where the µH † ηS * term is softly broken under Z 2 symmetry, and we expect µ to be of a rather small scale compared to the electroweak scale. We parametrize the scalar fields as where η 0 and S are complex scalars, v H 246 GeV is VEV of the SM Higgs, and w ± , z, and z are GeV is the mass of the SM Higgs. Here we define the mixing matrices as where c(s) a(α) is the short-hand notation of cos(sin) a(α) . While values s a > 0 could be chosen within experimental limits, we take s a = 0 in our numerical analysis for convenience as shown later.
Neutral Dirac Fermions: After the e − µ spontaneous breaking, the Dirac fermion mass matrix in basis of [N e , N µ , N τ ] T is found as: where M 1,2,3 is the mass eigenstate, and ψ is the mass eigenvector of N .

Neutral Majorana Heavier Fermions:
In a way similar to the Dirac fermion, the Majorana fermion mass matrix in the basis of [ν Re , ν Rµ ] T is found as: where M R 1,2 = M Neµ is the mass eigenstate, and Ψ is the mass eigenvector of ν R .

B. Active Neutrino Mass
The dominant contribution to the active neutrino mass matrix arises from the canonical seesaw model, but the Dirac mass matrix m D is given at one-loop level. Thus the neutrino mass is induced at two-loop level. Before formulating the neutrino sector, we evaluate the number of complex parameters. First of all, three components of y η can be real by phase redefinition for L Le,µ,τ , which implies that the phases of N Re,µ,τ and e R , µ R , τ R are fixed. Second, the two components of y S can also be real by the redefinition for ν Re,µ , which suggests that the phases of N Le,µ are fixed. Finally, one phase in M N can be real by the phase redefinition for N Lτ . Here we identify M τ to be real.
Thus, we have six phases in M N . The canonical seesaw is given by the following form: where m D is found as follows [39,40]: The neutrino mass matrix is then diagonalized by a unitary matrix U ν as Here we can identify U ν as the PMNS matrix [41] because of the diagonal mass matrix of the charged leptons, which is achieved by the U (1) e−µ gauge symmetry. Each of the mixings is then given by: (II. 13) In case of NH, we find that the neutrino mass eigenvalues and the effective neutrinoless double beta decay m ee are respectively given in terms of observables and phases as m ee ∆m sol sin 2 θ 12 cos 2 θ 13 e iα 21 + ∆m 2 atm sin 2 θ 13 e i(−2δ CP ) , (II. 15) where ∆m 2 atm and ∆m 2 sol are respectively atmospheric mass difference square and solar mass difference square which are observables [42]; therefore these three neutrino mass eigenvalues are uniquely determined. Here, we redefine the neutrino mass eigenstate as Then, s α (c α ) can be rewritten by which implies that s α is determined by the two parameters ∆m 2 atm and |m 3 | 2 . Also, ∆m 2 atm is fixed by Similar to the case of NH, we also find the neutrino mass eigenvalues and m ee in case of IH to be And s α and ∆m 2 sol are found by Here, we redefine the neutrino mass eigenstate as

C. Lepton Flavor Violations
Lepton Flavor Violations (LFVs) arise from the term y η at one-loop level, and their branching ratios are given by and these bounds give constraints on the related Yukawa couplings and masses in the loop. It is worthwhile to mention the muon anomalous magnetic moment ∆a µ . Although we have a new contribution to ∆a µ from the same term as LFVs, its sign is negative, which is opposite to the experimental result. Thus, we assume a different effect to dominantly cause the anomaly and do not discuss it further.

D. Z e−µ Gauge Boson
After the U (1) e−µ symmetry breaking, we find the massive Z e−µ gauge boson that is denoted by Z hereafter, and its mass m Z is given by where g is the gauge coupling of the U (1) e−µ symmetry and we neglect kinetic mixing for simplicity.
Gauge interactions among Z are given by Then each of the decay rates of Z is given by where we have assumed 2M 1 < m Z < M R 1,2 , M 2 , M 3 and M 1 is considered to be the DM in the next section. When the decay rate of Γ(Z → XX) can be negligible, the branching ratios are respectively found as Since Z couples to an electron and positron pair, we have to impose the following constraint which comes from LEP [44]: where we have adopted a conservative bound. Here, we briefly mention other possibilities to detect signatures at colliders in the future, for the case of the Z -mass being of the order of 100 GeV. First is the Large Hadron Collider (LHC), which can observe the mode qq . Second is the future International Linear Collider (ILC) which 46,47]. So far there is no analysis of LHC data for above channels, thus LHC provides no constraint on the model parameters.

III. PROPERTIES OF DARK MATTER PARTICLES
The DM candidate in this model is the lightest Dirac fermion X ≡ ψ 1 , and its mass given by m X ≡ M 1 . In this section we study with which model parameters DM consisting of X and X is viable, with the goal of showing the existence of an allowed region, leaving a complete scan of the whole possible parameter space for future work. First, we briefly discuss detectability by direct detection searches and the reason why we take s a = 0. Then, we explain the calculation of the DM relic density which is determined by gauge interaction via s-channel, and perform a numerical analysis to explore the region around the pole m X = m Z /2 which satisfies all discussed constraints. Finally, we discuss that by applying Breit-Wigner enhancement to our model, the annihilation cross section in the current Universe can be increased by a boost factor B compared to a generic thermally produced DM with velocity independent annihilation cross section.

A. Direct Detection
The latest bound on spin-independent scattering is reported by the XENON1T experiment, which gives an upper limit on the spin independent elastic DM-nucleon cross section σ: σ < 4.1 × 10 −47 cm 2 at m X = 30 GeV with 90% confidence level [4]. Our DM dominantly interacts with nuclei only via the mixing of s a at tree level arising from the terms f ( ) ϕ . Then, our scattering cross section is given by where µ nX ≡ m n m X /(m n + m X ), m n is the mass of neutron, and C ≈ 0.3 is determined by lattice simulation among DM and nucleon. The easiest way to evade this constraint is to assume s a to be zero. LHC results also favor s a to be small with an upper bound of s a 0.2 [48]. Thus, with the choice of s a = 0, no direct detection bounds need to be considered for our DM candidate. We leave exploring possible bounds for the case s a > 0 for future work.

B. Relic Density
In the following we discuss the relic density of DM. The relevant processes arise from Yukawa interaction via y η and kinetic interaction via g . In case of Yukawa interaction, the coupling is mainly restricted by µ → eγ, which is typically of the order 0.01, although it depends on its flavor structure. Then the cross section via y η is 10 −17 GeV −2 at most. Thus, Yukawa contribution is negligibly small compared to the cross section ∼ 10 −9 GeV −2 required to explain the relic density.
As a result, the dominant cross section to the relic density comes from kinetic interaction.
We make use of the micrOMEGAs package [49] to calculate the speed averaged cross section σv rel , and the relic density. micrOMEGAs is adapted to this model by defining the properties of the interaction mediated by Z in the form of a kinetic term simplified from Eq. (II.25) as follows: The model parameter space is thus effectively given by m X , m Z , g and a x , with a x taking values in the interval [0, 1].
The evolution of the DM abundance is given through the Boltzmann equation is the entropy density and H is the Hubble parameter, which are respectively given by Here g ≈ 107 is the total number of effective relativistic degrees of freedom, and the Planck mass GeV. Finally, the DM relic density is given by where Y ∞ is the final DM abundance [50][51][52][53]. Observed relic density at 2σ is given by Ref. [54] as In the numerical analysis, we adopt a rather relaxed value range, 0.11 Ωh 2 0.13, and the LEP constraint expressed in Eq. (II.29) is imposed.
The Breit-Wigner effect causes a higher DM annihilation rate than for a thermally produced where σv rel c corresponds to x = v −2 ≈ 10 6 , while σv rel CMB corresponds to x ≈ 10 12 [55].

C. Numerical Analysis
We have performed a numerical analysis to find the allowed region for obtaining the correct relic density of DM, where neutrino oscillation data is implicitly reproduced and the LFV constraint We also conduct a numerical analysis on LFVs for these parameter sets, finding that if the Yukawa couplings are below ∼ 10 −2 , the experimental limits given in Eq. (II.23) are not exceeded for any of the parameter sets. Given that the Yukawa couplings are independent from the parameters defining the DM properties, there is no constraint from LFV on the studied parameter space.
With micrOMEGAs, we calculate the speed averaged cross section for v = 10 −3 and v = 10 −6 to obtain B and B CMB for each parameter set respectively, taking σv rel th = 0.12 Ωh 2 × 3 × 10 −26 cm 3 s −1 to compare with the generic model yielding the same relic density. Figure 1 shows boost  for annihilation to µ − + µ + . While velocity dependence of the annihilation cross section was not considered for these results, it can be assumed that the annihilation cross section in the CMB formation era is decisive. With the limits in principle being inversely proportional to the energy injected into the thermal bath, the limit for the annihilation of X can be calculated as For example, at m X = 400 GeV, B CMB > 40.32 would be excluded, with the part of the parameter space excluded by this and corresponding limits for other values of m X indicated in Figure 1.
We find B CMB ≈ B for most of the studied parameter space, except for two regions at γ 10 −7 as shown in Figure 2, matching the results shown in Ref. [58]. While there is a region in which B is up to three orders of magnitude larger than B CMB , it is ruled out by the constraint on B CMB . agree well for the e − + e + spectrum, which is a prerequisite for the combined fitting of the CALET e − + e + spectrum and the AMS-02 e + -only spectrum without assuming an inherent systematic offset. Due to the systematic difference of the DAMPE e − + e + spectrum results [22] from both AMS-02 and CALET spectra, we chose not to consider them in our study.

A. Electron and Positron Flux from Annihilation in the Galactic Halo
To predict the shape of the spectral component from annihilation of X andX , the positron spectra (identical to electron spectrum due to the symmetry of the process) per annihilation in the electron and muon channels have been calculated with PYTHIA 8.2 [59], which in turn were used as input for the propagation calculation with DRAGON [60] to obtain the flux at Earth. For the local DM density, ρ 0 = 0.3 GeV/cm 3 is assumed, and the speed averaged annihilation cross section normalized to the value predicted for a thermal relic DM, σv rel th = 3×10 −26 cm 3 s −1 . The choice of the DM halo shape model has no strong impact on the spectrum as the propagation range of electrons is limited and discussed models agree around the position of the solar system [61]. A NFW parametrization [62] is used for the calculation. The flux for annihilation of X andX is composed according to the branching ratio from Eq. (II.28) as the sum of one third of the normalized flux for electron channel and one third of the flux for muon channel, with the annihilation to neutrinos not contributing.
For the propagation calculation, we consider two propagation models strongly distinct in diffusion zone height L and diffusion coefficient normalization D 0 , denoted Model A and Model B.
Model A comprises a gradual change in the slope of the diffusion coefficient with rigidity [63] according to  represents an additional potential for positive charge only at low energy following Ref. [73].
positron excess such as secondary production in dense clouds around SNRs [74][75][76][77][78][79][80]. While a DMonly explanation of the positron excess is also not ruled out, it requires specific conditions such as decaying dark matter yielding softer spectra than the electron and muon channel annihilation of our DM candidate [81][82][83].
The model used for describing the background spectra and fitted to electron and positron cosmicray data is the sum of the above mentioned components, with the electron spectrum written as and the positron spectrum as Due to their large energy loss in propagation, the spectrum of primary electrons depends on the distribution of individual SNR in the galactic neighborhood of the solar system, which is yet mostly unknown. As an effective model of the local (after propagation) primary electron spectrum from all contributing SNR, it is parametrized by a power law with a soft spectral break (normalization C e , index γ e , break position E b and index change ∆γ e are free fit parameters, softness s = 0.05 is fixed) at low energy, and a high-energy exponential cut-off at E cut d representing radiative energy loss of high energy electrons. E cut d is not well constrained by the measurement and therefore treated as a fixed nuisance parameter for which we consider values of 2 TeV, 4 TeV and 10 TeV.
The secondary positron (Φ s(e + ) ) and electron (Φ s(e − ) ) fluxes are taken from the output of numerical propagation calculation with DRAGON for the nuclei spectra from which the propagation conditions are derived. The propagation conditions are used consistently for calculation of fluxes from secondaries, pulsars and Dark Matter. With an initial scale factor C norm obtained from normalizing the proton flux to measurements of AMS-02 [65], a common rescaling factor (C s /C norm ) is included in the fit as free parameter to account for remaining uncertainties in secondary particle production.
in which the characteristic diffusion distance r dif is expressed as The distance to the Monogem pulsar r = 0.28 kpc is taken from the ATNF catalog [84], as well as its age T = 1.11 × 10 2 kyr and energy loss rateĖ = 3.81 × 10 34 erg s −1 . The initial rotation energy of the pulsar Q 0 = 1.48 × 10 48 erg is calculated as Q 0 =Ė T 2 /τ , where τ = 10 kyr is the assumed spin-down timescale [70], so that the spectrum scales with the acceleration efficiency η, which is a free parameter in the fitting. A common assumption is that the accelerated particles are trapped for some time in the pulsar wind nebula (PWN) forming around the pulsar, and released with the dissolution of the PWN. The release delay T r is thus subtracted from the age T of the pulsar to determine the diffusion time t dif , with T r scanned in steps of 1 kyr considering the range up to 100 kyr [68].
The model is fitted to the data of CALET [21] based on total flux Φ − e + Φ + e and data of AMS-02 [19] for E > 10 GeV based on Φ + e by minimizing the sum of χ 2 of both comparisons, with systematic uncertainties of both measurements taken into account. For the CALET measurement, the 1σ deviation ∆ (k,i) as a function of each data point's energy (E i ) is listed in the supplemental material of Ref. [21] for the systematic uncertainty associated with the following parts of the analysis: Normalization, tracking, charge selection, electron identification, Monte Carlo model dependence. A systematic shift of the data-points is performed as part of the fit function with weights w k as free parameters and the squared weight of each uncertainty is added to the total χ 2 of the fit as given by where i iterates over the data points and k over the different systematic uncertainty types. Systematic errors associated with the trigger and the boosted decision tree proton rejection are added quadratically to the statistical error. For the AMS-02 measurement, the error on mean energy σ E in each bin is translated into an error on flux σ J(E) using the power law index γ e + also shown in Ref. [19] via the relation σ J(E) = J(E)(σ E /E)(γ e + − 1).
The lower boundary of E > 10 GeV for the data points used in the fitting is chosen due to charge and time dependent solar modulation effects expected below this energy [73]. Solar modulation effects above this energy are calculated using the force field approximation with a modulation potential of Φ = 500 MV, common for both charge signs. To check the potential influence of the parameter on our results, Φ = 300 MV and Φ = 700 MV are used as alternative fixed values.
The results of the base-model fit are shown in Figure 4 for the two propagation models. For Model A, the best fit is obtained with T r = 9 kyr, while for Model B, T r = 56 kyr yields lowest χ 2 .
The best fit for both propagation models uses E cut d = 10 TeV, which is thus taken as the default case. The reduced χ 2 is in either case χ 2 /ndof ≈ 0.5, indicating that the base model already more than adequately describes the data. With the Geminga pulsar as source of the positron excess, the fit quality is significantly worse unless T r > 100 kyr, which is the reason why we chose the Monogem pulsar.

D. Limit on Annihilation cross section from CALET and AMS-02 Data
The predicted flux from DM annihilation is added to the base model as an additional component of Φ ex with varied boost factor B, and the change of χ 2 studied. To derive a limit on B, or equivalently the speed averaged annihilation cross section σv rel c , B is increased in steps until χ 2 exceeds the 95% CL threshold for the respective number of degrees of freedom [85].
To determine the precise value of B for which the 95% CL threshold is crossed, the scan is repeated from the last allowed value with a factor 10 smaller step size, down to a step size of 0.01. To avoid reporting a too stringent limit due to the fitting function having no unique minimum, the "Migrad" and "Simplex" minimizers of Minuit [86] are used in alternation with different starting points as explained in Ref. [24]. Multiplying the normalization cross section by the scale factor at which the 95% CL threshold is crossed yields the limit on cross section σv rel c . By performing this procedure with m X scanned in steps of 5 GeV up to 500 GeV, 25 GeV from 500 to 1 TeV, 50 GeV from 1 TeV to 2 TeV and 100 GeV above 2 TeV, limits depending on m X are calculated, which are shown in Figure 5.
It is found that the limit varies only slightly under change of the nuisance parameters Φ and E cut d .
In principle, these limits are subject to the modeling of the astrophysical background flux being a good representation of the actual spectrum, disregarding whether it is an correct interpretation e.g. whether or not the Monogem pulsar is indeed the dominating source of the positron excess.
However to judge the conservativeness of the limits, it should be considered that for the peaked DM signal to be hidden by structures of the background from multiple astrophysical sources, these structures would have to form a deficit in a rather narrow energy range which can be considered an implausible coincidence given the smoothness of the spectrum in general. Due to these reasons, we consider the limits rather conservative, however to estimate the utmost possible influence of the background variability, also limits without any assumption on the background were calculated, using the method described above, but with only excess of the flux from DM annihilation over the flux measured by CALET and AMS-02 contributing to χ 2 . They are also shown in Figure 5 for comparison.

E. Structures in the CALET Spectrum as Possible Dark Matter Signatures
The addition of the predicted DM flux and increase of B for limit calculation causes eventually an increase in χ 2 of the fit. However, it is found that the addition of the DM flux with B smaller than the limit value improves the fit compared to the base model with a pulsar extra source in two ranges of m X , corresponding to step-like structures in the CALET spectrum. Given the excellent energy resolution combined with detailed energy calibration [87] over the wide dynamic range [88] of CALET, it is permissible to assume that the measured structures are features of the physical spectrum and not measurement artifacts, thus warranting an interpretation. To quantify the significance of interpreting the spectral structures as a signature of the proposed DM candidate, the optimal value of B and associated maximal χ 2 reduction are determined depending on m X .
FIG. 5. Limit on the annihilation cross section as a function of m X , compared to limits for e ± and µ ± channels from γ-ray observation of dwarf galaxies with Fermi-LAT from the supplemental material of Ref. [10] multiplied with three to account for the branching fraction. The purple area shows the variation of the limit among all cases with nuisance parameters Φ and E cut d changed between 300 MV,500 MV,700 MV and 2 TeV,4 TeV,10 TeV respectively. The dashed gray line indicates the limit without any background modeling, with only excess over measured flux contributing to χ 2 . Also, the shown limit is the worst in Using an approach similar to the limit calculation, B is initially scanned in 20 steps between zero and the 95% CL limit value, and then the interval around the value with best χ 2 scanned in nested intervals.
The best-fit B and χ 2 improvement as a function of DM mass are shown in Figure 6. The  Figure 1, with part of the 2σ-region being below the constraint from CMB anisotropy. It is also not ruled out by overproduction of γ-rays as shown by the comparison with limits from dwarf-galaxy observation by Fermi-LAT in Figure 5.
With both propagation models, the fit also improves with addition of the DM signal above 1 TeV where another a step-like structure exists in the CALET spectrum. For the default case, the fit improvement is maximal at m X = 1350 GeV for both Model A and Model B, but it is much less significant than the improvement for m X = 390 GeV. Apart from the larger errors, the significance is low, since even for the direct annihilation to e − + e + , the DM signal spectrum is not localized (hard) enough to match the structure well, as the best-fit graphs in Figure 8 demonstrate. The structure at m X = 390 GeV is thus a better candidate for being a DM signature, despite the 1 TeV structure being visually more prominent. After calculating the expected signature of the DM candidate in electron and positron cosmic rays for two largely distinct propagation models, we performed a combined search in the measured e − + e + CALET and e + AMS-02 spectra on top of an astrophysical background model assuming a single young pulsar as the source of the positron excess. As outcome we presented limits on the annihilation cross section close to those from γ-ray observation with Fermi-LAT, as well as a possible association of structures in the CALET spectrum with a DM signature. The significance of the fit improvement by adding the DM signature to the base model exceeds the 2σ-level depending on the propagation model, with the best fit for M DM at 390 GeV with a value of B which is well within the range predicted by the Breit-Wigner enhancement. These results demonstrate the significance of the step-like structure itself, and while other interpretations are possible, for example by overlapping spectra from individual astrophysical sources [89], it is shown that the annihilation of the DM candidate from the model presented herein also provides a suitable explanation.