The inter-arrival time of a Poisson Process, $t$, conforms to the exponential distribution, so the probability density function for $t$ is $f(t)=\lambda e^{-\lambda t}$, $t>0$. ($\lambda$ is the arrival rate of the Poisson Process.)
For two adjacent arrivals $i$ and $j$ ($i$ is in front of $j$), if their inter-request time $t$ is less than a constant $\tau$ ($t<\tau$), then remove $j$. If the arrival immediately after $j$ is $k$, and the inter-request time between $j$ and $k$ is still less than $\tau$, we also remove $k$, so on and so forth.
In other words, arrivals close enough ($t<\tau$) are all removed except the beginning one.
Here is my question: What is the process that the rest arrivals conform to? If there is no well-known process that can describe this process, what is the best way to approximate it?
Edit: Except for the answer from @Did, I found another answer regarding this question.
Let $t_i$ denote the inter-arrival time of the original Poisson Process, and $t$ denote the inter-arrival time for the remaining process. Assume $t_1$, $t_2$,...,$t_n$ to be a series of continuous inter-arrival variables for the original Poisson Process, with $t_1$ starting from arrival $A_j$ and $t_n$ ending with arrival $A_k$.
Further assume that arrivals $A_j$ and $A_k$ to be two adjacent arrivals that will remain (i.e., no arrivals remain between $A_j$ and $A_k$). Then we have
$$t=\sum_{i=1}^{n-1}t_i+t_n$$
And $t_i<\tau~(i=1...n-1)$, $t_n>\tau$.
So we can derive the distribution density function $g(t)$ for $t$ is
$$g(t)=\sum_{n=1}^{\infty}f(t|n)P\{t<\tau\}^{n-1}(1-P\{t<\tau\})$$
where $f(t)=\lambda e^{-\lambda t}$ and $P\{t<\tau\}=1-e^{-\lambda \tau}$.
(Note:The above formula seems to be right at first glance, but I do not know how to interpret mathematically.)
Based on this formula, the explicit form for $g(t)$ is
$$g(t)=\sum_{n=1}^{\infty}(-1)^{n+1}\lambda^{-1}\frac{(t-n\tau)^{n-1}}{(n-1)!}u(t-n\tau)$$ where $u(t)$ is a step function with $u(t)=0$ when $t<0$ and $u(t)=1$ otherwise.
The first remaining arrival time is the first arrival time $T_1$ of the original process $(T_n)_{n\geqslant1}$ hence it is exponentially distributed with parameter $\lambda$. The rest of the remaining process is a homogenous point process whose interarrival time distribution has the density $$f_{\lambda,\tau}:t\mapsto \lambda g(t)\exp\left(-\lambda\int_0^tg(s)\mathrm ds\right),$$ where the function $g$ is defined, for every $t\geqslant0$, by $$g(t)=1-P(A_t),$$ and $A_t$ is the event that $t\lt\tau$ or that, in the original process $(T_n)$, there exists some integer $n$ such that $T_1\lt\tau$, $T_2-T_1\lt\tau$, ..., $T_n-T_{n-1}\lt \tau$ and $t-T_n\lt\tau$. We now characterize $P(A_t)$.
If $t\lt\tau$, $P(A_t)=1$. If $t\gt\tau$, conditioning on the time $t-s$ of the last event before $t$ in the original process, one gets $$P(A_t)=\int_0^\tau P(A_{t-s})\lambda\mathrm e^{-\lambda s}\mathrm ds,$$ that is, $$P(A_t)=\lambda\mathrm e^{-\lambda t}\int_{t-\tau}^t P(A_s)\mathrm e^{\lambda s}\mathrm ds.$$ Thus, $P(A_t)=1$ for every $t\lt\tau$, and, for every $t\gt\tau$, $$\frac{\mathrm d}{\mathrm dt}P(A_t)=-\lambda\mathrm e^{-\lambda\tau}P(A_{t-\tau}).$$ For every $k\geqslant0$ and $u$ in $(0,1)$, let $$a_k(u)=P(A_{(k+u)\tau}),$$ then $a_0(u)=1$ for every $u$ in $(0,1)$, and, for every $k\geqslant1$ and $u$ in $(0,1)$, $$a'_k(u)=-ca_{k-1}(u),\qquad c=\lambda\tau\mathrm e^{-\lambda\tau},$$ hence each function $a_k$ is (the restriction to the interval $(0,1)$ of) a polynomial function of degree $k$ whose coefficients depend on ($k$ and) $\lambda\tau$ only, each function $a_k$ is decreasing on $(0,1)$, and $a_k\to0$ when $k\to\infty$.
Edit: A heuristics for the behaviour of $P(A_t)$ when $t$ is large is to assume that the shape of $a_k$ converges when $k\to\infty$ in the sense that $$\lim_{k\to\infty}\frac{a_k(u)}{a_k(0)}=b(u).$$ This suggests that $a'_k(u)\approx a_k(0)b'(u)$, hence the differential equation linking $a_k$ to $a_{k-1}$ would read $a_k(0)b'(u)\approx-ca_{k-1}(0)b(u)$. Since $a_k(0)=a_{k-1}(1)$, this suggests that $a_{k}(0)/a_{k-1}(0)\to\alpha$ and $b(u)\approx\mathrm e^{-cu/\alpha}$, where $\alpha$ solves the identity $$\alpha=\mathrm e^{-c/\alpha}.$$ Since one solution is $\alpha=\mathrm e^{-\lambda\tau}$, after some algebraic simplifications, all this suggests that, at least when $t$ is large, $$P(A_t)\approx\mathrm e^{-\lambda t}.$$ Edit-edit: Explicit formulas (but with summations) are as follows. Consider the generating function $$A(u,x)=\sum_{k=0}^\infty a_k(u)x^k,$$ then the identities $a_0(u)=1$ and $a_k(0)=a_{k-1}(1)$ for every $k\geqslant1$ yield $$A(0,x)=1+xA(1,x),$$ while the differential equations $a'_0(u)=0$ and $a'_k(u)=-ca_{k-1}(u)$ for every $k\geqslant1$ yield $$\partial_uA(u,x)=-cxA(u,x).$$ Thus, for every $u$ in $(0,1)$, $$A(u,x)=A(0,x)\mathrm e^{-cxu}.$$ Using this for $u=1$ and comparing to the formula for $A(1,x)$, one gets $$A(0,x)=\frac1{1-x\mathrm e^{-cx}},$$ and, finally, $$A(u,x)=\frac{\mathrm e^{-cxu}}{1-x\mathrm e^{-cx}}=\sum_{n\geqslant0}x^n\mathrm e^{-cx(u+n)}.$$ Expanding each exponential on the RHS and collecting the $x^k$ terms yields, for every $k\geqslant0$, $$a_k(u)=\sum_{i=0}^k(-1)^i\frac{c^i}{i!}(u+k-i)^i,$$ that is, for every $t\geqslant0$, $$P(A_t)=\sum_{n=0}^{\lfloor t/\tau\rfloor}\frac{(-1)^n}{n!}\mathrm e^{-n\lambda\tau}(\lambda t-n\lambda\tau)^n=\sum_{n=0}^\infty\frac{(-1)^n}{n!}\mathrm e^{-n\lambda\tau}[(\lambda t-n\lambda\tau)_+]^n.$$