6 Eq. (13), n$_{˜}$ $_{τ}$/n$_{˜}$ χ ∼ exp( − δm/T ), since they follow the Boltzmann distribution before their freeze-out. Thus, the relic density of stau strongly depends on δm . Here, we comment on the dependence of the stau relic density n$_{˜}$ τ − on other parameters such as m$_{˜}$ $_{τ}$, θ$_{τ}$ , and γ$_{τ}$ . The number density of the negatively charged stau is expressed in terms of the total relic density N by n$_{˜}$ τ − = N 2(1 + e $^{δm/T$_{f}$$_{(ratio)}$}$) . (23) Here, the freeze-out temperature T$_{f}$$_{(ratio)}$ hardly depends on these pararameters. This is because the cross sec- tion of the exchange processes are changed by these pa- rameters at most by factors but not by orders of magni- tudes, and the T$_{f}$$_{(ratio)}$ depends logarithmically on 〈 σv 〉 as shown in Eq. (9). On the other hand, the total relic density N is proportional to m$_{˜}$ τ as in Eq. (22). The value of N is also affected by the left-right mixing θ$_{τ}$ as seen in Fig. 2 since the annihilation cross section de- pends on this parameter. In contrast, γ$_{τ}$ scarcely affects the relic density, since this parameter appears in the an- nihilation section through the cross terms of the contri- butions from the left-handed stau and the right-handed one, and such terms always accompany the suppression factor of m$_{τ}$/m$_{˜}$ τ compared to the leading contribution. Thus the relic number density of stau n$_{˜}$ τ strongly de- pends on m$_{˜}$ τ and θ$_{τ}$ while scarecely depends on γ$_{τ}$ . We comment on the generality of our method to calcu- late the density of exotic heavy particles that coannihi- late with other (quasi)stable particles: we calculate the total number density of these particles and then calcu- late the ratio among them by evaluating the exchange processes such as Eq. (10). This method of calculation can be found versatile in various scenarios including the catalyzed BBN and the exotic cosmological structure for- mation [31–33]. C. Long-lived stau and BBN B. Stau relic density at the BBN era FIG. 3: The evolution of the number density of negative charged stau. Each line attached [ δm ] shows the actual evo- lution of the number density of stau, while the one atattched [ δm (thermal)] shows its evolution under the equilibrium de- termined given by Eq. (13) and the total relic abundance. Yellow band represents the allowed region from the WMAP observation at the 2 σ level [1]. After the number density of stau freezes out, stau de- cays according to its lifetime [3], or forms a bound state with a nuclei in the BBN era. Their formation rate has been studied in literatures [14, 16, 17]. The bound states modify the predictions of SBBN, and make it possible to solve the $^{7}$Li problem via internal conversion processes in the bound state [14–16]. In Fig. 4 we show parameter regions that are consis- tent with the observed abundances of the DM and of the light elements. We calculate the relic density of stau by varying the value of δm with the values of m$_{˜}$ τ = 350 GeV, sin θ$_{τ}$ = 0.8, and γ$_{τ}$ = 0. With these parameters, the allowed region is shown inside the dotted oval. We see that there are allowed regions at Y$_{˜}$ τ − ∼ 10 − 13 – 10 − 12 for δm ≲ 130 MeV to solve the $^{7}$Li problem at 3 σ . On the other hand, it is found that the observational $^{6}$Li to $^{7}$Li ratio excludes Y$_{˜}$ τ − ≳ 10 − 15 and δm ≲ 100 MeV. We will explain this feature as follows. The number density evolution of stau is qualitatively understood as follows. As shown in Fig. 3, the freeze- out temperature of stau almost does not depend on δm . It is determined by the exchange processes Eq. (10), whose magnitude 〈 σv 〉 Y$_{˜}$ $_{τ}$Y$_{γ}$ is governed by the factor e − ( m$_{τ}$ − δm ) $^{/T}$, where m$_{τ}$ represents the tau lepton mass. The freeze-out temperature of the stau density T$_{f}$$_{(ratio)}$ is given by ( m$_{τ}$ − δm ) /T$_{f}$$_{(ratio)}$ ≃ 25 as in Eq. (9), since the cross section of the exchange process is of the same mag- nitude as weak processes. Thus T$_{f}$$_{(ratio)}$ hardly depends on δm . In contrast, the ratio of the number density be- tween stau and neutralino depends on δm according to Next, we solve the Boltzmann equations (18), (19), and (20) numerically, and obtain the ratio of the stau number density to the total number density of stau and neutralino. Fig. 3 shows the evolution of the number density of stau as a function of the universe tempera- ture. Here we took m$_{˜}$ τ = 350 GeV, sin θ$_{τ}$ = 0.8, and γ$_{τ}$ = 0 and chose δm = 10 MeV, 50 MeV, and 100 MeV as sample points. Each line attached [ δm ] shows the actual evolution of the number density of stau, while the one atattched [ δm (thermal)] shows its evolution under the equilibrium determined by Eq. (13) and the total relic abundance. Horizontal dotted line represents the relic density of DM, which is the total abundance calculated above. We took it as a initial condition of total value for the calulation of the number density ratio. Yellow band represents the allowed region from the WMAP observa- tion at the 2 σ level [1]. 10 MeV WMAP 10 MeV(thermal) 50 MeV 50 MeV(thermal) 100 MeV 100 MeV(thermal) DM abundance GLYPH<CC> GLYPH<CT>GLYPH<D1>GLYPH<D4> GLYPH<CT>GLYPH<D6>GLYPH<CP>GLYPH<D8>GLYPH<D9>GLYPH<D6>GLYPH<CT> GLYPH<B4>GLYPH<BZ>GLYPH<CT>GLYPH<CE>GLYPH<B5> GLYPH<CA>GLYPH<CT>GLYPH<D0>GLYPH<CX>  GLYPH<CS>GLYPH<CT>GLYPH<D2>GLYPH<D7>GLYPH<CX>GLYPH<D8> GLYPH<DD> GLYPH<D3>GLYPH<CU> GLYPH<D8>GLYPH<CW>GLYPH<CT> GLYPH<D7>GLYPH<D8>GLYPH<CP>GLYPH<D9>$^{GLYPH<CH>}$ GLYPH<DI> GLYPH<AS> GLYPH<BC>GLYPH<BA>GLYPH<BD> GLYPH<BC>GLYPH<BA>GLYPH<BC>GLYPH<BJ> GLYPH<BC>GLYPH<BA>GLYPH<BC>GLYPH<BG> GLYPH<BC>GLYPH<BA>GLYPH<BC>GLYPH<BD> GLYPH<BD>GLYPH<BC> GLYPH<A0> GLYPH<BD>GLYPH<BF> GLYPH<BD>GLYPH<BC> GLYPH<A0> GLYPH<BD>GLYPH<BE>