The mass flow rate of Poiseuille flow of rarefied gas through long ducts of two-dimensional cross-sections with arbitrary shape is critical in the pore-network modeling of gas transport in porous media. Here, for the first time, the high-order hybridizable discontinuous Galerkin (HDG) method is used to find the steady-state solution of the linearized Bhatnagar–Gross–Krook equation on two-dimensional triangular meshes. The velocity distribution function and its traces are approximated in piecewise polynomial spaces (of degree up to 4) on the triangular meshes and mesh skeletons, respectively. By employing a numerical flux that is derived from the first-order upwind scheme and imposing its continuity weakly on the mesh skeletons, global systems for unknown traces are obtained with fewer coupled degrees of freedom when compared to the original discontinuous Galerkin formulation. To achieve fast convergence to the steady-state solution, a diffusion-like equation for flow velocity, which is asymptotic-preserving into the fluid dynamic limit, is solved by the HDG simultaneously on the same meshes. The proposed HDG-synthetic iterative scheme is proved to be accurate and efficient. Specifically, for flows in the near-continuum regime, numerical simulations have shown that, to achieve the same level of accuracy, our scheme could be faster than the conventional iterative scheme by two orders of magnitude, also it is faster than the synthetic iterative scheme based on the finite difference discretization in the spatial space by one order of magnitude. In addition, the implicit HDG method is more efficient than an explicit discontinuous Galerkin gas kinetic solver, as well as the implicit discontinuous Galerkin scheme when the degree of approximating polynomial is larger than 2. The HDG-synthetic iterative scheme is ready to be extended to simulate rarefied gas mixtures and the Boltzmann collision operator.