ABSTRACT (original article): We consider the scenario of self-interacting dark matter(SIDM) with a light mediator in a model-independent way, which can alleviate two long-standing issues of the small scale cosmology namely cusp vs. core and too-big-to-fail. A Yukawa potential is chosen to achieve mediator exchange between DM particles as part of their self-interactions. The dynamics of self-interacting transfer cross-section are studied for a range of mediator mass(

m
Z
). Also, a relationship is established between the cross-section and DM particles’ relative velocity, which ensures the solution to the DM crisis at small scales. Our obtained numerical results are efficient compared to the earlier works in the context that a lesser number of l modes have been used by us to achieve the same level of accuracy in the cross-section calculations. For a better understanding of the SIDM parameter space, we perform an analytical analysis on the dependence of transfer cross-section over the other important SIDM parameters using a Hulthén potential which is similar in its behavior to Yukawa potential. A detailed evolution of particle dynamics using the Boltzmann equation and the effect of Sommerfeld enhancement on such calculations has also been studied here. We also provide a minimal anomaly-free leptophilic extension of the standard model, that can incorporate SIDM and its mediator candidate in the framework. CITATION (original article): Utkarsh Patel, Sudhanwa Patra (2022), Numerically analyzing self-interacting dark matter, arXiv:2204.11551. https://doi.org/10.48550/arXiv.2204.11551
In this work, I present the numerical routines that I have used in my work arXiv: 2204.11551 for the model independent analysis of a Self-Interacting Dark Matter (SIDM) scenario. SIDM provides an elegant solution for the two long standing issues of the small scale cosmology namely cusp vs. core
1
and too-big-to-fail
2
. In the study, a Yukawa potential is chosen to achieve mediator exchange between DM particles as part of their self-interactions. The dynamics of self-interacting transfer cross-section is then studied for a range of mediator mass
M
Z'
. The code for obtaining SIDD cross-section in this notebook makes use of the instructions that are written in arXiv: 1302.3898. Based on these codes, one would be able to reproduce the plots and results that are provided in arXiv: 2204.11551.
Introduction

A bit background on the topic

Defining the Scattering Cross-section

