A computationally efficient methodology for propagating the epistemic uncertainty in the median ground motion in probabilistic seismic hazard analysis is developed using the polynomial chaos (PC) approach. For this application, the epistemic uncertainty in the median ground motion for a specific scenario is assumed to be lognormally distributed and fully correlated across earthquake scenarios. In the hazard calculation, a single central ground‐motion model (GMM) is used for the median along with the epistemic standard error of the median for each scenario. A set of PC coefficients is computed for each scenario and each test ground‐motion level. The additional computation burden in computing these PC coefficients depends on the order of the approximation but is less than computing the median ground motion from one additional GMM. With the PC method, the mean and fractiles of the hazard due to the epistemic uncertainty distribution of the median ground motion are computed as a postprocess that is very fast computationally. For typical values of the standard deviation of epistemic uncertainty in the median ground motion ( natural log units), the methodology accurately estimates the epistemic uncertainty distribution of the hazard over the 1%–99% range. This full epistemic range is not well modeled with just a small number of GMM branches uses in the traditional logic‐tree approach. The PC method provides more accuracy, faster computation, and reduced memory requirements than the traditional approach. For large values of the epistemic uncertainty in the median ground motion, a higher order of the PC expansion may be needed to be included to capture the full range of the epistemic uncertainty.