4 of tau lepton. When the SM particles X , Y , and Z are in the thermal equilibrium, ( n$_{Y}$ n$_{Z}$ ) / ( n eq Y n eq Z ) = ( n$_{X}$n$_{Y}$ ... ) / ( n eq $_{X}$n eq Y ... ) = 1, and hence tau leptons are sufficiently produced through the inverse annihilation and/or the inverse decay processes as long as these inter- action rates are larger than the Hubble expansion rate. Therefore, whether tau leptons are in the thermal bath or not can be distinguished by comparing the Hubble ex- pansion rate H with the inverse annihilation rate of tau lepton 〈 σv 〉 n eq $_{X}$, and the inverse decay rate of tau lepton 〈 Γ 〉 . In other words, the inequality expression indicates that tau leptons are in the thermal bath. Fig. 1 shows 〈 Γ 〉 n eq $_{X}$, 〈 Γ 〉 , and H as a function of the thermal bath temperature. As shown in Fig. 1, the inverse decay rate of tau lepton is much larger than the Hubble expan- sion rate. Thus, we can conclude that tau leptons remain in the thermal bath still at the beginning of the BBN. D. Calculation of the number density ratio of stau and neutralino 〈 σv 〉 n eq X > H and/or 〈 Γ 〉 > H (16) FIG. 1: The inverse annihilation rate of tau leptons 〈 Γ 〉 n $_{X}$, the inverse decay rate of tau leptons 〈 Γ 〉 , and the Hubble expansion rate H as a function of the thermal bath tempera- ture.  10  -40  10  -30  10  -20  10  -10  0.02  0.06  0.1  0.14  0.18 Interaction rate and expansion rate (GeV) Temperature (GeV) tau inverse annihilation rate tau inverse decay rate Hubble expansion rate We are now in a position to calculate the number den- sity ratio of stau and neutralino. In this subsection, we will show a set of relevant Boltzmann equations. The right-hand side of the Boltzmann equations (Eqs. (2), (3), and (4)) depends only on temperature, and hence it is convenient to use temperature T instead of time t as independent variable. To do this, we reformu- late the Boltzmann equations by using the ratio of the number density to the entropy density s : One is the competition between the interaction rate of the exchange processes and the Hubble expansion rate, since when these interaction rates get smaller than the Hubble expansion rate, the relic density of stau would be frozen out. The other is whether tau leptons are in the thermal bath or not. The interaction rate of the exchange processes strongly depends on the number density of tau leptons. When tau leptons are in the thermal bath, the number density ratio between stau and neutralino are given by the thermal ratio, through the exchange processes. On the contrary, once tau leptons decouple from the thermal bath, the ratio cannot reach this value. To calculate the relic density of stau, we have to comprehend the temperature of tau lepton decoupling. To see whether tau leptons are in the thermal bath or not, we consider the Boltzmann equation for its number density, n$_{τ}$ , n$_{˜}$ τ n$_{χ}$ ≃ e − m$_{˜}$ τ/T e − m$_{˜}$ χ/T = exp ( − δm T ) , (13) Y$_{i}$ = n$_{i}$ s . (17) Consequently, we obtain the Boltzmann equations for the number density evolution of stau and neutralino where indices X , Y , and Z denote the SM particles, and 〈 Γ 〉 represents the thermal averaged decay rate dn$_{τ}$ dt + 3 Hn$_{τ}$ = − 〈 σv 〉 [ n$_{τ}$n$_{X}$ − n eq τ n eq X ( n$_{Y}$ n$_{Z}$ n eq Y n eq Z )] − 〈 Γ 〉 n$_{τ}$ − n eq τ ( n$_{X}$n$_{Y}$ ... n eq $_{X}$n eq Y ... )] , (14) 〈 σv 〉 = ∑ $_{X,Y,Z}$〈 σv 〉 τ X ↔ Y Z , 〈 Γ 〉 = ∑ $_{X,Y,...}$〈 Γ 〉 τ ↔ $_{XY...}$, (15) dY$_{˜}$ τ − dT = [ 3 HT g$_{∗}$ ( T ) ] − 1 [ 3 g$_{∗}$ ( T ) + T dg$_{∗}$ ( T ) dT ] s × { 〈 σv 〉 ˜ τ − γ → ˜ χτ − Y$_{˜}$ τ − Y$_{γ}$ − 〈 σv 〉 ˜ χτ − → ˜ τ − $_{γ}$Y$_{˜}$ $_{χ}$Y$_{τ}$ − + 〈 σv 〉 ˜ τ − τ $^{+}$→ ˜ $_{χγ}$Y$_{˜}$ τ − Y$_{τ}$ + − 〈 σv 〉 ˜ χγ → ˜ τ − τ + Y$_{˜}$ $_{χ}$Y$_{γ}$ } (18) dY$_{˜}$ τ + dT = [ 3 HT g$_{∗}$ ( T ) ] − 1 [ 3 g$_{∗}$ ( T ) + T dg$_{∗}$ ( T ) dT ] s × { 〈 σv 〉 ˜ τ $^{+}$γ → ˜ χτ + Y$_{˜}$ τ + Y$_{γ}$ − 〈 σv 〉 ˜ χτ $^{+}$→ ˜ τ + $_{γ}$Y$_{˜}$ $_{χ}$Y$_{τ}$ + + 〈 σv 〉 ˜ τ $^{+}$τ − → ˜ $_{χγ}$Y$_{˜}$ τ + Y$_{τ}$ − − 〈 σv 〉 ˜ χγ → ˜ τ $^{+}$τ − Y$_{˜}$ $_{χ}$Y$_{γ}$ } (19)