Computing eigenvalues and eigenvectors of large matrix and scaling of data

Question

I am trying to plot graphs 1(a1,b1,c1) from https://arxiv.org/pdf/2404.14971 by solving Eq.1. This gives a matrix. For the energy plot, we have to find the lowest 2 eigenvalues and find the difference between them for a range of h and some different matrix lengths. For localization length (Eq.2) and IPR (Eq.7) we have to follow the same process but in this case, we need an eigenvector corresponding to the lowest eigenvalue. Here is my code

H[deltaj_, \[Delta]value_, \[Phi]value_, deltah_, length_Integer] := 
  SparseArray[{{i_, 
      i_} -> (2*deltaj + \[Delta]value)*
       Cos[2 \[Pi] i Fibonacci[length - 1]/Fibonacci[length] + 
         2 \[Pi] \[Phi]value] + i deltah, {i_, j_} /; 
      Abs[i - j] == 1 -> -1}, {length, length}];(*The matrix from 1*)
Clear[index]; (*need for localization length*)
index[length_] := Range[1, length];
With[{deltaj = 1 (*deltaj=j, deltah=h in paper*)},
  h[\[Delta]value_, \[Phi]value_, deltah_, length_] = 
   H[deltaj, \[Delta]value, \[Phi]value, deltah, length];
  Clear[gs];
  gs[\[Delta]value_, \[Phi]value_?NumericQ, deltah_?NumericQ, length_,
     m_Integer /; m >= 1] := 
   gs[\[Delta]value, \[Phi]value, deltah, length, 
     m] = -Eigensystem[-h[\[Delta]value, N[\[Phi]value], N[deltah], 
           length], m, 
         Method -> {"Arnoldi", "Criteria" -> "RealPart", 
           MaxIterations -> 10^6}] // Transpose // Sort // Transpose;
  ];

(*calculating energy difference*)
Do[ \[Phi]value = RandomReal[WorkingPrecision -> 8]; 
 deltaE = Table[{deltah, 
    gs[0, \[Phi]value, deltah, length, 2][[1, 2]] - 
     gs[0, \[Phi]value, deltah, length, 2][[1, 1]]}, {deltah, 10^-8, 
    10^-4, 9*10^-8}, {length, {144, 233, 377, 610, 987}}], 1]

(*Plotting energy difference with h and system size*)
ListLogLogPlot[{deltaE[[All, 1]], deltaE[[All, 2]], deltaE[[All, 3]], 
  deltaE[[All, 4]], deltaE[[All, 5]]}, 
 PlotRange -> {{10^-8, 10^-4}, {10^-6, 10^-2}}, 
 PlotLegends -> {"L=144", "L=233", "L=377", "L=610", "L=987"}, 
 ImageSize -> 500, 
 PlotStyle -> {{Blue, Thickness[0.5]}, {Orange, 
    Thickness[0.5]}, {Green, Thickness[0.5]}, {Red, 
    Thickness[0.5]}, {Purple, Thickness[0.5]}}, 
 LabelStyle -> {24, Bold, Large, Black}, Frame -> True, 
 FrameStyle -> Thickness[0.003], 
 FrameTicksStyle -> {Black, Thickness[0.003]}, 
 FrameLabel -> {{Style["\[CapitalDelta]E", 
     FontFamily -> "Times New Roman", FontSlant -> "Italic", 
     FontWeight -> Bold, FontSize -> 30], 
    None}, {Style["h", FontFamily -> "Times New Roman", 
     FontSlant -> "Italic", FontWeight -> Bold, FontSize -> 30], 
    None}}]

(*calculating localized length*)
Do[ \[Phi]value = RandomReal[WorkingPrecision -> 8]; 
 locallength = 
  Table[{deltah, Sqrt[
    Total[(index[length] - 
        Total[index[length]*
          Abs[gs[0, \[Phi]value, deltah, length, 1][[2, 
            1]]]^2])^2 Abs[
       gs[0, \[Phi]value, deltah, length, 1][[2, 1]]]^2]]}, {deltah, 
    10^-8, 10^-4, 9*10^-8}, {length, {144, 233, 377, 610, 987}}], 1]

(*plotting localized length*)
ListLogLogPlot[{locallength[[All, 1]], locallength[[All, 2]], 
  locallength[[All, 3]], locallength[[All, 4]], 
  locallength[[All, 5]]}, PlotRange -> {{10^-8, 10^-4}, {5, 150}}, 
 PlotLegends -> {"L=144", "L=233", "L=377", "L=610", "L=987"}, 
 ImageSize -> 500, 
 PlotStyle -> {{Blue, Thickness[0.5]}, {Orange, 
    Thickness[0.5]}, {Green, Thickness[0.5]}, {Red, 
    Thickness[0.5]}, {Purple, Thickness[0.5]}}, 
 LabelStyle -> {24, Bold, Large, Black}, Frame -> True, 
 FrameStyle -> Thickness[0.003], 
 FrameTicksStyle -> {Black, Thickness[0.003]}, 
 FrameLabel -> {{Style["\[Xi]", FontFamily -> "Times New Roman", 
     FontSlant -> "Italic", FontWeight -> Bold, FontSize -> 30], 
    None}, {Style["h", FontFamily -> "Times New Roman", 
     FontSlant -> "Italic", FontWeight -> Bold, FontSize -> 30], 
    None}}]

(*calculating IPR*)
Do[ \[Phi]value = RandomReal[WorkingPrecision -> 8]; 
 ipr = Table[{deltah, 
    Total[Abs[
      gs[0, \[Phi]value, deltah, length, 1][[2, 1]]]^4]}, {deltah, 
    10^-8, 10^-4, 9*10^-8}, {length, {144, 233, 377, 610, 987}}], 1]
(*plotting IPR*)
ListLogLogPlot[{ipr[[All, 1]], ipr[[All, 2]], ipr[[All, 3]], 
  ipr[[All, 4]], ipr[[All, 5]]}, 
 PlotRange -> {{10^-8, 10^-4}, {0.05, 0.3}}, 
 PlotLegends -> {"L=144", "L=233", "L=377", "L=610", "L=987"}, 
 ImageSize -> 500, 
 PlotStyle -> {{Blue, Thickness[0.5]}, {Orange, 
    Thickness[0.5]}, {Green, Thickness[0.5]}, {Red, 
    Thickness[0.5]}, {Purple, Thickness[0.5]}}, 
 LabelStyle -> {24, Bold, Large, Black}, Frame -> True, 
 FrameStyle -> Thickness[0.003], 
 FrameTicksStyle -> {Black, Thickness[0.003]}, 
 FrameLabel -> {{Style["IPR", FontFamily -> "Times New Roman", 
     FontSlant -> "Italic", FontWeight -> Bold, FontSize -> 30], 
    None}, {Style["h", FontFamily -> "Times New Roman", 
     FontSlant -> "Italic", FontWeight -> Bold, FontSize -> 30], 
    None}}]

I need help with a couple of cases;

1. This code takes a very long time 20-30min to generate each graph if I take a small h like 10^(-8). If anyone can help me make it faster.
2. How can I scale them like fig. 1(a2,b2,c2) and find critical exponents like z,s v.
3. Originally I had to run the DO loop many times (maybe 100 times) (instead of a single time which I have done in the attached code ) for different random phi values and take an average of those data and then plot. How can I take an average of those data with every loop?
4. How can I run this code in the cluster? Does the code need any modification or a special name for the Mathematica file?