Topology optimization for time-dependent systems is computationally intensive because each design iteration requires solving and storing forward and adjoint solutions over the entire time horizon. We address this challenge with proposing a high-order space--time spectral element method for transient heat conduction, treating time as an additional dimension to solve all time steps simultaneously. High-order spectral element discretization in both space and time provides improved solution accuracy and convergence rates, yielding more accurate optimized structures with fewer space--time degrees of freedom (DoFs) respect to finite elements and finite difference methods. Our space--time formulation uses summation-by-parts (SBP) operators and imposes interface, boundary, and initial/terminal conditions weakly via simultaneous approximation terms (SAT), resulting in a stable monolithic space–time scheme on heterogeneous domains. Stability is proven under specific conditions on the SAT parameters, scaled with the spatial mesh resolution and material properties. We compute design sensitivities using a discrete space–time adjoint scheme that is dual-consistent with the primal SBP–SAT scheme. Dual consistency ensures that the discrete adjoint consistently approximates the continuous dual problem and, under standard smoothness assumption, yields superconvergent functional estimates. We validate the proposed approach by comparing the optimized designs against an independently obtained reference solution, and we present performance results including time-to-solution and accuracy-cost benchmarks against a conventional low-order time-marching approach and an alternative all-at-once solver. Our space--time spectral element method achieves high accuracy and maintains stability with significantly fewer DoFs, reducing time-to-solution and memory compared with an alternative all-at-once solver. This improvement makes large-scale topology optimization of unsteady thermal systems feasible in practice.