We use an angular momentum approach to study the Cassini states (CS) of large natural satellites such as the Galilean satellites and Titan. Unlike classical approaches where obliquity is the solution of a trigonometric equation, our approach allows us to identify not only the mean obliquity of satellites, but also their nutation in space as well as their polar motion (PM) with respect to the solid surface. We show that triaxiality has a significant effect on the mean obliquities of CS \(\textrm{I}\) (up to \(55\%\) for Titan), CS \(\textrm{II}\) and CS \(\textrm{IV}\) (up to 22 degrees for Callisto), but no effect on CS \(\textrm{III}\) . We assess the stability of the Cassini states over a wide range of free and forced precession frequency ratios and find that CS \(\textrm{I}\) and CS \(\textrm{III}\) are always stable. CS \(\textrm{II}\) and CS \(\textrm{IV}\) are only stable for relatively fast orbital precession and are therefore unstable for the current orbital parameters of the Galilean satellites and Titan. By solving the dynamic equations governing CS \(\textrm{I}\) and CS \(\textrm{III}\) for a fully rigid satellite without averaging the external torque over the mean anomaly, we find mean obliquities close to their time-constant classical counterparts. CS \(\textrm{I}\) nutations are at quasi semi-diurnal and diurnal periods, with amplitudes ranging from tens to thousands of mas, while the nutations for CS \(\textrm{III}\) are at quasi-diurnal, semi-diurnal or quarter-diurnal periods. We analytically describe one diurnal and one long-period CS \(\textrm{I}\) polar motion, the latter induced by the precession of the pericenter. A third, ter-diurnal PM also appears in CS \(\textrm{III}\) . CS \(\textrm{I}\) diurnal and long-period PMs are of the same order of magnitude for Io, Europa and Ganymede, while the long-period PM dominates the solution for Callisto and Titan. The analytical solution for the long-period PM obtained here by second-order developments correctly describes the Moon’s actual behavior and could help explain future observations of the rotation of large satellites. By extending the system of equations governing CS \(\textrm{I}\) , including gravitational and pressure couplings between misaligned layers, we predict the orientation of the spin axes of the outer shell, internal ocean and solid interior for an ocean-bearing body. We find five eigenmodes, including a free ocean nutation (FON) and an inner Chandler wobble (ICW). Far from resonance with an eigenmode, obliquity asymptotically tends toward values close to or slightly higher than that of a fully solid satellite. A large amplitude of nutation and polar motion can be obtained through a node precession/FON resonance, but also through resonance between the pericenter precession and the ICW. Depending on the thickness of the shell and the ocean, polar motion is dominated by its diurnal or long-period term, with amplitudes that can reach tens or even hundreds of meters.