  We explore theoretical uncertainties in the hydrodynamic description of
relativistic heavy-ion collisions by examining the full non-linear causality
conditions and quantifying the second-order transport coefficients' role on
flow observables. The causality conditions impose physical constraints on the
maximum allowed values of inverse Reynolds numbers during the hydrodynamic
evolution. Including additional second-order gradient terms in the
Denicol-Niemi-Moln\'{a}r-Rischke (DNMR) theory significantly shrinks the casual
regions compared to those in the Israel-Stewart hydrodynamics. For Au+Au
collisions, we find the variations of flow observables are small with and
without imposing the necessary causality conditions, suggesting a robust
extraction of the Quark-Gluon Plasma's transport coefficients in previous
model-to-data comparisons. However, sizable sensitivity is present in small
p+Au collisions, which poses challenges to study the small systems'
collectivity.
