There is an ongoing debate concerning the best rheological model for liquid flows in elastohydrodynamic lubrication (EHL). Due to the small contact area and high relative velocities of bounding solids, the lubricant experiences pressures in excess of 500 MPa and strain rates that are typically 10 5 − 10 7 s −1 . The high pressures lead to a dramatic rise in Newtonian viscosity ηN and the high rates lead to large shear stresses and pronounced shear-thinning. This paper presents detailed simulations of a model EHL fluid, squalane, using nonequilibrium molecular dynamics methods to extract the scaling of its viscosity with shear rate (10 5 − 10 10 s −1 ) over a wide range of pressure P (0.1 MPa to 1.2 GPa), and temperature T (150 − 373 K). Simulation results are consistent with a broad range of equilibrium and nonequilibrium experiments. At high T and low P , where ηN is low, the response can be fit to a power-law, as in the common Carreau model. Shear-thinning becomes steeper as ηN increases, and for ηN 1 Pa-s, shear-thinning is consistent with the thermally activated flow assumed by another common model, Eyring theory. Simulations for a bi-disperse Lennard-Jones (LJ) system show that the transition from Carreau to Eyring is generic. For both squalane and the LJ system, the viscosity decreases by only about a decade in the Carreau regime, but may fall by many orders of magnitude in the Eyring regime. Shear thinning is often assumed to reflect changing molecular alignment, but the alignment of squalane molecules saturates after the viscosity has dropped by only about a factor of three. In contrast, thermal activation describes shear thinning by six or more decades in viscosity. Changes in the diagonal elements of the stress tensor with rate and shear stress are also studied. * vjadhao@iu.edu † mr@jhu.edu [6,7]. A second reason is that EHL conditions present fundamental theoretical challenges. They are in the regime where fluids exhibit poorly understood glassy rheology, with rapid rises in η N with pressure, strongly nonlinear viscoelastic behavior, rapid shear thinning, and changes in molecular order. While there is a general agreement on the qualitative changes in shear stress σ, the underlying mechanisms of shear thinning and the functional form of the rheological response are hotly debated [4,6,7].Two phenomenological models at the center of this active debate are based on work by Eyring [8,9] and Carreau [10], respectively. The Eyring model assumes that flow occurs by thermal activation over a single characteristic energy barrier height and predicts that σ rises as log(γ) at highγ. In the Carreau model, σ increases as a power law, σ ∝γ n , with n typically around 0.5. As discussed below, this power law behavior can arise from several different physical pictures of flow, including thermal activation over a broad range of energy barrier heights.Both of these rheological models can be fit to existing rheometer data on simple fluids that extends up to rates of about 10 4 s −1 . The two models predict dramatically d...