We develop a robust structure theory for multiple ergodic averages of commuting transformations along Hardy sequences of polynomial growth. We then apply it to derive a number of novel results on joint ergodicity, recurrence and convergence. In particular, we prove joint ergodicity for (a) pairwise independent Hardy sequences and weakly mixing transformations, (b) strongly independent Hardy sequences and ergodic transformations, (c) strongly irrationally independent Hardy sequences and totally ergodic transformations. We use these joint ergodicity results to provide new recurrence results for multidimensional patterns along strongly independent Hardy sequences, showing for instance that all subsets of $\mathbb{Z}^2$ of positive upper density contain patterns of the form $$ (m_1, m_2),\; (m_1 + \lfloor n^{\sqrt{2}}\rfloor, m_2),\; (m_1, m_2 + \lfloor n^{\sqrt{2}} + n^{1/2}\rfloor).$$ Last but not least, we positively resolve the joint ergodicity classification problem for pairwise independent Hardy sequences, of which the aforementioned families are special cases. While building on recent technical advances (e.g. PET coefficient tracking schemes and joint ergodicity criteria), our work introduces a number of technical developments of its own. We construct a suitable generalization of Host-Kra and box seminorms that quantitatively control ergodic averages along Hardy sequences. We subsequently use them to obtain Host-Kra seminorm estimates for averages along all pairwise independent Hardy sequences. Furthermore, we develop an ergodic version of the quantitative concatenation argument that has recently found extensive use in combinatorics, number theory and harmonic analysis. Lastly, we obtain new simultaneous Taylor approximations for Hardy sequences, a crucial ingredient to deal with the aforementioned classes of Hardy sequences.