### Abstract:

Recent progress in the hydrodynamic simulation of heavy-ion collisions have characterized the fluctuating initial state and the viscous corrections to the corresponding collective flow. These fluctuations naturally explain the " ridge" and " shoulder" structure of the measured two-particle correlation functions at RHIC and the LHC. We introduce a cumulant expansion for analyzing the azimuthal fluctuations in the initial state. The cumulant definitions systematically describe the azimuthal anisotropy order by order. In particular, the dipole asymmetry $\epsilon_1$ appears at third order in the expansion, and the response to this initial fluctuation produces a radipity even dipole flow $v_1$, which has been subsequently confirmed by experiment. In addition, the cumulant expansion organizes the study of the nonlinear response to the initial conditions. The linear and nonlinear response coefficients to a given initial state were calculated with ideal and viscous hydrodynamic simulations. The collective flow is generated either linearly or nonlinearly, and the relative contribution of these two mechanisms to the observed flow pattern is calculated as a function of harmonic order, collision centrality, and the shear viscosity. For non-central collisions and high harmonic orders $\nge 4$, the nonlinear response is the dominant mechanism. This result is also seen in event-by-event hydrodynamic simulations. Using the cumulant expansion and the corresponding linear and nonlinear response coefficients, we determine the event plane correlations and compare to first measurements of this type. The observed event plane correlations are rooted in the initial state participant plane correlations, but a large fraction of the observed correlations are determined by harmonic mixing during the bulk expansion. Viscous corrections to the hydrodynamic formulation of collective flow are reflected in hydrodynamic equations of motion, as well as the correction to the distribution function at freeze-out, $\delta f(x,\vec p)$. Taking into consideration of the connection between kinetic theory and hydrodynamics, the consistent form of $\delta f(x,\vec p)$ is determined through second order in the gradient expansion, $\delta f(x,vec p) = \delta f_{(1)}(x,\vec p) + \delta f_{(2)}(x,\vec p) + \ldots$, The effect of $\delta f_{(2)}(x,\vec p)$ is found to be small for lower order harmonic flows $\nleq3 $, but is significant for the higher harmonics, $\ngeq 4$. In addition, the effect of $\delta f_{(2)}(x,\vec p)$ is small for nucleus-nucleus collisions at the LHC, but is more pronounced at RHIC and in small collision systems such as proton-nucleus collisions. $\delta f_{(2)}(x,\vec p)$ delineates the domain of applicability of viscous hydrodynamics, and systematically improves the hydrodynamic description of heavy ion collisions.