How can I improve the plotting speed of plot3D with region calculations?

Question

I'm looking to achieve the desired result as shown in the attached image. There may be various methods to obtain such a plot, but my current code seems to be running quite slowly. Unfortunately, my expertise only allows me to write code of this nature. I'm wondering if anyone can help me optimize my code for better performance or suggest alternative approaches to create this plot more efficiently. Thank you

Clear["Global`*"]
s0 = PauliMatrix[1];

s1 = PauliMatrix[2];
s2 = PauliMatrix[3];
s3 = PauliMatrix[4];
Ham =  kx s1 + ky s2 ;
DD = Eigenvalues[Ham];
VV = Eigenvectors[Ham];

region = RegionUnion[Disk[{0, 0}, 0.1], 
   RegionDifference[Disk[{0, 0}, 0.3], Disk[{0, 0}, 0.2]], 
   RegionDifference[Disk[{0, 0}, 0.5], Disk[{0, 0}, 0.4]], 
   RegionDifference[Disk[{0, 0}, 0.7], Disk[{0, 0}, 0.6]], 
   RegionDifference[Disk[{0, 0}, 0.9], Disk[{0, 0}, 0.8]], 
   RegionDifference[Disk[{0, 0}, 1.1], Disk[{0, 0}, 1]], 
   RegionDifference[Disk[{0, 0}, 1.3], Disk[{0, 0}, 1.2]], 
   RegionDifference[Disk[{0, 0}, 1.5], Disk[{0, 0}, 1.4]], 
   RegionDifference[Disk[{0, 0}, 1.7], Disk[{0, 0}, 1.6]], 
   RegionDifference[Disk[{0, 0}, 1.9], Disk[{0, 0}, 1.8]]];


region2 = RegionDifference[Disk[{0, 0}, 1], region];


P1 = Plot3D[DD[[1]], {kx, ky} \[Element] region, 
   PlotRange -> {-1.5, 1.5}, BoxRatios -> {1, 1, 1.1}, 
   ClippingStyle -> None, 
   ColorFunction -> (ColorData[{"SunsetColors", "Reverse"}][#3] &), 
   Mesh -> None, Axes -> True, Boxed -> False, MaxRecursion -> 2, 
   PlotPoints -> 20];
P2 = Plot3D[DD[[2]], {kx, ky} \[Element] region2, 
   PlotRange -> {-1.5, 1.5}, BoxRatios -> {1, 1, 1.1}, 
   ClippingStyle -> None, 
   ColorFunction -> (ColorData["SunsetColors"][#3] &), Mesh -> None, 
   Axes -> True, Boxed -> False, MaxRecursion -> 2, PlotPoints -> 20];
P3 = Graphics3D[{Red, Sphere[{0, 0, 0}, 0.1]}];
Show[P1, P2, BoxRatios -> {1, 1, 1.2}]

Answer 1

You can get the desired result much faster using the options MeshFunctions, Mesh and MeshShading as follows:

Plot3D[DD[[1]], {kx, -2, 2}, {ky, -2, 2}, PlotRange -> {-1.5, 1.5}, 
 BoxRatios -> {1, 1, 1.1}, ClippingStyle -> None, 
 ColorFunction -> (ColorData[{"SunsetColors", "Reverse"}][#3] &), 
 Axes -> True, Boxed -> False, MaxRecursion -> 2, PlotPoints -> 20, 
 MeshFunctions -> {#3 &}, Mesh -> {Range[-2, 2, .1]}, 
 MeshShading -> {Automatic, None}]

Do the same for DD[[1]] and DD[[2]] with association color schemes:

Show[
   Plot3D[DD[[#]], {kx, -2, 2}, {ky, -2, 2}, 
     PlotRange -> {-1.5, 1.5}, 
     BoxRatios -> {1, 1, 1.1}, 
     ClippingStyle -> None, 
     ColorFunction -> ({x, y, z} |-> ColorData[#2][z]), 
     Axes -> True, 
     Boxed -> False,
     MaxRecursion -> 2, 
     PlotPoints -> 20, 
     MeshFunctions -> {#3 &}, 
     Mesh -> {Range[-2, 2, .1]}, 
     MeshShading -> {Automatic, None}] & @@@ 
   {{1, {"SunsetColors", "Reverse"}}, 
    {2, "SunsetColors"}}, 
  PlotRange -> All]

Aside:

If you have to work with regions, you might consider Annulus to define your region:

region2 = Apply[RegionUnion] @
   Prepend[Disk[{0, 0}, .1]] @
   Map[Annulus[{0, 0}, {#, # + .1}] &, Range[.2, 1.9, .2]];

Region @ region2

Answer 2

We can BoundaryDiscretizeRegion the regions region and region2 to speed up the 3D plot.

region = BoundaryDiscretizeRegion[region];
region2 = BoundaryDiscretizeRegion[region2];

That is

Clear["Global`*"]
s0 = PauliMatrix[1];

s1 = PauliMatrix[2];
s2 = PauliMatrix[3];
s3 = PauliMatrix[4];
Ham = kx s1 + ky s2;
DD = Eigenvalues[Ham];
VV = Eigenvectors[Ham];

region = 
  RegionUnion[Disk[{0, 0}, 0.1], 
   RegionDifference[Disk[{0, 0}, 0.3], Disk[{0, 0}, 0.2]], 
   RegionDifference[Disk[{0, 0}, 0.5], Disk[{0, 0}, 0.4]], 
   RegionDifference[Disk[{0, 0}, 0.7], Disk[{0, 0}, 0.6]], 
   RegionDifference[Disk[{0, 0}, 0.9], Disk[{0, 0}, 0.8]], 
   RegionDifference[Disk[{0, 0}, 1.1], Disk[{0, 0}, 1]], 
   RegionDifference[Disk[{0, 0}, 1.3], Disk[{0, 0}, 1.2]], 
   RegionDifference[Disk[{0, 0}, 1.5], Disk[{0, 0}, 1.4]], 
   RegionDifference[Disk[{0, 0}, 1.7], Disk[{0, 0}, 1.6]], 
   RegionDifference[Disk[{0, 0}, 1.9], Disk[{0, 0}, 1.8]]];

region2 = RegionDifference[Disk[{0, 0}, 1], region];

region = BoundaryDiscretizeRegion[region];
region2 = BoundaryDiscretizeRegion[region2];

P1 = Plot3D[DD[[1]], {kx, ky} ∈ region, 
   PlotRange -> {-1.5, 1.5}, BoxRatios -> {1, 1, 1.1}, 
   ClippingStyle -> None, 
   ColorFunction -> (ColorData[{"SunsetColors", "Reverse"}][#3] &), 
   Mesh -> None, Axes -> True, Boxed -> False, MaxRecursion -> 2, 
   PlotPoints -> 20];
P2 = Plot3D[DD[[2]], {kx, ky} ∈ region2, 
   PlotRange -> {-1.5, 1.5}, BoxRatios -> {1, 1, 1.1}, 
   ClippingStyle -> None, 
   ColorFunction -> (ColorData["SunsetColors"][#3] &), Mesh -> None, 
   Axes -> True, Boxed -> False, MaxRecursion -> 2, 
   PlotPoints -> 20];
P3 = Graphics3D[{Red, Sphere[{0, 0, 0}, 0.1]}];
Show[P1, P2, BoxRatios -> {1, 1, 1.2}]