Significance of "virtual" zeroes/poles introduced by the Pade expansion of $e^{sT}$ Let $s$ be a complex parameter and $T$ a real. The function $e^{sT}$ is entire. It can also be expanded in a Pade approximant via
$$
e^{sT} = \frac{e^{sT/2}}{e^{-sT/2}} = \frac{1+ sT/2 + \frac{(sT/2)^2}{2!}+ \cdots}{1 - sT/2 + \frac{(sT/2)^2}{2!} - \cdots},
$$
which upon truncation to first order yields the bilinear mapping
$$
\frac{2+sT}{2-sT}.
$$
Approximations of the exponential function by bilinear mappings are found frequently in digital signal processing, linearized analysis of delay systems, and other applications in systems theory, where they are known as "Tustin's method" or sometimes "Tustin's approximation". 
In all of these applications, poles and zeros introduced into a system transfer function are important factors which determine system stability, transient characteristics, and frequency response. The character of a system, say $e^{sT}H(s)$, where $H(s)$ is an arbitrary rational transfer function, may in fact be changed considerably if the exponential factor is approximated by Tustin's method, despite the fact that numerically the two are close when the latter is defined and nonzero, since the introduction of the pole-zero combination affects relative stability margins and the shape of the system's frequency response. 
Clearly the pole-zero combination is "virtual" in the sense that they do not arise from physical considerations, but it is unclear to me what their exact significance is. 
 A: In my comment I mistakenly assumed that your system is like
$$\dot{x}(t) = a x(t-T) + u(t)$$
which has an infinite number of poles that can be described with Lambert-W function. But instead your system is like
$$\dot{x}(t) = a x(t) + b u(t-T)$$
which becomes like 
$$\dot{x}(t) = a x(t) + b k x(t-T)$$
for the closed loop system, which also has an infinite number of poles.
But you are right, for the open loop system, there cannot be additional zeros/poles when there is input or output delay. Assuming a causal transfer function (hence $e^{-sT}$ multiplier) the delay effect can be obtained with a non-minimum phase system, because it oscillates around zero (depending on the order) before it fades out and non-delay part takes over. See the figure in page 31 of http://engineering.nyu.edu/mechatronics/Control_Lab/Criag/Craig_RPI/2002/Week2/First-Order_Process_Time_Delay_2002.pdf to see what I mean.
A: This is not an answer but it is too long for a comment.
A Padé approximant is not the ratio of Taylor series. Taking the case of $e^x$, what you wrote is
$$e^x=\frac {1+\frac x 2 +\frac {x^2} 8 } {1-\frac x 2 +\frac {x^2} 8 } \tag 1$$ while the $[2,2]$ Padé approximant is
$$e^x=\frac {1+\frac x 2 +\frac {x^2} {12} } {1-\frac x 2 +\frac {x^2} {12} } \tag 2$$  Developed as series
$$e^x-\frac {1+\frac x 2 +\frac {x^2} 8 } {1-\frac x 2 +\frac {x^2} 8 } =\frac {x^\color{red}{3}} {24}$$
$$e^x-\frac {1+\frac x 2 +\frac {x^2} {12} } {1-\frac x 2 +\frac {x^2} {12} }= \frac {x^\color{red}{5}} {720}$$
Consider also the norms
$$\Phi_1=\int_{-1}^{+1} \Bigg[e^x-\frac {1+\frac x 2 +\frac {x^2} 8 } {1-\frac x 2 +\frac {x^2} 8 } \Bigg]^2\, dx=1.625 \times 10^{\color{red}{-3}}$$
$$\Phi_2=\int_{-1}^{+1} \Bigg[e^x-\frac {1+\frac x 2 +\frac {x^2} {12} } {1-\frac x 2 +\frac {x^2} {12} }\Bigg]^2\, dx=1.255 \times 10^{\color{red}{-6}}$$
