Find the identity sandpile
Mathematica, 177 157 135 133 bytes
Colorize[f=BlockMap[⌊{l={0,1,0},1-l,l}#/4⌋~Total~2+#[[2,2]]~Mod~4&amp;,#~ArrayPad~1,{3,3},1]&amp;~FixedPoint~#&amp;;k=Table[6,#,#];f[k-f@k]]&amp;

Takes a number n. The output is the identity sandpile. 0 is black, 1 is light gray, 2 is magenta, and 3 is blue-gray.
Sadly, Mathematica does not have a builtin for this...
Uses the algorithm stated in Scott Corry and David Perkinson's paper.
Takes 91.7 seconds on my 5-year-old laptop to calculate the 50x50 identity sandpile. I am confident that a reasonable modern desktop computer is more than 50% faster. (I also have a way faster code at the end).
Explanation
f= ...

Define function f (the input is a sandpile matrix): a function that ...
BlockMap[ ... ]~FixedPoint~#&amp;

... repeats the BlockMap operation until the output does not change. BlockMap operation: ...

#~ArrayPad~1

... pad the input array with one layer of 0 ...
{3,3},1

... partition it into 3x3 matrices, with offset 1 ...
⌊{l={0,1,0},1-l,l}#/4⌋~Total~2+#[[2,2]]~Mod~4&amp;

... and for each partition, add the number of sand grains toppled onto the center cell and the center cell value mod 4.
i.e. the output of f is the stabilized version of the input.

k=Table[6,#,#]

Define k as an n by n array of 6s.
f[k-f@k]]

Compute f(k - f(k)).
Colorize[ ... ]

Apply colors to the result.
Faster version (142 bytes)
Colorize[r=RotateLeft;a=ArrayPad;f=a[b=#~a~1;b+r[g=⌊b/4⌋,s={0,1}]+g~r~-s+r[g,1-s]+r[g,s-1]-4g,-1]&amp;~FixedPoint~#&amp;;k=Table[6,#,#];f[k-f@k]]&amp;

Same code, but uses built-in list rotation instead of BlockMap. Computes n=50 in 4.0 seconds on my laptop.