  Quasisymmetry (QS) is a hidden symmetry of the magnetic field strength, $B$,
that confines charged particles effectively in a three-dimensional toroidal
plasma equilibrium. Here, we show that QS has a deep connection to the
underlying symmetry that makes solitons possible. Our approach uncovers a
hidden lower dimensionality of $B$ on a magnetic flux surface, which could make
stellarator optimization schemes significantly more efficient. Recent numerical
breakthroughs have yielded configurations with excellent volumetric QS and
surprisingly low magnetic shear. Our approach elucidates why the magnetic shear
is low in these configurations. Furthermore, we deduce an upper bound on the
maximum toroidal volume that can be quasisymmetric and verify it for the
Landreman-Paul precise quasiaxisymmetric (QA) stellarator configuration. In the
neighborhood of the outermost surface, we show that the B approaches the form
of the 1-soliton reflectionless potential. We present three independent
approaches to demonstrate that quasisymmetric $B$ is described by well-known
integrable systems such as the Korteweg-de Vries (KdV) equation. The first
approach is weakly nonlinear multiscale perturbation theory, which highlights
the crucial role that magnetic shear plays in QS. We show that the
overdetermined problem of finding quasisymmetric vacuum fields admits solutions
for which the rotational transform is not free but highly constrained. We
obtain the KdV equation (and, more specifically, Gardner's equation for certain
choices of parameters). Our second approach is non-perturbative and based on
ensuring single-valuedness of $B$, which directly leads to its Painlev\'e
property shared by the KdV equation. Our third approach uses machine learning,
trained on a large dataset of numerically optimized quasisymmetric
stellarators. We robustly recover the KdV (and Gardner's) equation from the
data.
