Mathematics of the Virtual Displacement So I'm pretty certain this question has been asked to death here, but I still can't find a good explanation of a very particular aspect of the virtual displacements in physics. 
Background 
For clarity, let me list (vaguely) what I do understand, and am not requesting an explanation for 


*

*"What" are virtual displacements? 

*"What" does  D'Alambert's principle mean? 

*Why are virtual displacements necessary in the derivation of the Lagrangian? 

*Why are virtual displacements perpendicular to constraint forces? 


The real question I guess would be 


*

*How are virtual displacements perpendicular to constraint forces? 


Let me be more specific. As a note, there are various references for this, Goldstein for example, but in the interest of being concise I'm going to particularly reference Fetter and Walecka's Continuum book, page 51. 
Setup
We have some $f_j (x_1,\dots,x_n,t) = c_j $ for $j = 1,\dots,k$ constraint equations, where $k$ reflects the number of (holonomic?) constraints to which the system is subject. These $k$ constraint equations may be either (a) explicitly time dependent (rheonomous) or (b) implicitly time dependent (scleronomous). 
From these $k$ equations, we may identify $n-k$ generalized coordinates, $q_\sigma$, with which we can parameterize the positions of the particles in the system according to 
$$x_i = x_i(q_1,q_2,\dots,q_{n-k},t)$$
Problem 
Now, and this makes sense, any infinitesimal displacement in the system yields a change in any $x_i$ coordinate 
$$dx_i = \sum_{\sigma=1}^{n-k} \frac{\partial x_i}{\partial q_\sigma}dq_\sigma + \frac{\partial x_i}{\partial t}dt \tag{eq. 1}$$
Now, what makes a lot less sense is why a variation of the variable $x_i$ yields 
$$\delta x_i = \sum_{\sigma=1}^{n-k} \frac{\partial x_i}{\partial q_\sigma}\delta q_\sigma \tag{eq. 2}$$
Now, unless it comes with a (reasonably) rigorous mathematical analysis, I will not accept the following answer


*

*Variations are instantaneous, i.e. occur in zero time. 

*In other words, $dt = 0$. 


Qualitatively, sure, that answer is fine. But quantitatively, I don't see how it's well grounded within the framework of calculus. In essence, it's too "handwavy" an argument for me. 
What I am really interested in is how the requirement that constraint forces and virtual displacements must be perpendicular leads to the definition (eq.2) of virtual displacements. Accordingly, answers to any of the following are what I am looking for


*

*Quantitatively, what steps must you take to go from (eq.1) to (eq.2)?  

*Why, if time were not standing still, would a virtual displacement not be perpendicular to a constraint force?   

*Virtual displacements and infinitesimals are identical in scleronomous systems. In a rheonomous system, how does the definition of a virtual displacement ensure that it is always perpendicular to a time dependent constraint force? How is this different from the way an infinitesimal displacement might behave? 

 A: Here is one way to imagine a variation. You have a path $\vec y(t)=(y_1(t), ..., y_n(t)),$ where $y_i(t) = x_i(q_1(t),q_2(t),\dots,q_{n-k}(t),t).$ So the idea is that for every time $t$ you have some $q$s and they (together with the $t$) give you all your $y_i$ hence give you your $\vec y.$  Seems basic but the idea is that given a $t$ you have an $\vec y.$
Now for a variation you can consider a second path which we will call $\vec Y (t)$ but we also want it to satisfy the constraints so we will consider some functions $Q_j(t)$ so that we can get $\vec Y(t)=(Y_1(t), ..., Y_n(t)),$ where $Y_i(t) = x_i(Q_1(t),Q_2(t),\dots,Q_{n-k}(t),t).$ And that notation might seem confusing but think about how you really have functions $q_j(t)$ and functions $Q_j(t)$ and they combine with the time $t$ to give you the real coordinates in the unconstrained space.
Now we can talk about a variation. And the variation is all about a whole family of paths. instead of using $q_j(t)$ or $Q_j(t)$ you can consider the function $q^\delta_j(t)=(1-\delta)q_j(t)+\delta Q_j(t)$ and then we get a path for every value of $\delta$ and in fact the corresponding path is $\vec y^\delta(t)=(y^\delta_1(t), ..., y^\delta_n(t)),$ where $y^\delta_i(t) = x_i(q^\delta_1(t),q^\delta_2(t),\dots,q^\delta_{n-k}(t),t).$ 
Basically for each value of $\delta$ you have interpolated the $q$s to be between the original path $q^0=q$ and the other path $q^1=Q.$ And the only real difference is that we transform from the space of $q$ and $t$ to the space of $x$s.
So now when you do a variation you imagine the paths corresponding to different values of $\delta,$ subtract the two paths and divide by the difference in $\delta.$ In fact you can look at $\delta \vec y(t)$ as a short hand for $\frac{d}{d\delta}\vec y^\delta(t).$ You are asking about how the coordinates in the final n dimensional space vary as you change $\delta.$
Since you are really asking how the whole path changes you are doing this for each time. So you can think of $\delta y_i$ as $\left.\frac{d}{d\delta}g_i(\delta)\right|_{\delta=0}$ where $g_i(\delta)=y^\delta_i(t).$ And $\delta q_j$ as $\left.\frac{d}{d\delta}h_j(\delta)\right|_{\delta=0}$ where $h_j(\delta)=q^\delta_j(t).$
So you are finding something new for each $t$ not $q(t)$ not $Q(t)$ and not even $q^\delta$ but finding $\delta y_i=\left.\frac{d}{d\delta}g_i(\delta)\right|_{\delta=0}$ in terms of $\delta q_j=\left.\frac{d}{d\delta}h_j(\delta)\right|_{\delta=0}.$ Which itself is something you find out for each time. What it tells you is how the path was changing at that time and the direction and the how "fast" in the change is in the variation space $\delta.$
The idea is that by taking the derivative with respect to $\delta$ at $\delta=0$ then you as deal with path close to the original (since you evaluate at $\delta=0$) and you are seeing how they vary (change) i.e. $d/d\delta.$
