Impact of inhomogeneous reionization on post-reionization 21 cm intensity mapping measurement of cosmological parameters

21 cm intensity mapping (IM) has the potential to be a strong and unique probe of cosmology from redshift of order unity to redshift potentially as high as 30. For post-reionization 21 cm observations, the signal is modulated by the thermal and dynamical reaction of gas in the galaxies to the passage of ionization fronts during the Epoch of Reionization. In this work, we investigate the impact of inhomogeneous reionization on the post-reionization 21 cm power spectrum and the induced shifts of cosmological parameters at redshifts $3.5 \lesssim z \lesssim 5.5$. We make use of hydrodynamics simulations that could resolve small-scale baryonic structure evolution to quantify HI abundance fluctuation, while semi-numerical large box 21cmFAST simulations capable of displaying inhomogeneous reionization process are deployed to track the inhomogeneous evolution of reionization bubbles. We discussed the prospects of capturing this effect in two post-reionization 21 cm intensity mapping experiments: SKA1-LOW and PUMA. We find the inhomogeneous reionization effect could impact the HI power spectrum up to tens of percent level and shift cosmological parameters estimation from sub-percent to tens percent in the observation of future post-reionization 21 cm intensity mapping experiments such as PUMA, while SKA1-LOW is likely to miss this effect at the redshifts of interest given the considered configuration. In particular, the shift is up to 0.0206 in the spectral index $n_s$ and 0.0192 eV in the sum of the neutrino masses $\sum m_\nu$ depending on the reionization model and the observational parameters. We discuss strategies to mitigate and separate these biases.