Solve for the correspondence between $\gamma$ and $l$ for the case of $h$ maximum

Question

The analytic expression for $h$ is the following function on $l$, $p$, $\gamma$: $$ h=2^{\left| l \right|+2}p!\left( \left| l \right|+p \right) !\gamma ^2\left( \sum_{m=0}^p{\frac{\left( -2 \right) ^m\varGamma \left( m+\left| l \right|/2+1 \right)}{\left( p-m \right) !\left( \left| l \right|+m \right) !m!\left( 1+\gamma ^2 \right) ^{m+\left| l \right|/2+1}}} \right) ^2 $$

h[l_, p_, \[Gamma]_] := 
  2^(Abs[l] + 2) p! (Abs[l] + 
      p)! \[Gamma]^2 (Sum[((-2)^m Gamma[
        m + Abs[l]/2 + 
         1])/((p - m)! (Abs[l] + m)! m! (1 + \[Gamma]^2)^(m + 
          Abs[l]/2 + 1)), {m, 0, p}])^2;

When $p=0$:

We can find the relationship between $\gamma$ and $l$ for the condition that $h$ is at its maximum value.

$Assumptions = \[Gamma] > 0;
Reduce[D[h[l, 0, \[Gamma]], \[Gamma]] == 0, \[Gamma]] // FullSimplify

[Gamma] == 1/Sqrt[1 + Abs[l]]

i.e. $\gamma =\frac{1}{\sqrt{1+\left| l \right|}}$

But as $p$ increases, the expression for $h$ becomes complicated, and perhaps when $h$ is maximum, it is impossible to find the correspondence between $\gamma$ and $l$. Although I can calculate the result by numerical solution, I ideally can solve the correspondence between $\gamma$ and $l$, for example, when $p=1,2,3...$.

For example: when $l=1$, $p=2$.

Plot[h[1, 2, \[Gamma]], {\[Gamma], 0, 10}, PlotRange -> All]
NMaximize[h[1, 2, \[Gamma]], {\[Gamma], 0, 10}]

{0.162661, {[Gamma] -> 0.311767}}

Answer 1

Function

maxi[p_Integer] :=Block[
{sol = NMaximize[h[l, p , \[Gamma] ], {l,\[Gamma]  }]},
Join[{p,l , \[Gamma] } /. sol[[2]], {sol[[1]]}]

]

evaluates the maximum of h depending on integer p.

Table[maxi[k], {k, 0, 5}] 

Answer 2

Clear["Global`*"];

h[l_, p_, γ_] := 2^(Abs[l] + 2) p! (Abs[l] + 
            p)! γ^2 (Sum[((-2)^m Gamma[m + Abs[l]/2 +         
           1])/((p - m)! (Abs[l] + m)! m! (1 + γ^2)^(m + 
                      Abs[l]/2 + 1)), {m, 0, p}])^2;

maxh[l_, 0] = Module[
   {g = SolveValues[{D[h[l, 0, γ], γ] == 0, γ > 0}, γ][[1]] // 
     Simplify}, {h[l, 0, g] // Simplify, {γ -> g}}];

maxh[l_?NumericQ, p_Integer?Positive] := 
 Maximize[{h[l, p, γ] // Simplify, γ > 0}, γ] // RootReduce

Exact inputs provides exact results

{maxh[1, 2], maxh[2, 3]}

or approximately,

% // N

(* {{0.162661, {γ -> 0.311767}}, {0.0996705, {γ -> 1.66016}}} *)

Tabulating,

Prepend[
  Flatten[
    Table[t = maxh[l, p]; {l, p, γ /. t[[2]], t[[1]]},
     {l, 0, 2, 0.25}, {p, 0, 4}],
    1] /. {rat_Rational :> N[rat], r_Root :> N[r]},
  Style[#, 14] & /@
   {"l", "p", "γ", Subscript["h", "max"]}] //
 Grid[#, Frame -> All] &