Parallellization over perturbation modes does not scale optimally with many threads

[hersle]

SymBoltz parallellizes integration of independent perturbation modes over the number of threads that Julia is run with. Ideally, this should cut the runtime by almost the number of threads. But:

for n in 1 2 4 8 16; do
julia -t $n -E "using SymBoltz; using BenchmarkTools; M = ΛCDM(K = nothing); pars = SymBoltz.parameters_Planck18(M); prob = CosmologyProblem(M, pars); ks = range(1, 1000, length = 200); @btime sol = solve(prob, ks)"
done

  17.703 s (1953557 allocations: 321.27 MiB)
  11.434 s (1953569 allocations: 321.27 MiB)
  6.637 s (1953579 allocations: 321.27 MiB)
  3.483 s (1953599 allocations: 321.28 MiB)
  2.419 s (1953639 allocations: 321.28 MiB)

The scaling looks decent up to 8x, but could it be improved, particularly for more threads?

This was on a machine with 20 virtual cores:

julia> versioninfo()
Julia Version 1.11.7
Commit f2b3dbda30a (2025-09-08 12:10 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 20 × 12th Gen Intel(R) Core(TM) i7-12800H
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, alderlake)
Threads: 1 default, 0 interactive, 1 GC (on 20 virtual cores)

[hersle]

for n in 1 5 10 15 20; do
julia -t $n -E "using LinearAlgebra: BLAS; BLAS.set_num_threads(1); using SymBoltz; using BenchmarkTools; M = ΛCDM(K = nothing); pars = SymBoltz.parameters_Planck18(M); prob = CosmologyProblem(M, pars); ks = range(1, 1000, length = 200); @btime sol = solve(prob, ks)"
done

  11.413 s (1953557 allocations: 321.27 MiB)
  4.574 s (1953584 allocations: 321.27 MiB)
  3.553 s (1953609 allocations: 321.28 MiB)
  3.103 s (1953634 allocations: 321.28 MiB)
  3.669 s (1953659 allocations: 321.29 MiB)

This was with a charger connected, unlike the first example.