Standard error with maximum likelihood estimation ASTM E2238 describes the prediction of the largest inclusion in a polished steel surface of 150,000 mm² based on a Gumbel distribution of the largest inclusions measured in each of 24 polished microsections (150 mm² each). The location and scale parameters $\lambda$ and $\delta$ of the distribution are determined using the maximum likelihood method. Then, section 6.9.5 of the standard states:

The standard error, SE, for any inclusion of length $x$ based upon the
ML method is
$SE(x)=\delta_{ML}\cdot \sqrt{(1.109 + 0.514\cdot y+0.608\cdot y^{2})/n}$

where $y = (x-\lambda )/ \delta$ and $n$ is the number of measurements (24 in this standard).
(I can't upload the complete ASTM standard for copyright reasons, but there are a number of openly accessible sources that use this method, e.g. here, p. 22 here or eq. (6) here. They all use, without further comment, this formula).
My question is, where do these coefficients $1.109$, $0.514$ and $0.608$ come from?
It seems to me that the answer should be obvious to the professional as everyone uses that formula and no-one wonders about its orgin, so I'd like apologize in advance if this is a dumb question but I'm not a statistician. A reference to any textbook or other source would be highly welcomed.
Additional information as requested in comment: The SE formula applies to all values (including the largest value determined by the procedure) and is used to draw the confidence bands $95\%CI=\pm 2\cdot SE(x)$ as in fig 1. Yes, the cdf used is $F(x) = e^{-e^{-(x-\lambda)/\delta}}$. The largest inclusion is calculated for a return period of $T = 1000$ (i.e. $A_{ref}=1000 A_{0}$ where $A_{0}=150 mm^{2}$) by $-\delta_{ML}\cdot ln(-ln((T-1)/T))+\lambda_{ML}$.
 A: Let $\hat{\delta}_{ML}$ and $\hat{\lambda}_{ML}$ be the maximum likelihood estimates of $\delta$ and $\lambda$.
Also, let $\gamma$ denote the Euler-Mascheroni constant (approximately 0.5772).
The asymptotic covariance matrix of $\left(\hat{\lambda}_{ML},\hat{\delta}_{ML} \right)'$ is:
$$\frac{\delta^{2}}{n \pi^2}\begin{bmatrix}6(\gamma-1)^2+\pi^2 & 6(1-\gamma)\\6(1-\gamma) & 6\end{bmatrix}$$
That is found from the inverse of the Fisher Information matrix.
The Fisher information for a single observation $X$ is found by finding the expected value of the second derivatives $\frac{d^2}{d\lambda^2},\frac{d^2}{d\lambda d\delta},\frac{d^2}{d\delta^2}$ of $\log f(x)$; then finding the expected value of those where you change $x$ to a random variable $X$ having that Gumbel distribution. Then, take the negative of that matrix of expected values of second derivatives. Those are not trivial to find, but I did it in Mathematica. The Fisher information adds from independent samples, so you multiply by $n$ to get the Fisher information for the entire sample.
fx = D[E^(-E^(-(x - lam)/del)), x];
InputForm[
 Inverse[-n Assuming[Element[lam, Reals] && del > 0, 
    Integrate[
      {{D[Log[fx], lam, lam], D[Log[fx], lam, del]},
       {D[Log[fx], lam, del], D[Log[fx], del, del]}} fx, 
         {x, -Infinity, Infinity}]]]]

returns
{{(del^2*(6*(-1 + EulerGamma)^2 + Pi^2))/(n*Pi^2), 
  (6*del^2*(1 - EulerGamma))/(n*Pi^2)}, 
 {(6*del^2*(1 - EulerGamma))/(n*Pi^2), (6*del^2)/(n*Pi^2)}}

Now, to find the approximate variance of $x=\hat{\delta}y+\hat{\lambda}$, you use
$$(1,y)\frac{\hat{\delta}_{ML}^{2}}{n \pi^2}\begin{bmatrix}6(\gamma-1)^2+\pi^2 & 6(1-\gamma)\\6(1-\gamma) & 6\end{bmatrix}(1,y)'$$
This is equal to $\frac{\hat{\delta}_{ML}^{2}}{n}\left(\frac{6-12\gamma+6\gamma^2+\pi^2}{\pi^2}+\frac{12-12\gamma}{\pi^2}y+\frac{6}{\pi^2}y^2 \right)\approx\frac{\hat{\delta}_{ML}^{2}}{n}(1.10866+0.514044y+0.607927y^2)$
The estimated standard error is the square root of the estimated variance:
$$\hat{\delta}_{ML}\sqrt{\frac{1.10866+0.514044y+0.607927y^2}{n}}$$
