Large-eddy simulation of an underexpanded sonic jet injection into supersonic crossflows is performed to obtain insights into key physics of the jet mixing. A high-order compact differencing scheme with a recently developed localized artificial diffusivity scheme for discontinuity-capturing is used. Progressive mesh refinement study is conducted to quantify the broadband range of scales of turbulence that are resolved in the simulations. The simulations aim to reproduce the flow conditions reported in the experiments of Santiago and Dutton [Santiago, J. G., and Dutton, J. C, "Velocity Measurements of a Jet Injected into a Supersonic Crossflow," Journal of Propulsion and Power, Vol. 132, 1997, pp. 264-273] and elucidate the physics of the jet mixing. A detailed comparison with these data is shown. Statistics obtained by the large-eddy simulation with turbulent crossflow show good agreement with the experiment, and a series of mesh refinement studies shows reasonable grid convergence in the predicted mean and turbulent flow quantities. The present large-eddy simulation reproduces the large-scale dynamics of the flow and jet fluid entrainment into the boundary-layer separation regions upstream and downstream of the jet injection reported in previous experiments, but the richness of data provided by the large-eddy simulation allows a much deeper exploration of the flow physics. Key physics of the jet mixing in supersonic crossflows are highlighted by exploring the underlying unsteady phenomena. The effect of the approaching turbulent boundary layer on the jet mixing is investigated by comparing the results of jet injection into supersonic crossflows with turbulent and laminar crossflows.