Abstract Finite heat flux often exists in space and astrophysical plasmas, which can be a free energy source for heat flux instability. The solar wind is a well-known example of such plasmas and a number of previous studies have investigated the characteristics of heat flux instability in the context of solar wind. In the literature there exists some uncertainties regarding the properties of heat flux instability. While some linear theories predict the association of the heat flux instability with right-hand polarized whistler waves, other studies argue for left-hand polarized unstable modes. The present study investigates the nonlinear development of initially unstable left-hand heat flux mode by means of particle-in-cell simulation. It is found that while the early phase is characterized by the left-hand polarization, in agreement with linear theory, as the wave amplitude becomes high and the instability enters the nonlinear phase, the dominant wave mode gradually switches over to the right-hand polarized waves. Such a behavior is related to the pitch angle scattering of the heat flux carrying electrons by nonlinear interaction with large-amplitude waves. The present study shows that the heat flux instability generally requires nonlinear treatment such that characterizing its behavior with linear theories may not always be adequate.