For any \(d \geqslant 2\) , we prove that there exists an integer \(n_0(d)\) such that there exists an \(n \times n\) magic square of \(d^\text {th}\) powers for all \(n \geqslant n_0(d)\) . In particular, we establish the existence of an \(n \times n\) magic square of squares for all \(n \geqslant 4\) , which settles a conjecture of Várilly-Alvarado. All previous approaches had been based on constructive methods and the existence of \(n \times n\) magic squares of \(d^\text {th}\) powers had only been known for sparse values of n. We prove our result by the Hardy-Littlewood circle method, which in this setting essentially reduces the problem to finding a sufficient number of disjoint linearly independent subsets of the columns of the coefficient matrix of the equations defining magic squares. We prove an optimal (up to a constant) lower bound for this quantity.