How to find the expectation $mathbb{E}left[ a mathcal{Q} left( sqrt{b } gamma right) right]$?

It's not too hard to set this up. I had no luck with the more complex four parameter GammaDistribution but fortunately you only mentioned the two parameter version:

q[x_] = 1/(2 π) Integrate[Exp[-u^2/2], {u, x, ∞}]
G = GammaDistribution[κ, θ];
gpdf = PDF[G, y]
result = a*Expectation[q[Sqrt[b] y], y [Distributed] G]

Result:

2^(-3 - κ/2) a b^(-(1/2) - κ/
  2) θ^(-1 - κ) (2 Sqrt[2] Sqrt[
    b] θ HypergeometricPFQRegularized[{(1 + κ)/
      2, κ/2}, {1/2, (2 + κ)/2}, 1/(
     2 b θ^2)] - κ HypergeometricPFQRegularized[{(
      1 + κ)/2, (2 + κ)/2}, {3/2, (3 + κ)/2}, 1/(
     2 b θ^2)])

Let's do a quick check to make sure it returned something reasonable. We'll generate some gamma distributed random numbers, apply the function and get the mean. Then also use these constants with our result obtained earlier to confirm it's a close match to this numerical experiment:

qn[x_?NumericQ] := 1/(2 π) NIntegrate[Exp[-u^2/2], {u, x, ∞}]
With[{κ = 1, θ = 2, a = 4, b = 3},
 rands = RandomVariate[GammaDistribution[κ, θ], 5000];
 Mean[a*qn[Sqrt[b]*#] & /@ rands]
]
(* result: 0.155478 *)

N[result /. {κ -> 1, θ -> 2, a -> 4, b -> 3}]
(* result: 0.15502 *)

Looks about right!

---

Your version of the Q-function has 1/(2 π). I think this should be 1/Sqrt[2 π] instead. In which case the result changes:

q[x_] = 1/Sqrt[2 π] Integrate[Exp[-u^2/2], {u, x, ∞}]
G = GammaDistribution[κ, θ];
gpdf = PDF[G, y]
result = a*Expectation[q[Sqrt[b] y], y [Distributed] G]

Result:

2^(-(5/2) - κ/2) a b^(-(1/2) - κ/
  2) Sqrt[π] θ^(-1 - κ) (2 Sqrt[2] Sqrt[
    b] θ HypergeometricPFQRegularized[{(1 + κ)/
      2, κ/2}, {1/2, (2 + κ)/2}, 1/(
     2 b θ^2)] - κ HypergeometricPFQRegularized[{(
      1 + κ)/2, (2 + κ)/2}, {3/2, (3 + κ)/2}, 1/(
     2 b θ^2)])