What exactly is a virtual displacement in classical mechanics? I'm reading Goldstein's Classical Mechanics and he says the following:

A virtual (infinitesimal) displacement of a system refers to a change in the configuration of the system as the result of any arbitrary infinitesimal change of the coordinates $\delta \mathbf{r}_i$, consistent with the forces and constraints imposed on the system at the given instant $t$. The displacement is called virtual to distinguish it from an actual displacement of the system ocurring in a time interval $dt$, during which the forces and constraints may be changing.

Then he discusses virtual work and so on. Now, I can't grasp what this thing of virtual really is. By this text there's a diference between one infinitesimal change and one virtual change and I really don't get what this virtual really is.
Also, this is based on infinitesimals. How can this be expressed rigourously without refering to infinitesimals? I tried looking on Spivak's Physics for Mathematicians where he considers these virtual displacements as tangent vectors to a certain manifold, but I'm not sure this is the most "standard" way to do it rigorously.
 A: In short: virtual displacement is "pretend you are moving, but don't really move". In other words - you move by such a small amount that you don't change the state of the system - but it gives you insight (through work done etc) in what would happen if you did move.
In other words - if the system is really moving, you can look at an interval $dt$ to see how little it moved in that time. That's an "infinitesimal" motion. With virtual motion, you pretend you moved by $dx$ - but not because the system is in motion, but just imagining that you made the tiniest motion (in finite time - so there is no velocity, $\frac{dx}{dt}=0$)
A: Let $Q$ denote the set of all possible configurations of the system (the configuration manifold).  Consider a point $q_0\in Q$.  For the sake of conceptual clarity, and to make contact with physics notation, let's work in some local coordinate patch around $q_0$.  
Suppose that $q_0$ represents the position of the system under consideration at time $t_0$.  At a given time $t$ later, the system will be at some position say $q(t)$ that is determined by the evolution equations (the Euler-Lagrange equations if we are doing Lagrangian mechanics), and the quantity 
\begin{align}
q(t) - q(t_0) = q(t) - q_0
\end{align}
would be the displacement of the system after a time $t-t_0$.  Suppose, instead we consider some other curve $\gamma(s)$ in the configuration space which starts at the point $s_0$;
\begin{align}
    \gamma(s_0) = q_0,
\end{align}
and suppose that we compute the displacement
\begin{align}
    \gamma(s) - \gamma(s_0) = \gamma(s) - q_0
\end{align}
that would result from moving along this other curve of our choosing.  We call this displacement the virtual displacement after a "time" $s-s_0$ corresponding to moving along the curve $\gamma$.  It's called virtual because it is the displacement in the position of the system that would occur if the system were to move along the curve $\gamma$ of our choosing  -- a "virtual" curve as opposed to the "real" curve along which the system travels according to the Lagrangian evolution of the system.
Note. As Qmechanic suggested in the comments, I used the parameter $s$ for the virtual curve $\gamma$ instead of $t$ to emphasize that moving along that curve does not correspond to time-evolution, but rather any curve of our choosing.
Now what about virtual "infinitesimal" displacements?  Well, recall that the term "infinitesimal" in physics essentially always refers to "first order" approximations, see, e.g. this SE post:
Rigorous underpinnings of infinitesimals in physics
So when we are discussing a virtual infinitesimal displacement, what we have in mind is taking the virtual displacement $\gamma(s) - q_0$, Taylor expanding it to first order in $s$, and extracting only the first order term.  Let's do this:
\begin{align}
    \gamma(s) - q_0 = \gamma(s_0) + \dot\gamma(s_0) (s-s_0) + O((s-s_0)^2) - q_0
\end{align}
Using the fact that $\gamma(s_0) = q_0$, we see that the Taylor expansion of the virtual displacement is
\begin{align}
    \gamma(s) - q_0 = \dot\gamma (s_0) (s-s_0) + O((s-s_0)^2),
\end{align}
and now we notice that to first order in $s$, the size of the virtual displacement is controlled by the coefficient of $s-s_0$, namely $\dot\gamma(s_0)$.  In other words, virtual infinitesimal displacements (meaning we just keep the first order contribution in $s-s_0$), are determined by the velocity vector of the chosen "virtual curve" at $s_0$.   But if you've taken a differential geometry course, then you know that velocities of curves on a manifold are simply tangent vectors to that manifold!
So virtual infinitesimal displacements can be associated with tangent vectors to the configuration manifold.  The intuition to keep in mind here as that a virtual displacement just tells us how far we would get away from a certain point on the manifold if we were to travel on a certain curve of our choosing that may not coincide with the actual motion of the system determined by time evolution.  The "infinitesimal" part and identifying this part with tangent vectors comes simply from considering what happens only to first order.
A: You may look into this. It has elaborate works on understanding what virtual displacement is. 
 https://www.researchgate.net/publication/2174249_On_Virtual_Displacement_and_Virtual_Work_in_Lagrangian_Dynamics 
A: As I understand, it must be a displacement in generalized coordinates. If they are orthogonal space coordinates they are not virtual.
