I am trying to prove that there exists a unique orthonormal matrix $T \in \mathbb{R}^{p \times p}$ (up to permutation and sign changes) maximizing the quartimax criterion as $$ \underset{T}{\arg\max} \,\, \sum_{i=1}^n \sum_{j=1}^p (AT)_{ij}^4 $$
where $A \in \mathbb{R}^{n \times p}$ has orthonormal columns ($A^TA = I)$ for $n > p$.
I checked that the quartimax criterion is maximized when each row of $AT$ has at most one nonzero element (simple structure) or when $A$ has full rank with $n=p$.
I'm unsure how to proceed for the general case where $A$ does not have that simple structure with $n>p$.
$ \def\bbR#1{{\mathbb R}^{#1}} \def\o{{\tt1}}\def\p{\partial} \def\LR#1{\left(#1\right)} \def\BR#1{\Big(#1\Big)} \def\bR#1{\big(#1\big)} \def\trace#1{\operatorname{Tr}\LR{#1}} \def\skew#1{\operatorname{Skew}\LR{#1}} \def\cay#1{\operatorname{Cayley}\LR{#1}} \def\vc#1{\operatorname{vec}\LR{#1}} \def\qiq{\quad\implies\quad} \def\grad#1#2{\frac{\p #1}{\p #2}} \def\hess#1#2#3{\frac{\p^2 #1}{\p #2\,\p #3}} \def\c#1{\color{red}{#1}} \def\CLR#1{\c{\LR{#1}}} \def\fracLR#1#2{\LR{\frac{#1}{#2}}} \def\gradLR#1#2{\LR{\grad{#1}{#2}}} $For typing convenience, define the variables $$\eqalign{ M &= AT \\ M^{\odot 3} &= M\odot M\odot M \\ Q &= 4\,A^TM^{\odot 3} \\ }$$ where the symbol $(\odot)$ denotes Hadamard products/powers.
It will also be convenient to use $J$ to denote the all-ones matrix (which is the same size as $A)\,$ and a colon to denote matrix inner product $$\eqalign{ A:M &= \sum_{i=1}^n\sum_{j=1}^p A_{ij}M_{ij} \;=\; \trace{A^TM} \\ A:A &= \|A\|^2_F \\ }$$ Introduce an unconstrained matrix variable $U\in\bbR{p\times p}$ and use it to construct $T$ $$\eqalign{ S &= \skew{U} &\doteq \fracLR{U-U^T}2 \\ T &= \cay{S} &\doteq \LR{I+S}^{-1}\LR{I-S} \\ S &= \cay{T} \\ }$$ These definitions have some interesting consequences $$\eqalign{ \LR{I-S}^{-1}T &= \LR{I+S}^{-1} \\ I+T &= \LR{I+S}^{-1}\BR{(I+\c{S})+(I-\c{S})} \\ &= 2\;\LR{I+S}^{-1} \\ dT &= -2\;\LR{I+S}^{-1}\,dS\,\LR{I+S}^{-1} \\ }$$ Now we are ready to calculate the gradient $\bR{{\rm wrt}\:U}$ of the quartimax function $$\eqalign{ \phi &= J:M^{\odot 4} \\ d\phi &= J:\LR{4M^{\odot 3}\odot dM} \\ &= {4 M^{\odot 3}}:\LR{A\;dT} \\ &= \c{Q}:dT \\ &= -2\,Q:\LR{I+S}^{-1}\,dS\,\LR{I+S}^{-1} \\ &= -2\,\LR{I-S}^{-1}Q\LR{I-S}^{-1}:\skew{dU} \\ &= -2\,\skew{\LR{I-S}^{-1}Q\LR{I-S}^{-1}}:{dU} \\ \grad{\phi}{U} &= -2\,\skew{\LR{I-S}^{-1}Q\LR{I-S}^{-1}} \\ &= {\LR{I+S}^{-1}Q^T\LR{I+S}^{-1}\;-\;\LR{I-S}^{-1}Q\LR{I-S}^{-1}} \\ &= \LR{I-S}^{-1}\LR{TQ^TT-Q}\LR{I-S}^{-1} \\ }$$ Setting the gradient to zero yields the first order condition (FOC) $$\eqalign{ Q &= TQ^TT \\ }$$ Since $(Q,T)$ are functions of $U$ this represents a complicated function of $U$ which must be solved for the optimal matrix $U_*$ which will satisfy the FOC. This will be very difficult.
Conversely, a numerical solution will be very simple. Since you have an exact expression for the gradient $\gradLR{\phi}{U}$ you can employ any gradient ascent algorithm to calculate the optimal $U_*$
Analyzing the numerical solutions for an ensemble of $(n,p,A)$ parameters might reveal something about the structure of $T$.