a[k_, n_] := (2^(-1 - n^k) \[Pi] Gamma[n^k])/Gamma[1 + n^k/2]^2
q[n_] := Module[{img},
  img = Rasterize[
    DensityPlot[
     Arg[a[x + I y, n]], {x, -\[Pi], \[Pi]}, {y, -\[Pi], \[Pi]},
     PlotPoints -> 50, WorkingPrecision -> 50,
     PerformanceGoal -> "Quality", MaxRecursion -> 3,
     ColorFunction -> (ColorData[{"DarkRainbow", {-Pi, Pi}}][#] &),
     ColorFunctionScaling -> False, Frame -> False,
     PlotRangePadding -> 0, ImageMargins -> 0], "Image",
    RasterSize -> 512];
  Plot3D[Abs[a[x + I y, n]], {x, -2, 2}, {y, -2, 2},
   PlotStyle -> Texture[img], PlotPoints -> 50, Mesh -> None,
   Lighting -> "Neutral", WorkingPrecision -> 50,
   PerformanceGoal -> "Quality", ClippingStyle -> White,
   MaxRecursion -> 3, ExclusionsStyle -> None,
   BoundaryStyle -> Directive[Thin, Black],
   AxesLabel -> TraditionalForm /@ {Re[z], Im[z]}, ImageSize -> 600]
  ]