For an alternative bound (which I just saw was referred to by Sandeep in the comments while revising this), you can adapt the standard argument for the expectation of the maximum of Gaussians.
I do this below, where $M = \max_{i\in\mathcal{N}}h_j$.
I will also assume that $h_j$ has shape $k_j>0$, scale $\theta_j>0$, and pdf $f_j(x) = \frac{1}{\Gamma(k_j)\theta_j^{k_j}}x^{k_j-1}\exp(x/\theta_j)$.
One then gets that
\begin{align*}
\exp(t\mathbb{E}[M]) \leq\mathbb{E}[\exp [tM]]=\mathbb{E}[\max_{j\in\mathcal{N}} \exp(th_j)] \leq \sum_{j\in\mathcal{N}}\mathbb{E}[\exp(th_j)] \leq \sum_{j\in\mathcal{N}} (1-\theta_jt)^{-k_j},
\end{align*}
valid for $t < \frac{1}{\max_{j\in\mathcal{N}}\theta_j}$.
It follows that for such $t$ we get that
$$
\mathbb{E}[M] \leq \frac{1}{t}\ln(\sum_{j\in\mathcal{N}}(1-\theta_jt)^{-k_j}).
$$
Noting that this is an instance of the LogSumExp function, and using the linked upper bound, we get that
$$\label{1}\tag{1}
\forall t\in(0, (\max_{j\in\mathcal{N}}\theta_j)^{-1}): \mathbb{E}[M] \leq \frac{1}{t}\left(\max_{j\in\mathcal{N}}(k_j\ln\frac{1}{1-\theta_jt}) + \ln|\mathcal{N}|\right).
$$
Under the approximation
$$
\ln(\frac{1}{1-\theta_jt}) = \ln(\sum_{i = 0}^\infty (\theta_j t)^i) \stackrel{*}{\approx} \ln(1+\theta_j t) \leq \theta_jt,
$$
the above bound simplifies to
$$
\mathbb{E}[M] \leq \max_{j\in\mathcal{N}}\mathbb{E}[H_j] + \frac{\ln|\mathcal{N}|}{t}.
$$
Choosing $t$ maximally, we get the bound
$$
\mathbb{E}[M] \leq \max_{j\in\mathcal{N}}\mathbb{E}[h_j] + \max_{j\in\mathcal{N}}\theta_j\ln|\mathcal{N}|.
$$
The approximation of $\ast$ is still not justified.
It is of course formally false.
That being said, it is easy to show that something close to it is true.
For example, it is easy to see that
$$
\sum_{i = 0}^\infty x^i = 1 + x + x^2\sum_{i = 0}^\infty x^i = 1 + x + \frac{x^2}{1-x} \leq 1 + x + \epsilon(1+x) = (1+\epsilon)(1+x)
$$
provided $x^2 \leq \frac{\epsilon}{1+\epsilon}$, or in our setting $t \leq \frac{1}{\sqrt{1+\epsilon^{-1}}\max_{j\in\mathcal{N}}\theta_j}$.
Under this choice of $t$, one gets the bound
$$\label{2}\tag{2}
\mathbb{E}[M] \leq \max_{j\in\mathcal{N}}\mathbb{E}[h_j] + \max_{j\in\mathcal{N}}k_j\ln (1+\epsilon) + \sqrt{1+\epsilon^{-1}}\max_{j\in\mathcal{N}}\theta_j\ln|\mathcal{N}|,
$$
where $\epsilon\in(0,1)$.
I wouldn't be surprised if there was a better bound one could use for $\ast$ than the one I used above, but hopefully Eq. \eqref{2} is better for your purposes than ChatGPT's answer.
In particular, it has the benefit of scaling with $\ln|\mathcal{N}|$, rather than $O(\sqrt{|\mathcal{N}|})$, which may be preferable for $|\mathcal{N}|$ large.
Edit:
One can instead use the inequality
$$\tag{3}\label{3}
\ln(1/(1-x)) \leq x + x^2,
$$
valid for $x\in(0, 0.683803\dots)$. I do not have an exact value for the right endpoint, so will instead use that $2/3 \leq 0.683803\dots$ is relatively close, and Eq. \eqref{3} is valid for $x\in[0,2/3]$.
If one chooses $\epsilon = 4/5$ such that the inequality $\ln(1/(1-x)) \leq (1+\epsilon)x$ is valid on the same interval $[0, 2/3]$, it is simple to see that the approximation $x+x^2$ is much tighter, see here.
Substituting this into Eq. \eqref{1}, we get that
$$
\mathbb{E}[M] \leq \max_{j\in\mathcal{N}} (\mathbb{E}[h_j] + t\mathsf{Var}[h_j]) + \frac{\ln|\mathcal{N}|}{t},
$$
for $t\in(0, \min(2/3, (\max_{j\in\mathcal{N}}\theta_j)^{-1})]$.
So, if $2/3 \leq (\max_{j\in\mathcal{N}}\theta_j)^{-1}$, we may choose $t = 2/3$ to get the bound
$$
\mathbb{E}[M] \leq \max_{j\in\mathcal{N}}(\mathbb{E}[h_j] + (2/3)\mathsf{Var}[h_j]) + \frac{3}{2}\ln|\mathcal{N}|.
$$
If $2/3 > (\max_{j\in\mathcal{N}\theta_j})^{-1}$, we instead get that
$$
\mathbb{E}[M] \leq \max_{j\in\mathcal{N}}(\mathbb{E}[h_j] + \frac{\mathsf{Var}[h_j]}{\max_{i\in\mathcal{N}}\theta_i}) + \max_{i\in\mathcal{N}}\theta_i\ln|\mathcal{N}| \leq 2\max_{j\in\mathcal{N}}\mathbb{E}[h_j] + \max_{i\in\mathcal{N}}\theta_i\ln|\mathcal{N}|.
$$
So, this should be better than the previously-derived bound.
Perhaps there are even better variants of $\ast$ one can leverage, I don't know.
