2 dn$_{˜}$ τ + dt + 3 Hn$_{˜}$ τ + = − ∑ i ∑ $_{X,Y}$〈 σv 〉 ˜ τ $^{+}$i ↔ XY [ n$_{˜}$ τ + n$_{i}$ − n eq ˜ τ $_{+}$n eq i ( n$_{X}$n$_{Y}$ n eq $_{X}$n eq Y )] − i =˜ τ + $_{X,Y}${ 〈 σ $^{′}$v 〉 ˜ τ $^{+}$X → iY [ n$_{˜}$ τ + n$_{X}$ ] − 〈 σ $^{′}$v 〉 iY → ˜ τ $^{+}$X [ n$_{i}$n$_{Y}$ ] } (3) dn$_{˜}$ χ dt + 3 Hn$_{˜}$ χ = − ∑ i ∑ $_{X,Y}$〈 σv 〉 ˜ χi ↔ XY [ n$_{˜}$ $_{χ}$n$_{i}$ − n eq ˜ χ n eq i ( n$_{X}$n$_{Y}$ n eq $_{X}$n eq Y )] − i = ˜ χ $_{X,Y}${ 〈 σ $^{′}$v 〉 ˜ χX → iY [ n$_{˜}$ $_{χ}$n$_{X}$ ] − 〈 σ $^{′}$v 〉 iY → ˜ χX [ n$_{i}$n$_{Y}$ ] } . (4) II. FORMALISM FOR THE CALCULATION OF RELIC DENSITY OF STAU AT THE BBN ERA the yeild value of the DM and δm to be around 0 . 1 GeV, respectively. In this work, we improve our previous analysis by cal- culating Y$_{˜}$ τ, FO with taking the relic abundance of DM into account. The outline of this paper is as follows. In section II, we review the formalism for the calculation of the relic density of stau, and derive the Boltzmann equations of stau and neutralino for the calculation. In section III, we present the numerical results of the cal- culation, and prove the parameter space for solving the $^{7}$Li problem. Section IV is devoted to a summary and discussion. Here n and n eq represent the actual number density and the equilibrium number density of each particle, and H is the Hubble expansion rate. Index i denotes stau and neutralino, and indices X and Y denote SM particles. Note that if relevant SM particles are in thermal equi- librium, n$_{X}$ = n eq $_{X}$, n$_{Y}$ = n eq Y , and ( n$_{X}$n$_{Y}$ /n eq $_{X}$n eq Y ) = 1 then these equations are reduced into a familiar form. 〈 σv 〉 and 〈 σ $^{′}$v 〉 are the thermal averaged cross sections, which is defined by In this section, we prepare to calculate the relic density of stau at the BBN era. Firstly, in subsection II A, we briefly review the Boltzmann equations for the number density of stau and neutralino based on the thermal relic scenario. Then, in subsection II B, we discuss the number density evolution of stau and neutralino quantitatively. In subsection II C, we investigate the significant processes for the calculation of the relic density of stau. Finally, we obtain the Boltzmann equations for the relic density of stau in a convenient form. A. Boltzmann equations for the number density evolution of stau and neutralino 〈 σv 〉 12 → 34 ≡ g$_{12}$ ∫ d 3 p $_{1}$d 3 p 2 f$_{1}$f$_{2}$ ( σv )$_{12}$$_{→}$$_{34}$ ∫ d 3 p $_{1}$d 3 p 2 f$_{1}$f$_{2}$ = g$_{12}$ d 3 p $_{1}$d 3 p 2 f$_{1}$f$_{2}$ ( σv )$_{12}$$_{→}$$_{34}$ n eq 1 n eq 2 , (5) where f is the distribution function of a particle, v is the relative velocity between initial state particles, and g$_{12}$ = 2(1) for same (different) particles 1 and 2. In this work, we assume that all of the supersymmetric particles except for stau and neutralino are heavy, and therefore do not involve them in the coannihilation processes. } The first line on the right-hand side of Eqs. (2), (3), and (4) accounts for the annihilation and the inverse annihilation processes of the supersymmetric particles ( ij ↔ XY ). Here index j denotes stau and neutralino. As long as the R-parity is conserved, as shown later, the final number density of neutralino DM is controlled only by these processes. The second line accounts for the exchange processes by scattering off the cosmic ther- mal background ( iX ↔ jY ). These processes exchange stau with neutralino and vice versa, and thermalize them. Consequently, the number density ratio between them is controlled by these processes. Instead, these processes leave the total number density of the supersymmetric particles. Note that in general, although there are terms which account for decay and inverse decay processes of stau (˜ τ ↔ ˜ χXY... ) in the Boltzmann equations, we omit them. It is because we are interested in solving the $^{7}$Li In this subsection, we show the Boltzmann equations of stau and neutralino and briefly review their quantitative structure based on [26] (see also a recent paper[27]). We are interested in the relic density of stau in the coannihilation scenario. In this scenario, stau and neu- tralino are quasi-degenerate in mass and decouple from the thermal bath almost at the same time [2]. Thus the relic density of stau is given by solving a coupled set of the Boltzmann equations for stau and neutralino as si- multaneous differential equation. For simplicity, we use the Maxwell-Boltzmann statistics for all species instead of the Fermi-Dirac for fermions and the Bose-Einstein for bosons, and assume T invariance. With these simpli- fications, the Boltzmann equations of them are given as follows dn$_{˜}$ τ − dt + 3 Hn$_{˜}$ τ − = − ∑ i ∑ $_{X,Y}$〈 σv 〉 ˜ τ − i ↔ XY [ n$_{˜}$ τ − n$_{i}$ − n eq ˜ τ $_{−}$n eq i ( n$_{X}$n$_{Y}$ n eq $_{X}$n eq Y )] − i =˜ τ − $_{X,Y}${ 〈 σ $^{′}$v 〉 ˜ τ − X → iY [ n$_{˜}$ τ − n$_{X}$ ] − 〈 σ $^{′}$v 〉 iY → ˜ τ − X [ n$_{i}$n$_{Y}$ ] (2)