◼
  • The required value of self interaction cross-section 0.10
    2
    cm
    g
    <
    σ
    T
    M
    DM
    < 10
    2
    cm
    g
    at dwarf scale to solve the core-cusp issue, is around 13-15 orders of magnitude larger than the typical WIMP scale cross-section. It is already known that standard freeze-out annihilation or production channels fail to provide such large cross-section values. Thus, a new gauge portal interaction involving a light mediator could serve the purpose. The key parameters on which the cross-section value of SIDM can depend on, are:
    • M
    DM
    ,
    M
    Z'
    ,
    α
    DM
    ,
    v
    DM
    . Here,
    α
    DM
    =
    2
    g
    DM
    4π
    ,
    v
    DM
    is the DM particle velocity, and
    g
    DM
    is the coupling strength of the new force mediator (Z’) with DM.
    ◼
  • The details of this cross-section are governed by non-perturbative calculations. Also, we are required to have a interaction cross-section whose strength die out as we move towards the larger scales of the universe, because non self-interacting WIMP dark matter simulation results matches excellently with the observations at large scales. The simplest motivation for such a cross-section can be drawn from the nucleonic interactions, where the interaction cross-section is velocity dependent. In such scenarios, Sommerfeld effect plays a vital role with the idea that the wave-function of the DM gets modified due to multiple exchange of light mediators before the main interaction, leading to large enhancement of the scattering cross-section. This mediator exchange ,as motivated from nucleonic interactions, is realised by a Yukawa-type potential, given as:
    • V(r)=±
    α
    DM
    r
    -
    M
    Z'
    r
    e
    (
    1
    )
    ◼
  • In the simple 2-particle scenario, the scattering of DM particles and their wave-function behaviour can be studied by solving the Schrodinger equation in the presence of external potential. This potential has already been chosen to be a Yukawa-type. The usual parameter to focus is the differential scattering cross-section
    dσ
    dΩ
    however if the mediators are light (<<1GeV), then the more relevant quantity is the transfer cross-section which is the usual differential cross-section weighted by a (1-Cos θ) parameter. Thus, transfer cross-section is given as:
    • σ
    T
    =
    dσ
    dΩ
    (1-Cosθ)dΩ
    (
    2
    )
    Here, θ is the scattering angle .

    Various Regimes

    Researchers have developed analytical solutions for the transfer cross-section under specific conditions and assumptions. The overall picture can be categorized in three different regimes. Here, we discuss those regimes.

    Born Regime

    ◼
  • The transfer cross-section
    σ
    T
    when computed perturbatively (i.e. in the Born approximation) keeping only the dominant t-channel contributions in the interactions, is of the form
    3
    :
    • σ
    T
    =
    2π
    2
    M
    Z'
    2
    β
    ln(1+
    2
    R
    )-
    2
    R
    1+
    2
    R
    
    (
    3
    )
    Here, β =
    2
    α
    DM
    *
    M
    Z'
    M
    DM
    *
    2
    v
    DM
    and R ≡
    M
    DM
    *
    v
    DM
    M
    Z'
    . This above equation (3) works well only within the limit
    α
    DM
    *
    M
    DM
    M
    Z'
    <<1
    and is called as the Born regime.

    Classical Regime

    ◼
  • The parameter space where Born approximation loses validity is characterized by non-perturbative corrections. In here, a considerable amount of interactions occur outside the Debye sphere and thus they provide a substantial amount of contribution to the transfer cross-section. Within this regime, analytical formula has been derived in the literature only for the cases where quantum mechanical effects are sub-dominant
    4
    . This is marked by the condition,
    M
    DM
    *
    v
    DM
    M
    Z'
    >>1
    and is called as the Classical Regime. For an attractive Yukawa-type potential, the expressions for
    σ
    T
    for various β ranges as given in reference
    5
    are:​
    • σ
    T
    =
    8π
    2
    M
    Z'
    2
    β
    1+1.5
    1.65
    β
    
    -1
    10
    ⩽β⩽
    3
    10
    π
    2
    M
    Z'
    2
    lnβ+1-
    1
    2
    -1
    ln
    β
    β⩾
    3
    10
    (
    4
    )

    Resonant Regime

    ◼
  • Upon analyzing the expressions in Eqs. (3) & (4), it can be seen that a vast amount of parameter space remains unexplored. This region can be expressed as being governed by both non-perturbative and quantum mechanical effects. The authors in reference
    6
    have focused broadly on this region, which they call as the “resonant regime” indicating the occurrence of quantum mechanical resonances in the value of
    σ
    T
    due to the Sommerfeld effect from the potential. There is no analytical expression available in the literature for
    σ
    T
    in this parameter space region. Thus, the only possibility is to numerically solve the Schrodinger equation to obtain the transfer cross-section. The techniques for numerically solving this two particle Schrodinger equation have been given in references
    7-8
    . The codes here have been developed following the steps in those works. In our analysis
    9
    , a slight change in the methodology has been implemented where the obtained values of
    x
    i
    and
    x
    m
    for s-wave mode is used as initial conditions for l>0 cases. This makes the obtained value of
    σ
    T
    saturate faster i.e. the desired
    σ
    T
    is obtained with lesser l modes.
    • Numerical Implementation
    In this chapter, I present the numerical methods for evaluating the scattering cross-section (
    σ
    T
    ). I will first layout the common variables that are required as a part of initializations for the data structures used ahead along with defining common parameter values. Then, I will present the analysis for the three regimes separately in three sections. Wherever required or possible, I will try to provide the supporting text for the code. As I am still learning Mathematica, so the code is not optimised well so one might find better and more efficient methods to solve the same problem.

    Common Variables:

    In[]:=
    SetDirectory@NotebookDirectory[];​​Off[InterpolatingFunction::dmval]
    In[]:=
    Array[δl,200,0];​​tab0={};​​tab1={};​​tab2={};​​tab3={};
    In[]:=
    v=N@
    1000*10
    3*
    8
    10
    ​​mx=200;​​αx=
    -2
    10
    ;
    Out[]=
    0.0000333333

    Born Regime:

    The code below stores the value of scaled
    σ
    T
    to the list tab3 along with the corresponding mphi value based on the formula given in Eq. (3). The variable mphi used hear is the mass of mediator particle (
    M
    Z'
    in our case).
    In[]:=
    TableAppendTotab2,mphi,
    8π*αx*αx
    mx*mx*v*v*v*v
    *Log1+
    mx*mx*v*v
    mphi*mphi
    -
    (mx*mx*v*v)
    (mphi*mphi)+(mx*mx*v*v)
    *4*
    -28
    10
    (mx*1.8*
    -24
    10
    ),{mphi,NestWhileList[If[#==10*10^Floor[Log10[Abs[#]]],#+10^Floor[Log10[Abs[1*#]]],#+10^Floor[Log10[Abs[0.1*#]]]]&,0.001,#<=8.6&,1,∞,-1]}​​
    Out[]=
    {{{0.001,160513.}},
    ⋯345⋯
    ,{{0.001,160513.},{0.0011,150249.},{0.0012,140968.},{0.0013,132520.},{0.0014,124785.},{0.0015,117671.},{0.0016,111101.},{0.0017,105013.},{0.0018,99354.8},{0.0019,94082.8},{0.002,89159.8},
    ⋯326⋯
    ,{7.7,1.58844×
    -8
    10
    },{7.8,1.50862×
    -8
    10
    },{7.9,1.43334×
    -8
    10
    },{8.,1.36325×
    -8
    10
    },{8.1,1.29776×
    -8
    10
    },{8.2,1.23556×
    -8
    10
    },{8.3,1.17713×
    -8
    10
    },{8.4,1.12158×
    -8
    10
    },{8.5,1.0704×
    -8
    10
    },{8.6,1.02126×
    -8
    10
    }}}
    Full expression not available
    (
    original memory size:
    5.8 MB)

    Classical Regime:

    In[]:=
    beta=
    2αx*mphi
    mx*v*v
    *
    0
    10
    Out[]=
    90000.mphi
    The code below checks the condition on β parameter given in Eq. (4), and accordingly append the value of scaled
    σ
    T
    to the list tab3 along with the corresponding mphi value. The NestWhileList function iterates over values of mphi from 0.01 to 100, incrementing mphi by powers of 10 based on its magnitude.
    In[]:=
    TableIfbeta<=
    3
    10
    ,AppendTotab3,mphi,
    8π
    mphi*mphi
    *
    beta*beta
    1+(1.5*
    1.65
    beta
    )
    *4*
    -28
    10
    (mx*1.8*
    -24
    10
    ),AppendTotab3,mphi,
    π
    mphi*mphi
    *Log[beta]+1-
    1
    2
    -1
    (Log[beta])
    *Log[beta]+1-
    1
    2
    -1
    (Log[beta])
    *4*
    -28
    10
    (mx*1.8*
    -24
    10
    ),{mphi,NestWhileList[If[#==10*10^Floor[Log10[Abs[#]]],#+10^Floor[Log10[Abs[1*#]]],#+10^Floor[Log10[Abs[0.1*#]]]]&,0.001,#<=8.6&,1,∞,-1]}
    Out[]=
    {{{0.001,89.8914}},{{0.001,89.8914},{
    ⋯1⋯
    }},
    ⋯344⋯
    ,{{0.001,89.8914},{0.0011,76.8149},{0.0012,66.5448},{0.0013,58.314},{0.0014,51.6037},{0.0015,46.0524},{0.0016,41.4013},{0.0017,37.4609},{0.0018,34.0899},{0.0019,31.1808},
    ⋯327⋯
    ,{7.7,0.0000122279},{7.8,0.0000119377},{7.9,0.0000116581},{8.,0.0000113883},{8.1,0.000011128},{8.2,0.0000108767},{8.3,0.0000106341},{8.4,0.0000103996},{8.5,0.000010173},{8.6,9.95383×
    -6
    10
    }}}
    Full expression not available
    (
    original memory size:
    5.8 MB)

    Numerical Regime:

    Attractive potential case (V = -
    α
    DM
    r
    -
    M
    Z'
    r
    e
    )

    In[]:=
    a=
    v
    2*αx
    ;​​b=
    αx*mx
    mphi
    ;​​mu=
    mx
    2
    ;​​k=mu*v;
    In[]:=
    Tableflag=0;sigmaT=0;​​l=0;​​Table​​xi=Quiet[xi0/.NSolve[100*xi0==b,xi0][[1]]];s1=NDSolvey''[x]+
    2
    a
    -
    l(l+1)
    2
    x
    +
    1
    x
    Exp
    -x
    b
    y[x]0,y[xi]==1,y'[xi]
    (l+1)
    xi
    ,y,{x,xi,xm2+1},MaxStepSizexi*0.1;​​s2=NDSolvey''[x]+
    2
    a
    -
    l(l+1)
    2
    x
    +
    1
    x
    Exp
    -x
    b
    y[x]0,y[xi]==1,y'[xi]
    (l+1)
    xi
    ,y',{x,xi,xm2+1},MaxStepSizexi*0.1;​​t1=Table[{x,First[y[x]/.s1]},{x,xi,xm2+1,xi}];​​t2=Table[{x,First[y'[x]/.s2]},{x,xi,xm2+1,xi}];​​f1=Quiet@Interpolation[t1];​​f2=Quiet@Interpolation[t2];​​bl=(xm2*f2[xm2])f1[xm2]-1;​​δl[l]=ArcTan[((a*xm2*D[SphericalBesselJ[l,x2],x2]/.x2a*xm2)-bl*SphericalBesselJ[l,a*xm2])/((a*xm2*D[SphericalBesselY[l,x2],x2]/.x2a*xm2)-bl*SphericalBesselY[l,a*xm2])];​​(*Print["
    x
    max
    =",xm2,"\t","
    x
    i
    =",xi,"\t","l=",l,"\t","δl[l]=",δl[l]];​​Print["Percentage Change=",((δl[l]-flag)/δl[l])*100];​​*)If[(Abs[(δl[l]-flag)/δl[l]])*100<1,Return[Null,Table]];​​flag=δl[l];​​count=xm2,xm2,Quietxm/.NSolvea*a==10*
    Exp
    -xm
    b
    
    xm
    ,xm[[1]],20000,1;​​flag=0;​​xm2=count;​​Table​​s1=NDSolvey''[x]+
    2
    a
    -
    l(l+1)
    2
    x
    +
    1
    x
    Exp
    -x
    b
    y[x]0,y[xi]==1,y'[xi]
    (l+1)
    xi
    ,y,{x,xi,xm2+1},MaxStepSizexi*0.1;​​s2=NDSolvey''[x]+
    2
    a
    -
    l(l+1)
    2
    x
    +
    1
    x
    Exp
    -x
    b
    y[x]0,y[xi]==1,y'[xi]
    (l+1)
    xi
    ,y',{x,xi,xm2+1},MaxStepSizexi*0.1;​​t1=Table[{x,First[y[x]/.s1]},{x,xi,xm2+1,xi}];​​t2=Table[{x,First[y'[x]/.s2]},{x,xi,xm2+1,xi}];​​f1=Quiet@Interpolation[t1];​​f2=Quiet@Interpolation[t2];​​bl=(xm2*f2[xm2])f1[xm2]-1;​​δl[l]=ArcTan[((a*xm2*D[SphericalBesselJ[l,x2],x2]/.x2a*xm2)-bl*SphericalBesselJ[l,a*xm2])/((a*xm2*D[SphericalBesselY[l,x2],x2]/.x2a*xm2)-bl*SphericalBesselY[l,a*xm2])];​​(*Print["
    x
    max
    =",xm2,"\t","
    x
    i
    =",xi,"\t","l=",l,"\t","δl[l]=",δl[l]];​​Print["Percentage Change=",((δl[l]-flag)/δl[l])*100];*)​​sigmaT=sigmaT+
    4π
    2
    k
    (l)*(Sin[δl[l]-δl[l-1]])*(Sin[δl[l]-δl[l-1]]);​​(*If[(Abs[(sigmaT-flag)/sigmaT])*100<1,Return[Null,Table]];*)​​flag=sigmaT;​​,{l,1,40,1};​​(*Print"Conversion= 1
    -2
    GeV
    =4
    -28
    x10
    2
    cm
    ","\t","1GeV=1.8
    -24
    x10
    g";*)​​(*Print"
    m
    ϕ
    =",mphi,"\t","Sigma(​
    -2
    GeV
    ​)=",sigmaT,"\t","Sigma(​
    2
    cm
    ​)=",sigmaT*4*
    -28
    10
    ,"\t","mχ(GeV)=",mx,"\t","mχ(grams)=",(mx*1.8*
    -24
    10
    ),"\t","
    Sigma
    mχ
    ​​
    2
    cm
    /g=",sigmaT*4*
    -28
    10
    (mx*1.8*
    -24
    10
    ),"\n\n";*)​​AppendTotab0,mphi,sigmaT*4*
    -28
    10
    (mx*1.8*
    -24
    10
    )​​(*WriteStringfile,mphi,"\t\t",NumberFormsigmaT*4*
    -28
    10
    (mx*1.8*
    -24
    10
    ),ExponentFunction(Null&),"\n"*)​​,{mphi,0.001,0.01,0.0001}

    Plotting

    ◼
  • Please feel free to contact me over utkarshp@iitbhilai.ac.in for any issues/queries or discussions surrounding this notebook’s work. Also, please cite this work if you find these codes useful for your study.

    • CITE THIS NOTEBOOK

    Scattering cross-section analysis for self-interacting dark matter​
    by Utkarsh Patel​
    Wolfram Community, STAFF PICKS, April 22, 2024
    ​https://community.wolfram.com/groups/-/m/t/3162442