We present redshift space two-point ($\xi$), three-point ($\zeta$) and reduced three-point (Q) correlations of Ly$\alpha$ absorbers (Voigt profile components having HI column density, $N_{\rm HI}>10^{13.5}$cm$^{-2}$) over three redshift bins spanning $1.7< z<3.5$ using high-resolution spectra of 292 quasars. We detect positive $\xi$ up to 8 $h^{-1}$ cMpc in all three redshift bins. The strongest detection of $\zeta =1.81\pm 0.59$ (with Q$=0.68\pm 0.23$), is in $z=1.7-2.3$ bin at $1-2h^{-1}$ cMpc. The measured $\xi$ and $\zeta$ values show an increasing trend with $N_{\rm HI}$, while Q remains relatively independent of $N_{\rm HI}$. We find $\xi$ and $\zeta$ to evolve strongly with redshift. Using simulations, we find that $\xi$ and $\zeta$ seen in real space may be strongly amplified by peculiar velocities in redshift space. Simulations suggest that while feedback, thermal and pressure smoothing effects influence the clustering of Ly$\alpha$ absorbers at small scales, i.e $<0.5h^{-1}$ cMpc, the HI photo-ionization rate ($\Gamma_{\rm HI}$) has a strong influence at all scales. The strong redshift evolution of $\xi$ and $\zeta$ (for a fixed $N_{\rm HI}$-cutoff) is driven by the redshift evolution of the relationship between $N_{\rm HI}$ and baryon overdensity. Our simulation using best-fitted $\Gamma_{\rm HI}(z)$ measurements produces consistent clustering signals with observations at $z\sim 2$ but under-predicts the clustering at higher redshifts. One possible remedy is to have higher values of $\Gamma_{\rm HI}$ at higher redshifts. Alternatively the discrepancy could be related to non-equilibrium and inhomogeneous conditions prevailing during HeII reionization not captured by our simulations.
Comment: 25 pages, 17 figures; Accepted in MNRAS; Comments welcomed