  Quantum fluctuations plays an essential role in forming the collective flow
of hadrons observed in relativistic heavy-ion collisions. Event-by-event
fluctuations of the collective flow can arise from various sources, such as the
fluctuations in the initial geometry, hydrodynamic expansion, hadronization,
and hadronic evolution of the nuclear matter, while the exact contribution from
each source is still an open question. Using a (3+1)-dimensional relativistic
hydrodynamic model coupled to a Monte-Carlo Glauber initial condition,
Cooper-Frye particlization and a hadronic transport model, we explore the
system size and shape dependences of the collective flow fluctuations in Au+Au,
Cu+Au, and O+O collisions at $\sqrt{s_\mathrm{NN}}=200$~GeV. The particle
yields, mean transverse momenta, 2-particle and 4-particle cumulant elliptic
flows ($v_2\{2\}$ and $v_2\{4\}$) from our calculation agree with the currently
existing data from RHIC. Different centrality dependences of the flow
fluctuations, quantified by the $v_2\{4\}/v_2\{2\}$ ratio, are found for
different collision systems due to their different sizes and shapes. By
comparing $v_2\{4\}/v_2\{2\}$ between different hadron species, and comparing
$v_2\{4\}/v_2\{2\}$ to the initial state geometric fluctuations quantified by
the cumulant eccentricity ratio $\varepsilon_2\{4\}/\varepsilon_2\{2\}$, we
find that while the initial state fluctuations are the main source of the $v_2$
fluctuations in large collision systems, other sources like nonlinear
hydrodynamic response, hadronization, and hadronic afterburner can
significantly affect the $v_2$ fluctuations in small systems.
