The following text is written by Christian Boyer and is coming from his web site www.multimagie.com/indexengl.htm , Thanks for his authorization to use his text.

Multimagic formulas

    How do we know what should be the theoretical values of sums of rows, columns and diagonals of the multimagic squares?
    Start with the single magic square. It is known, and easy to demonstrate (…except for an 8-year old child like Gauss who is supposed to have demonstrated at this age at school…), that the sum of the first integer numbers from 1 to N is:
                 N(N+1)/2
    A n-th order magic square contains all the numbers from 1 to n². So, the sum of the first integer numbers from 1 to n² is:
                 n²(n²+1)/2
    Since this square has n rows, it is sufficient to divide by n in order to know the magic sum S1 of each row, column or diagonal, so:
          S1 = n(n²+1)/2
    For the bimagic square, the sum of the squares of the first integers from 1 to N is:
                            N(N+1)(2N+1)/6
    Replacing N by n², then dividing by n, we get the sum S2 of each row, column or diagonal, so:
          S2 = n(n²+1)(2n²+1)/6 = S1*(2n²+1)/3
  For the trimagic square, the sum of the cubes:
          S3 = n * S1²
    For the tetramagic, Fermat gave as early as 1636, in a letter sent to Roberval, the solution for the sum of the integers from 1 to N to the 4th-power. Here is the Fermat's formula after having replaced N by n², then dividing by n:
              S4 = ((4n² + 2)*n*S1² - S2) / 5
   Replacing (4n² +2) * S1 by 6S2, we can write an easier formula than Fermat:
          S4 = S2 * (6n * S1 - 1) / 5 
    For the pentamagic (sum of 5th-power):
          S5 = (3n*S2² - S3) / 2
 
    If we prefer to develop Sp, it gives the following generous table, which permits us to calculate the sums of the future hexa, hepta or octomagic squares:

    In order to calculate Sp when the n-th order p-multimagic square contains the integers from 0 to n²-1 (instead of 1 to n²), it is sufficient to put a - sign (instead of the + red ) in the 2nd column above. You may remark that the first column of Sp is equal to (1/(p+1)) n^(2p+1). This feature has been known since the middle of the XVIIth century, thanks to Pascal who has dedicated a Treatise on the subject. Jacques Bernoulli later took an interest in this question of the sum of integer powers. After his Ars conjectandi published after his death at the beginning of the XVIIIth century, the numbers used in this development are now called Bernoulli numbers (1/6, -1/30, 1/42, -1/30, 5/66, …).
    Here are the formulas for some n-orders of p-multimagic squares, containing numbers from 1 to n² . Here we will find the sums of the Pfeffermann's 8th-order bimagic square, Pfeffermann's 9th-order bimagic square, and William Benson's 32-th order trimagic square:

And here are the formulas of some n-orders of p-multimagic squares, containing numbers from 0 to n²-1. Here we will find the sums of our record 512th-order tetramagic and 1024th-order pentamagic squares:

Bibliography
    For more details about sum of powers, see the excellent chapter 14 Sommation des puissances numériques, book II of the Théorie des Nombres by Edouard Lucas , Librairie Blanchard, Paris.