3 This is because the cross sections of the exchange pro- cesses are in the same order of magnitude as that of the annihilation and the inverse annihilation, but the num- ber density of the SM particles is much larger than that of the supersymmetric particles which is suppressed by the Boltzmann factor. Thus even if the total number density of stau and neutralino is frozen out at the tem- perature T$_{f}$ , each number density of them continue to evolve through the exchange processes. Thus, to calculate the relic density of stau, we have to follow the two-step procedures. As a first step, we calculate the total relic density of the supersymmetric particles by solving the Eq. (7). We use the publicly available program micrOMEGAs [28] to calculate it. The second step is the calculation of the number density ratio of stau and neutralino. The second step is significant for calculating the relic density of stau at the BBN era, and hence we will discuss it in detail in the next subsection. problem by the long-lived stau, and the whole intention of this work is to search parameters which can provide the solution for the $^{7}$Li problem. Hence we assume that the stau is stable enough to survire until the BBN era, and focusing on the mass difference between stau and neutralino is small enough to make it possible. B. The evolution of the number density of stau and neutralino C. The exchange processes and Lagrangian for describing them For N , that is the number density of the neutralino, we get the Boltzmann equation by summing up Eqs. (2), (3), and (4), After the freeze-out of the total number density of stau and neutralino, each of them is exchanged through the following processes ˜ τγ ←→ ˜ χτ χγ ττ. (10) N = ∑ i n$_{i}$ . (6) dN dt + 3 HN = −〈 σv 〉 sum [ NN − N $^{eq}$Neq ] (7) 〈 σv 〉 sum ≡ ∑ i = ˜ χ, ˜ τ ∑ $_{X,Y}$〈 σv 〉 ˜ χi ↔ XY . (8) Notice that although there are other exchange processes via weak interaction (for example, ˜ τW ↔ ˜ χν$_{τ}$ , ˜ τν$_{τ}$ ↔ ˜ χW , and so on), we can omit them. This is because the number density of W boson is not enough to work these processes sufficiently due to the Boltzmann factor suppression, and the final state W boson is kinematically forbidden when the thermal bath temperature is less than T$_{f}$ . These processes (Eq. (10)) are described by the Lagrangian Notice that the terms describing the exchange processes in each Boltzmann equations cancel each other out. Solv- ing the Eq. (7), we obtain N and find the freeze out temperature of the total number density of all the super- symmetric particles T$_{f}$ by using the standard technique [26]: Now we will discuss the number density evolution of stau. To obtain the relic density of stau, we solve a cou- pled set of the Boltzmann equations, (2), (3), and (4), as simultaneous differential equation. Each Boltzmann equation contains the contributions of the exchange pro- cesses. These processes exchange stau with neutralino and vice versa. At the temperature T$_{f}$ , the interaction rate of the exchange processes is much larger than that of the annihilation and the inverse annihilation processes. Here, g and m$_{˜}$ χ are the internal degrees of freedom and the mass of neutralino, respectively. The Planck mass m$_{pl}$ = 1 . 22 × 10 19 GeV, and g$_{∗}$ are the total number of the relativistic degrees of freedom. Consequently, we see that 4 GeV ≲ T$_{f}$ ≲ 40 GeV for 100 GeV ≲ m$_{˜}$ χ ≲ 1000 GeV. where e is the electromagnetic coupling constant, P$_{L}$ and P$_{R}$ are the projection operators, and l ∈ { e, µ } . g$_{L}$ and g$_{R}$ are the coupling constants, given by where g is the SU (2)$_{L}$ gauge coupling constant, and θ$_{W}$ is the Weinberg angle. The evolution of the stau number density is governed only by the exchange processes (Eq.(10)) after the freeze- out of the total relic density of stau and neutralino. When we calculate it, we should pay attention to two essential points relevant to the exchange processes. g$_{L}$ = g √$_{2 cos}$$_{θ}$$_{W}$ sin θ$_{W}$ cos θ$_{τ}$ , g$_{R}$ = √$_{2}$$_{g}$ cos θ$_{W}$ sin θ$_{W}$ sin θ$_{τ}$ e $^{iγ$_{τ}$}$, (12) L = ˜ τ $^{∗}$˜ χ $^{0}$( g$_{L}$P$_{L}$ + g$_{R}$P$_{R}$ ) τ − ie (˜ τ $^{∗}$( ∂$_{µ}$ ˜ τ ) − ( ∂$_{µ}$ ˜ τ $^{∗}$)˜ τ ) A µ + h.c. , (11) m$_{˜}$ χ T$_{f}$ = ln 0 . 038 g m$_{pl}$ m$_{˜}$ χ 〈 σv 〉 g 1 / 2 ∗ ( m$_{˜}$ $_{χ}$/T$_{f}$ ) ≃ 25 . (9) In this subsection, we discuss the evolution of the num- ber density of each species. Firstly, we discuss the num- ber density evolution of neutralino DM. Since we have as- sumed R-parity conservation, all of the supersymmetric particles eventually decay into the LSP neutralino. Thus its final number density is simply described by the sum of the number density of all the supersymmetric particles :