Thank you again. I replaced sum and diff tanh functions with unitstep function but still solutions didn't satisfy nonlinear equations.

In[1]:= d = 1;
T = 0.0001 d;
\[Epsilon]f0 = - 0.102001 d;
p = 0.875;
q0 = 0.5;
v = 0.5 d;
ds[\[Epsilon]_] := 3/(4 Sqrt[2] d) Sqrt[(\[Epsilon] + d)/d];

(*we calculate renormalized position of the f-level, chemical \
potential and slave boson amplitude by solving three nonlinear \
equations*)

\[Zeta][\[Epsilon]_, \[Mu]_] := \[Epsilon] - \[Mu];  \
(*\[Epsilon]=k^2/(2 m)-d*)
Rk[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
 Sqrt[(\[Epsilon]f - \[Zeta][\[Epsilon], \[Mu]])^2 + (2 v a)^2];
Enkp[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  1/2 (\[Epsilon]f + \[Zeta][\[Epsilon], \[Mu]] + 
     Rk[\[Epsilon], \[Mu], \[Epsilon]f, a]);
Enkm[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  1/2 (\[Epsilon]f + \[Zeta][\[Epsilon], \[Mu]] - 
     Rk[\[Epsilon], \[Mu], \[Epsilon]f, a]);

In[11]:= fwp[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  1/2 (1 - Tanh[Enkp[\[Epsilon], \[Mu], \[Epsilon]f, a]/(2 T)]);
fwm[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  1/2 (1 - Tanh[Enkm[\[Epsilon], \[Mu], \[Epsilon]f, a]/(2 T)]);
diff[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  1 - UnitStep[\[Epsilon] - (a^2 + 
       4 \[Epsilon]f \[Mu])/(4 \[Epsilon]f)];
sum[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  1 - UnitStep[\[Epsilon] - (a^2 + 
       4 \[Epsilon]f \[Mu])/(4 \[Epsilon]f)];
eq1int[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  diff[\[Epsilon], \[Mu], \[Epsilon]f, a]/
  Rk[\[Epsilon], \[Mu], \[Epsilon]f, a];
eq2int[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, 
   a_] := (\[Zeta][\[Epsilon], \[Mu]] - \[Epsilon]f) \
eq1int[\[Epsilon], \[Mu], \[Epsilon]f, a];

eq3int[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  sum[\[Epsilon], \[Mu], \[Epsilon]f, a];

fint1[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  ds[\[Epsilon]] eq1int[\[Epsilon], \[Mu], \[Epsilon]f, a];
fint2[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  ds[\[Epsilon]] eq2int[\[Epsilon], \[Mu], \[Epsilon]f, a];
fint3[\[Epsilon]_, \[Mu]_, \[Epsilon]f_, a_] := 
  ds[\[Epsilon]] eq3int[\[Epsilon], \[Mu], \[Epsilon]f, a];

In[21]:= eqn1[\[Mu]_?NumericQ, \[Epsilon]f_?NumericQ, 
  a_?NumericQ] := \[Epsilon]f - \[Epsilon]f0 - 
  v^2 NIntegrate[
    fint1[\[Epsilon], \[Mu], \[Epsilon]f, a], {\[Epsilon], -d, d}, 
    MaxRecursion -> 20, 
    Method -> {GlobalAdaptive, MaxErrorIncreases -> 8000}]
eqn2[\[Mu]_?NumericQ, \[Epsilon]f_?NumericQ, a_?NumericQ] := 
 2 (q0 - a^2 ) - p - 
  NIntegrate[
   fint2[\[Epsilon], \[Mu], \[Epsilon]f, a], {\[Epsilon], -d, d}, 
   MaxRecursion -> 20, 
   Method -> {GlobalAdaptive, MaxErrorIncreases -> 8000}]
eqn3[\[Mu]_?NumericQ, \[Epsilon]f_?NumericQ, a_?NumericQ] := 
 p - NIntegrate[
   fint3[\[Epsilon], \[Mu], \[Epsilon]f, a], {\[Epsilon], -d, d}, 
   MaxRecursion -> 20, 
   Method -> {GlobalAdaptive, MaxErrorIncreases -> 8000}]

In[24]:= mysol = 
 FindRoot[{eqn1[\[Mu], \[Epsilon]f, a], eqn2[\[Mu], \[Epsilon]f, a], 
   eqn3[\[Mu], \[Epsilon]f, a]}, {{\[Mu], -0.42628}, {\[Epsilon]f, 
    0.035940}, {a, -0.37342}}, WorkingPrecision -> 16, 
  AccuracyGoal -> 15, MaxIterations -> Infinity]


During evaluation of In[24]:= FindRoot::jsing: Encountered a singular Jacobian at the point {\[Mu],\[Epsilon]f,a} = {0.845381,0.196724,-0.730453}. Try perturbing the initial point(s).

Out[24]= {\[Mu] -> 0.845381, \[Epsilon]f -> 0.196724, a -> -0.730453}

In[25]:= {eqn1[\[Mu], \[Epsilon]f, a], eqn2[\[Mu], \[Epsilon]f, a], 
  eqn3[\[Mu], \[Epsilon]f, a]} /. mysol

Out[25]= {0.0630848, -0.292974, -0.125}

I would say

2 UnitStep[\[Epsilon] - (a^2 + 
       4 \[Epsilon]f \[Mu])/(4 \[Epsilon]f)] - 2

To check it:

Manipulate[
 Plot[Tanh[
    2500 (\[Epsilon] + \[Epsilon]f - \[Mu] - 
       Sqrt[a^2 + (-\[Epsilon] + \[Epsilon]f + \[Mu])^2])] - 
   Tanh[2500 (\[Epsilon] + \[Epsilon]f - \[Mu] + 
       Sqrt[a^2 + (-\[Epsilon] + \[Epsilon]f + \[Mu])^2])] - (
     2 UnitStep[\[Epsilon] - (a^2 + 
           4 \[Epsilon]f \[Mu])/(4 \[Epsilon]f)] - 
     2), {\[Epsilon], -1, 1}],
 {{\[Mu], -.2152}, -1, 1},
 {{\[Epsilon]f, 1.2474}, 0, 2},
 {{a, .1}, -1, 1}]

Thank you for the reply.

The one you gave me was for sum of Tanh functions.

Could you please give me an expresssion for difference of Tanh functions as well?

Tanh[2500 (\[Epsilon] + \[Epsilon]f - \[Mu] - Sqrt[
     a^2 + (-\[Epsilon] + \[Epsilon]f + \[Mu])^2])]- 
 Tanh[2500 (\[Epsilon] + \[Epsilon]f - \[Mu] + Sqrt[
     a^2 + (-\[Epsilon] + \[Epsilon]f + \[Mu])^2])]