In the paper, a system of hyperbolic linear Bhatnagar–Gross–Krook equations with single relaxation time for simulation of linear diffusion is presented. The lattices based on orthogonal, antitropic and null velocities are considered. Cases of parametric and non-parametric coefficients are considered. The Chapman–Enskog method is applied to the derivation of the linear diffusion equation from the proposed system. The expression for the diffusion coefficient with dependency on relaxation time is obtained. The structure of the solution of the proposed system is analyzed. Lattice Boltzmann equations and finite-difference-based lattice Boltzmann schemes are considered as discrete approximations for the proposed system. The existence of the fictitious numerical diffusion is demonstrated. Stability with respect to initial conditions is analyzed by investigation of wave modes. The necessary stability condition is written as a physical consistent condition of positivity of relaxation time. The sufficiency of this condition is demonstrated analytically and numerically for all considered lattices in parametric and nonparametric cases. It is demonstrated, that the dispersion of the solutions takes place at some values of relaxation time and parameter. The dispersion exists due to hyperbolicity of the proposed system. The presence of the dispersion may lead to effects, which are not typical for the solution of linear diffusion equation. But it is demonstrated, that the effect of the dispersion may take place only in the cases, when image parts of the frequencies of wave modes are equal to null. In the case of positive image parts this effect is damped and solutions decrease at the same manner, as a solution of the diffusion equation. Presented system may be considered as a basis for the construction of lattice Boltzmann equations and lattice Boltzmann schemes of various accuracy orders.