  Valuable information on the dynamics of expanding fluids can be inferred from
the response of such systems to perturbations in their initial geometry. We
apply this technique in high-energy $^{96}$Ru+$^{96}$Ru and $^{96}$Zr+$^{96}$Zr
collisions to scrutinize the expansion dynamics of the quark-gluon plasma,
where the initial geometry perturbations are sourced by the differences in
deformations and radial profiles between $^{96}$Ru and $^{96}$Zr, and the
collective response is captured by the change in anisotropic flow $V_n$ between
the two collision systems. Using a transport model, we analyze how the
nonlinear coupling between lower-order flow harmonics $V_2$ and $V_3$ to the
higher-order flow harmonics $V_4$ and $V_5$, expected to scale as
$V_{4\mathrm{NL}}=\chi_4 V_2^2$ and $V_{5\mathrm{NL}}=\chi_5 V_2V_3$, gets
modified as one moves from $^{96}$Ru+$^{96}$Ru to $^{96}$Zr+$^{96}$Zr systems.
We find that these scaling relations are valid to high precision: variations of
order 20\% in $V_{4\mathrm{NL}}$ and $V_{5\mathrm{NL}}$ due to differences in
quadrupole deformation, octupole deformation, and nuclear skin modify
$\chi_{4}$ and $\chi_5$ by about 1--2\%. Percent-level deviations are however
larger than the expected experimental uncertainties and could be measured.
Therefore, collisions of isobars with different nuclear structures are a unique
tool to isolate subtle nonlinear effects in the expansion of the quark-gluon
plasma that would be otherwise impossible to access in a single collision
system.
