What is the position of the maximal value of this bell-shaped function? Consider the following function :
$$\tag{1}
f(v) = v^{\frac{d}{2}} \int_v^{\infty} u^{\alpha \,-\, \smash{\frac{d}{2}} \,-\, 1} \; e^{-\, \alpha \, u} \; du,
$$
where $d \le 6$ and $\alpha > 0$ are two positive constants (parameters).  Notice the lower limit of the integral : $v$ is a variable.  The plot of this function (for $0 \le v < \infty$) shows an almost bell-shaped curve, so it has a single maximal value.
Now, I would like to find the position $v = v_0(\alpha, d)$ of the maximal value of this function.  I need an analytical expression for $v_0(\alpha, d)$, probably an approximation.
$$\tag{2}
v_0(\alpha, d) \approx \; ?
$$
From the graph of $f(v)$, I know that $v_0 \propto d$ (maybe with some exponent).
Take note that the derivative of function (1), set to 0, give this relation :
$$\tag{3}
f(v_0) \equiv f_{\text{max}} = \frac{2}{d} \; v_0^{\alpha} \; e^{-\, \alpha \, v_0},
$$
where $v_0 \equiv v_0(\alpha, d)$ is the position of the max value of (1).
Someone knows a method to find the function (2) ?

EDIT : Function (1) describes the "deformation" of a black body luminosity caused by the expansion of space in a cosmology model.  From Wien's and Planck's laws, $\alpha$ should be around 3 (depending on the presence of gaz and dust).  $d$ describes the kind of fluid contained in the cosmological model.  We have $d = 3$ for a dust filled universe, and $d = 4$ for a radiation universe.  $d = 0$ for an empty universe with a cosmological constant.  Since $\frac{d}{2} + 1$ gives 2.5 for dust and 3 for radiation, I suspect that $\alpha$ should be close to $\frac{d}{2} + 1$ (while it is an independant parameter).  In this special case, it is easy to explicitely evaluate the integral in (1) and we get the special case
$$\tag{4}
v_0(d) = \frac{d}{d + 2}, \quad \text{if $\alpha = \frac{d}{2} + 1$}.
$$
 A: This is a partial solution only and needs a correction for more precision in the result.
I wrote a Mathematica code to show the curve of function (1), with sliders to interactively change the two parameters $\alpha$ and $d$.  When $\alpha = \frac{d}{2} + 1$, the exact solution to the problem is
$$\tag{1}
v_0(d) = \frac{d}{d + 2}, \quad \text{if $\alpha = \frac{d}{2} + 1$}.
$$
When $\alpha \ne \frac{d}{2} + 1$, the Mathematica graphics shows that we still have $v_0 \approx \frac{d}{d \,+\, 2}$, for all values in the ranges $0 \le \alpha \le 10$ and $0 \le d \le 6$.
When $\alpha$ is very different than $\frac{d}{2} + 1$, formula (1) misses the true value by a small amount only.  So it needs a small correction like
$$\tag{2}
v_0(\alpha, d) \approx \frac{d}{d + 2} \; g(\alpha, d),
$$
where $g(\alpha, d) \approx 1 + \delta(\alpha, d)$, or maybe something like
$$\tag{3}
v_0(\alpha, d) \approx \frac{d}{d + 2} + h(\alpha, d),
$$
where $h(\alpha, d)$ is a small correction.
Any idea about how to find the correction $\delta(\alpha, d)$ or  $h(\alpha, d)$?
Here's the Mathematica code for you to see how good is the partial solution (1) above :
Clear["Global`*"]

function[v_, a_, d_] := v^(d/2) NIntegrate[u^(a - d/2 - 1) Exp[-a u], {u, v, Infinity}]

Xpic[a_, d_] := d/(d + 2)
Ypic[a_, d_] := function[d/(d + 2), a, d]

Pic[a_, d_] := Graphics@{Black, PointSize -> 0.01, Point[{Xpic[a, d], Ypic[a, d]}]}

Manipulate[
    Show[{Plot[{function[v, a, d]}, {v, 0, 4}, PlotStyle -> Directive[Red, Thick]], Pic[a, d]},
    Frame -> True,
    AspectRatio -> 1,
    PlotRange -> Automatic,
    GridLines ->  Automatic,
    GridLinesStyle -> Directive[Dashed, GrayLevel[0.75]],
    ImageSize -> 500
    ],
    {{a, 3, Style["\[Alpha]", 10]}, 0.5, 10, 0.01, ImageSize -> Large},
    {{d, 3, Style["d", 10]}, 0, 6, 0.01, ImageSize -> Large},
    ControlPlacement -> Bottom,
    FrameMargins -> None
]

A: Taking logarithmic derivatives, the maximum is attained when
$$ \frac{d}{2v} - \frac{\mathrm{e}^{-v \alpha} \alpha (v \alpha)^{\alpha - (d/2 + 1)}}{\Gamma(\alpha - d/2, v\alpha)} = 0  $$
The form of the first summand suggests finding the (Laurent) series expansion of the second term around $v = \infty$.  (I.e., find a Taylor series in powers of $1/v$.)  This series is
$$  \alpha - \frac{\alpha - (d/2-1)}{v} + \frac{2\alpha -(d+2)}{2 \alpha v^2} + ... $$
Solving $\dfrac{d}{2v} - \left( \alpha - \frac{\alpha - (d/2-1)}{v} \right) = 0$ for $v$, we get $\frac{\alpha - 1}{\alpha}$.  If we substitute $\frac{d}{2}+1$ for $\alpha$ in this, we get $\frac{d}{d+2}$, as expected from your prior work.
Solving $\dfrac{d}{2v} - \left( \alpha - \frac{\alpha - (d/2-1)}{v} + \frac{2\alpha -(d+2)}{2 \alpha v^2} \right) = 0$ for $v$, we get two roots: 
$$  v_\pm = \frac{\alpha - 1 \pm \sqrt{2d + (\alpha - 5)(\alpha-1)}}{2\alpha}  \text{.}  $$
If we substitute $\frac{d}{2}+1$ for $\alpha$ in these, we get $0$ for $v_-$ (corresponding to the zero at $v=0$ for certain parameter values) and $\frac{d}{d+2}$ for $v_+$, as expected.  This suggests $v_+$ is an improved estimate for your $v_0(\alpha,d)$.
(Note:  I haven't graphed, bounded, or sanity checked the above, and it looks like I won't have time for several days.  I wouldn't be very surprised to find I have a sign error (or worse) in transcription.  I'm considering making this Community Wiki so others can improve it while I'm otherwise engaged.  If someone finds a meaningful improvement, correction, or unrecoverable defect, say so in comments and I'll flip the switch (or withdraw the answer).)
