Return to Answer

I could perhaps track down the problem a little better: it seems to be the SphericalHarmonicY at angles [Theta] below Pi/4 and above Pi - Pi/4, which finally provide extremely high values in dat and hides the details of the other values (at least for lmax = 64).

dat = ParallelTable[ fieldN[\[Theta], \[Phi]], {\[Theta], Pi/4., Pi - Pi/4., Pi/128.}, {\[Phi], 0, 2 Pi, 2 Pi/256.} ]; 

I could perhaps track down the problem a little better: it seems to be the SphericalHarmonicY at angles θ below Pi/4 and above Pi - Pi/4, which finally provide extremely high values in dat and hides the details of the other values (at least for lmax = 64).

dat = ParallelTable[ fieldN[θ, ϕ], {θ, Pi/4., Pi - Pi/4., Pi/128.}, {ϕ, 0, 2 Pi, 2 Pi/256.} ]; 

EDIT

I could perhaps track down the problem a little better: it seems to be the SphericalHarmonicY at angles [Theta] below Pi/4 and above Pi - Pi/4, which finally provide extremely high values in dat and hides the details of the other values (at least for lmax = 64).

For example the table:

dat = ParallelTable[ fieldN[\[Theta], \[Phi]], {\[Theta], Pi/4., Pi - Pi/4., Pi/128.}, {\[Phi], 0, 2 Pi, 2 Pi/256.} ]; 

seems to me to produce a reasonable chart.

A part from efficiency, I noticed that with the definition below, lmax >= 64 works:

Clear[field]; field[\[Theta]_, \[Phi]_] := Chop@ Total[ Table[ alms[l, m] SphericalHarmonicY[l, m, \[Theta], \[Phi]], {l, 0, lmax}, {m, -l, l} ] , 2 ]; nn = 4.; dat = ParallelTable[ field[\[Theta], \[Phi]], {\[Theta], 0, Pi, Pi/nn}, {\[Phi], 0., 2 Pi, 2 Pi/nn/2} ]; Round[Re[dat], 0.01] 

{ {0.51, 0.51, 0.51, 0.51, 0.51, 0.51, 0.51, 0.51, 0.51},

{0.02, -1.54, -1.48, -0.47, -2.29, 0.61, 2.05, 1.82, 0.02},

{0.24, -1.61, 0.44, 0.32, -0.55, 0.65, -1.12, -0.08, 0.24},

{-0.34, 0.63, -0.54, 2.08, -0.72, 1.09, -2.09, -1.3, -0.34},

{0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3} }

which looks better and produces a reasonable plot.

So, I tested also:

Clear[field]; field[\[Theta]_, \[Phi]_] := Sum[ alms[l, m] SphericalHarmonicY[l, m, \[Theta], \[Phi]], {l, 0, lmax}, {m, -l, l} ]; 

and this too seems to me to work.

Consequently, it appears to me that Compile over Sum is doing some inappropriate manipulation.

A part from efficiency, I noticed that with the definition below, lmax >= 64 works:

Clear[field]; field[θ_, ϕ_] := Chop@ Total[ Table[ alms[l, m] SphericalHarmonicY[l, m, θ, ϕ], {l, 0, lmax}, {m, -l, l} ] , 2 ]; nn = 4.; dat = ParallelTable[ field[θ, ϕ], {θ, 0, Pi, Pi/nn}, {ϕ, 0., 2 Pi, 2 Pi/nn/2} ]; Round[Re[dat], 0.01] 

{ {0.51, 0.51, 0.51, 0.51, 0.51, 0.51, 0.51, 0.51, 0.51},

{0.02, -1.54, -1.48, -0.47, -2.29, 0.61, 2.05, 1.82, 0.02},

{0.24, -1.61, 0.44, 0.32, -0.55, 0.65, -1.12, -0.08, 0.24},

{-0.34, 0.63, -0.54, 2.08, -0.72, 1.09, -2.09, -1.3, -0.34},

{0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3} }

which looks better and produces a reasonable plot.

So, I tested also:

Clear[field]; field[θ_, ϕ_] := Sum[ alms[l, m] SphericalHarmonicY[l, m, θ, ϕ], {l, 0, lmax}, {m, -l, l} ]; 

and this too seems to me to work.

Consequently, it appears to me that Compile over Sum is doing some inappropriate manipulation.