Shor's algorithm - why doesn't the final collapse of the auxiliary qubits cripple the computation?
In the factoring algorithm, there are three kinds of qubits. In the OP's notation, there are "input qubits", which start in a superposition of all possible values, and which you eventually take the Fourier transform of. There are "value qubits", in which you compute the function $y^a \pmod{N}$, where $a$ is the value in the input qubits. And there are "auxiliary qubits", which you use as workspace to help do this computation.
In order to make the factoring algorithm work properly, you need to reset all the auxiliary qubits, which started as $|0\rangle$ at the beginning of the computation, to $|0\rangle$ at the end of the computation. This is called "uncomputing" these qubits. (Actually, you can set them to anything you please as long as it is a constant independent of the workings of the algorithm.) Theorems about reversible classical computation ensure that it is possible to do this.
If you reset the auxiliary qubits to $|0\rangle$, then if the environment, or somebody, measures them, nothing is revealed about the computation, and the computation is not "crippled". If you forget to reset them to $|0\rangle$, you probably won't get the right answer, whether or not anybody measures them.
Update: I originally thought the question was referring to the "value" qubits when the asker said "auxilliary". This answer explains why you don't need to measure the value qubits. For the actual auxilliary qubits, which are used as workspace while computing the value qubits, it should also be okay to measure them afterwards but only because a proper circuit uncomputes them back to 0.
After the value qubits are computed (the ones storing $B^k \text{ mod } R$, not the ones you used as helpers while computing that value!), you can just throw them out. You don't need to measure them or protect them or care for them. Just drop them on the floor. Nothing anyone does to them can hurt the remaining part of the computation. See this walkthrough of Shor's algorithm.
Let's do a simple example via my simulator Quirk. We will initialize a uniform superposition of qubits, and then compute their parity onto an auxiliary qubit (click the image to manipulate the circuit in the simulator):
The two green boxes are showing a representation of the density matrix of the top three qubits. We can show that information without disturbing the system because this is a simulator.
Before the parity computation, the qubits are entirely coherent. Afterwards, some of the off-diagonal indicators have disappeared (become zero). This indicates a partial loss of coherence. The states with an even number of ones have been decohered from the states that have an odd number of ones.
Now let's try to use the auxiliary qubit to "mess with" the top three qubits. If we succeed, the density matrix display will show something different. First thing to try is measurement:
Nothing's different.
Maybe we measured along the wrong axis? Let's rotate the qubit before measuring it:
Still no change!
In fact, no matter what we do to the bottom qubit, we can't change the density matrix of the top three qubits. Not without some kind of operation that crosses between them, or some kind of conditioning (e.g. only consider the subset of states where the measurement of the bottom qubit returned a particular result).
If you find this hard to believe, I recommend simply messing around in Quirk for a while, trying to make the densities of the top three qubits change by operating only on the bottom qubit.
Another way to confirm that it doesn't matter if you measure the auxilliary qubits is to just do the algebra and check.
The initial state is:
$$|\psi_0\rangle = |0\rangle_{\text{main}} \otimes |0\rangle_{\text{aux}} = |0\rangle_{\text{all}}$$
Then we Hadamard transform the main register:
$$|\psi_1\rangle = H_{\text{main}} |\psi_0\rangle = \sum_{k=0}^{2^n-1} |k\rangle_{\text{main}} \otimes |0\rangle_{\text{aux}}$$
Note that I am ignoring normalizing factors. In the end my argument is going to be based on the proportional size of various cases, instead of absolute size, so this is fine.
Then we pick a random base $B$ apply the modular exponentiation operation, which adds B^k mod R into the auxiliary register where k is the computational basis value of the main register. On an actual machine we'd use some temporary workspace to implement this operation, but that all gets cleaned up so here we only care about the effect on the main and aux registers:
$$M = \Big[ \text{aux} \text{ += } B^{\text{main}} \text{ mod } R \Big]$$
$$|\psi_2\rangle = M \cdot |\psi_1\rangle = \sum_{k=0}^{2^n-1} |k\rangle_{\text{main}} \otimes |B^k \text{ mod } R\rangle_{\text{aux}}$$
Now we can rewrite $k$ to be in terms of the unknown period of $B^k \text{ mod } R$. We'll use $k = l \cdot m + s$ where $l$ is the period, $s$ is an iteration variable for the offset between 0 and $l$, and $m$ is an iteration variable. With that in mind, we rewrite $|\psi_2\rangle$ as:
$$|\psi_2\rangle = \sum_{m=0}^{\;\;\lceil 2^n / l \rceil-1\;\;} \sum_{s=0}^{\text{min}(l, 2^n-lm)-1} |lm+s\rangle_{\text{main}} \otimes |B^{lm+s} \text{ mod } R\rangle_{\text{aux}}$$
Note that $B^{lm+s} = B^{s} \pmod{R}$. Also note that the complicated boundary conditions on $s$ can be simplified by approximating our actual sum with a sum that goes up to the first multiple of $l$ after $2^n$. This is a good approximation as long as $2^n >> l$, which is true since $n$ is chosen such that $2^n > R^2$ and we know that $R > l$. Anyways, after applying that simplification and approximation we get:
$$|\psi_2\rangle \approx \sum_{m=0}^{\;\;\lceil 2^n / l \rceil-1\;\;} \sum_{s=0}^{l-1} |lm+s\rangle_{\text{main}} \otimes |B^{s} \text{ mod } R\rangle_{\text{aux}}$$
Because the boundary condition of $s$ doesn't depend on $m$ anymore, we can rearrange the order of the sums. Which gives us:
$$|\psi_2\rangle \approx \sum_{s=0}^{l-1} \left( |B^{s} \text{ mod } R\rangle_{\text{aux}} \otimes \sum_{m=0}^{\lceil 2^n / l \rceil-1} |lm+s\rangle_{\text{main}} \right)$$
Now we apply the inverse Fourier transform operation to the main register. Notice that it can be moved from the outside of the sum to the inside:
$$\begin{align} |\psi_3\rangle &= \text{QFT}^{\dagger}_{\text{main}} \cdot |\psi_2\rangle \\ &\approx \text{QFT}^{\dagger}_{\text{main}} \cdot \sum_{s=0}^{l-1} \left( |B^{s} \text{ mod } R\rangle_{\text{aux}} \otimes \sum_{m=0}^{\lceil 2^n / l \rceil-1} |lm+s\rangle_{\text{main}} \right) \\ &= \sum_{s=0}^{l-1} \left( |B^{s} \text{ mod } R\rangle_{\text{aux}} \otimes \sum_{m=0}^{\lceil 2^n / l \rceil-1} \text{QFT}^{\dagger}_{\text{main}} \cdot |lm+s\rangle_{\text{main}} \right) \end{align}$$
Then expand the definition of the QFT into a sum over a variable $j$, and move that sum outward:
$$\begin{align} |\psi_3\rangle &\approx \sum_{s=0}^{l-1} \left( |B^{s} \text{ mod } R\rangle_{\text{aux}} \otimes \sum_{m=0}^{\lceil 2^n / l \rceil-1} \;\; \sum_{j=0}^{2^n-1} |j\rangle_{\text{main}} \cdot \text{exp}(i\tau \cdot 2^{-n} \cdot j \cdot (lm+s)) \right) \\ &= \sum_{s=0}^{l-1} \left( |B^{s} \text{ mod } R\rangle_{\text{aux}} \otimes \sum_{j=0}^{2^n-1} |j\rangle_{\text{main}} \sum_{m=0}^{\lceil 2^n / l \rceil-1} \text{exp}(i\tau \cdot 2^{-n} \cdot j \cdot (lm+s)) \right) \\ &= \sum_{s=0}^{l-1} \sum_{j=0}^{2^n-1} \left( |B^{s} \text{ mod } R\rangle_{\text{aux}} \otimes |j\rangle_{\text{main}} \sum_{m=0}^{\lceil 2^n / l \rceil-1} \text{exp}(i\tau \cdot 2^{-n} \cdot j \cdot (lm+s)) \right) \end{align}$$
Now we're going to measure the main register. The probability of getting the result $r$ is the total squared magnitude of states where the first register is $r$. Algebraically:
$$\begin{align} P(r) &= \sum_{a,b | a=r} \Big| (\langle a |_{\text{main}} \otimes \langle b |_{\text{aux}}) \cdot | \psi_3 \rangle \Big|^2 \\ &= \sum_{b} \Big| (\langle r |_{\text{main}} \otimes \langle b |_{\text{aux}}) \cdot | \psi_3 \rangle \Big|^2 \\ &\approx \sum_{b} \left| \langle r |_{\text{main}} \langle b |_{\text{aux}} \cdot \sum_{s=0}^{l-1} \sum_{j=0}^{2^n-1} |B^{s} \text{ mod } R\rangle_{\text{aux}} |j\rangle_{\text{main}} \sum_{m=0}^{\lceil 2^n / l \rceil-1} \text{exp}(i\tau \cdot 2^{-n} \cdot j \cdot (lm+s)) \right|^2 \end{align}$$
Because all of our basis kets are perpendicular, any summands that fail to satisfy $b=B^s \pmod{R}$ and $r=lm+s$ will be zero'd out. The remaining terms have bras and kets that exactly match, giving an inner product of 1. I'll do this over a few steps because it simplifies the sum considerably:
$$\begin{align} P(r) &\approx \sum_{b} \left| \langle r |_{\text{main}} \langle b |_{\text{aux}} \cdot \sum_{s=0}^{l-1} \sum_{j=0}^{2^n-1} |B^{s} \text{ mod } R\rangle_{\text{aux}} |j\rangle_{\text{main}} \sum_{m=0}^{\lceil 2^n / l \rceil-1} \text{exp}(i\tau \cdot 2^{-n} \cdot j \cdot (lm+s)) \right|^2 \\ &= \sum_{b} \left| \sum_{s=0}^{l-1} \sum_{j=0}^{2^n-1} \langle r |_{\text{main}} \langle b |_{\text{aux}} \cdot |B^{s} \text{ mod } R\rangle_{\text{aux}} |j\rangle_{\text{main}} \sum_{m=0}^{\lceil 2^n / l \rceil-1} \text{exp}(i\tau \cdot 2^{-n} \cdot j \cdot (lm+s)) \right|^2 \\ &= \sum_{b} \left| \sum_{s=0}^{l-1} \sum_{j=0}^{2^n-1} \langle r | j\rangle_{\text{main}} \langle b | B^{s} \text{ mod } R\rangle_{\text{aux}} \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot j \cdot (lm+s)) \right|^2 \\ &= \sum_{s=0}^{l-1} \left| \sum_{j=0}^{2^n-1} \langle r | j\rangle_{\text{main}} \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot j \cdot (lm+s)) \right|^2 \\ &= \sum_{s=0}^{l-1} \left| \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot r \cdot (lm+s)) \right|^2 \end{align}$$
Now we're getting somewhere. The next thing to do is get rid of that pesky $s$. Factor the $s$ component out of the inner sum, which lets you factor it out of the squared-magnitude, at which point you realize it contributes nothing and the sum can turn into a multiplication by $l$:
$$\begin{align} P(r) &\approx \sum_{s=0}^{l-1} \left| \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot r \cdot (lm+s)) \right|^2 \\ &= \sum_{s=0}^{l-1} \left| \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot r \cdot lm) \cdot \exp(i\tau \cdot 2^{-n} \cdot r \cdot s) \right|^2 \\ &= \sum_{s=0}^{l-1} \left| \exp(i\tau \cdot 2^{-n} \cdot r \cdot s) \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot r \cdot lm) \right|^2 \\ &= \sum_{s=0}^{l-1} \big| \exp(i\tau \cdot 2^{-n} \cdot r \cdot s) \big|^2 \left| \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot r \cdot lm) \right|^2 \\ &= \sum_{s=0}^{l-1} \left| \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot r \cdot lm) \right|^2 \\ &= l \cdot \left| \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot rl \cdot m) \right|^2 \end{align}$$
Nearly there. To make the structure of the sum blatant, we extract a variable $\omega = \exp(i\tau rl / 2^{n})$:
$$\begin{align} P(r) &\approx l \cdot \left| \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot rl \cdot m) \right|^2 \\ &= l \cdot \left| \sum_{m=0}^{\lceil 2^n / l \rceil-1} \omega^{m} \right|^2 \text{ where } \omega = \exp(i\tau r l / 2^{n}) \end{align}$$
The inner sum will be largest when all of its terms point in the same direction, i.e. when $\omega \approx 1$. Which means $\exp(i\tau rl / 2^{n}) \approx 1$, which in turn means that $2^{-n} r l$ is nearly an integer $d$. Rewrite $2^{-n} r l \approx d$ and you get:
$$r \approx 2^n \cdot \frac{d}{l}$$
In other words, if the period is $l$ then the values you are most likely to measure are placed near multiples of $2^n / l$. In practice you recover $l$ by solving "which possible multiple was my measurement nearest to?".
I leave it as an exercise for the reader to work out exactly how much more likely you are to measure values of $r$ that give values of $\omega$ close to 1.