1) The answer depends on what you mean by rigorous; obviously the 1-d derivation on wikipedia is not exactly rigorous. It also depends on what you would like to use as a micro-physical starting point, and how much detail you desire. Hundreds of text books have been written on deriving diffusive laws from (quantum) kinetic theory, linear response theory, etc.
2) The most basic argument is still one of the best. If the concentration is inhomogeneous we expect to see a current, and $\vec\nabla u$ is the only vector available. Then
$$
\vec\jmath = - D\vec\nabla u + O(\nabla^2)\, ,
$$
which is Fick's law. The modern view of this is called effective (field) theory: Simply write down the most general expression for the current as a power series in gradients of the thermodynamic variables.
3) The standard kinetic theory derivation goes roughly like this: Consider a chemical potential $\mu$ conjugate to $u$. Then the distribution function is
$f_0\sim \exp(-\mu(x)/T)$, where we have allowed for gradients in the concentration. If the inhomogeneity is weak, we can seek solutions of the Boltzmann equation of the form $f_0+\delta f$. Insert into the Boltzmnann equation and expand. We get
$$
\frac{1}{T}\vec{v}\cdot\vec\nabla\mu \,f_0 = -\frac{\delta f}{\tau}
$$
where I write the linearized collision operator in terms of its smallest
eigenvalue, the inverse collision time $1/\tau$. ($v$ is the quasi-particle velocity). I can directly solve for $\delta f =-\frac{\tau}{T}\vec{v}\cdot\vec\nabla\mu\, f_0$. Now compute the current generated by $\delta f$:
$$
\vec\jmath = \int d^3p \,\vec{v}\delta f = - \frac{\tau}{3mT}\left(\int d^3p
\,v^2 f_0\right) \vec\nabla\mu
$$
Use thermodynamic identities to relate $\nabla\mu$ to $\nabla u$, and compute the integral in terms of the total pressure or energy of the gas. Presto, we get Fick's law, plus an explicit result for $D$.
