Chemical potential in canonical partition function I'm a bit confused on the interpretation of the chemical potential in a canonical ensemble (a system which can only exchange energy with a reservoir but not particles). 
Here is what I think I know:
As far as I understand, when one deals with a system that can exchange energy and particles with a reservoir one is dealing with a grand canonical ensemble. At the end of the day this means that we are either working with Boltzmann's factors, $e^{-\epsilon \beta}$ or with Gibb's factors, $e^{-(\epsilon -\mu)\beta}$, depending on whether or not the system exchanges particles with the reservoir. A way to see the Boltzmann's factor is as the special case of $\mu=0$ of the Gibb's factor.
One way of deriving these distributions (the one I've seen) is to start with something like this:
$$
\frac{P(s_2)}{P(s_1)} = \frac{\Omega_R(s_2)}{\Omega_R(s_1)}= \frac{e^{S_R(s_2)/k_b}}{e^{S_R(s_1)/k_b}}= e^{(S_R(s_2)-S_R(s_1))/k_b}
$$
now one invoques the thermodynamic identity $dU=TdS-PdV+\mu dN$, solves for $dS$ and discard the constant terms like $dV=0$ or $dN=0$ for the canonical case and substitute above. For the grand canonical case $dN\neq 0$ and one yields the Gibb's factor.
Here is my question:
I've seen that one can define the chemical potential through the canonical free energy $F=-k_bT \log (Z)$ by 
$$
\mu = \left( \frac{\partial F}{\partial N}\right)_{T,V}.
$$
What is the meaning of this $\mu$? We were dealing with a system that couldn't exchange particles with the reservoir and yet it has a non zero chemical potential like a Grand canonical ensemble!
 A: Let's consider a more general setting, because this is something very general (the answer for your specific setting is then given at the bottom). Suppose that one of the thermodynamic variables associated to your system is some extensive quantity $X$.
There are 2 ensembles that you can consider (I use $V$ for the volume; we need it, or another extensive quantity, to take the thermodynamic limit and define densities):


*

*In the first one, $X$ is fixed to some value $x_0 V$ ($x_0$ is thus the density of $X$ per unit volume). The partition function is $Z_V(x_0)$ and the corresponding thermodynamic potential is $F(x_0) = (1/V)\log Z_V(x_0)$.

*In the second one, $X$ is allowed to vary. The partition function is $Q_V(\lambda) = \sum_x e^{-\lambda x V} Z_V(x)$. Here, $\lambda$ is the intensive thermodynamic variable conjugate to $X$. Let us denote by $G(\lambda)=(1/V)\log Q_V(\lambda)$ the corresponding thermodynamic potential.


Then, in large systems, the sum over $x$ in the partition function $Q_V(\lambda)$ is usually dominated by just a few terms (justifying this requires steepest descent or basic large deviations theory and relies on the convexity properties of $F$):
$$
Q_V(\lambda) = \sum_x e^{-V\lambda x} Z_V(x) = \sum_x e^{V(-\lambda x+F(x))} \approx e^{V \max_x (-\lambda x + F(x))}.
$$
Now the maximum can be found by differentiating (assuming smoothness and strict convexity, that is, no phase transition): we get that the unique value $x$  that realizes the maximum is such that
\begin{equation}
\lambda = \frac{\partial F}{\partial x}.\tag{$\star$}
\end{equation}
We immediately deduce the standard relation (Legendre transform) between the thermodynamic potentials $F$ and $G$:
$$
G(\lambda) = \max_x (-\lambda x + F(x)) = -\lambda x_0 + F(x_0),
$$
provided $\lambda=\lambda(x_0)$ is chosen as the solution to ($\star$). This is the equivalence of ensembles (at the level of thermodynamic potentials). This shows you that your "$X$ can vary" ensemble leads to the same thermodynamic behavior as the "$X=x_0$" ensemble, provided
that you choose $\lambda$ in such a way that the maximum above is reached at $x_0$.

The point I want to make is that this works for each pair of conjugate thermodynamic quantities:


*

*chemical potential $\mu$ vs particle numbers $N$;

*inverse temperature $\beta$ vs internal energy $E$;

*magnetic field $h$ vs magnetization $M$, etc.


Of course, for each such pair, there are conventions for the definition of the corresponding thermodynamic potentials that differ from what I used (additional signs, maybe some prefactor $kT$, etc.). These have only been introduced to match the original definitions in thermodynamics.
The latter predating statistical mechanics, several choices of definitions are rather unfortunate and make notations more intricate than they should be (for instance, life would be easier if one was using $-\beta$ rather than $T$ for measuring the temperature: many formulas would look simpler, it would be clear why one cannot reach $0$ temperature, why "negative temperatures" are hotter than infinite ones, etc.).
In any case, apart form these details, the structure of the partition functions/thermodynamic potentials is the same as I used above.
For instance, using your notations:
\begin{align}
Z_{\rm g.c.}(\mu)
&= \sum_{s} e^{-\beta(E(s)-\mu N(s))} \\
&= \sum_N e^{\beta\mu N} \sum_{s: N(s)=N} e^{-\beta E(s)} \\
&= \sum_N e^{\beta\mu N} Z_{\rm can}(N),
\end{align}
where I wrote $Z_{\rm g.c.}(\mu)$ for the grand canonical partition function (with chemical potential $\mu$) and $Z_{\rm can}(N)$ for the canonical partition function (with $N$ particles) and $s$ for the microstates.

Note that the equation ($\star$) is the analogue of the one you wrote. 
So, to answer directly the question you asked: "What is the meaning of this $\mu$?", it is the value of the chemical potential that is required for the equivalence of the canonical and grand-canonical ensembles to hold. 
Note: All this is explained much better and with more detail, but still on an informal level, in the first chapter (Introduction) of our book.
A: By a "fixed" value of $N$ here we mean that for the system at equilibrium, the value of $N$ is not allowed to fluctuate. So yes, $N$ is fixed for a canonical ensemble. 
If you allow particle number to change to new value so that the system finds a new equilibrium state with a new fixed value of $N$, the system has a new and different canonical ensemble. But you can still ask: How did the free energy of the system change when we changed the particle number?  The answer to this question is equal to the chemical potential of the system.  
A: Some excellent answers so far. I just wanted to reiterate that the Helmholtz potential is a state function of the volume, number of particles and temperature, and this definition holds without any connection to an ensemble. Hence, the chemical potential formula you quote is always true and in the special case of your initial and final states being equilibrium NVT distributions that only differ in the number of particles, you can calculate the free energy difference between them solely by integrating the chemical potential over the particle difference.
