"No theoretical proof of the Pauli's Exclusion Principle can be given as yet and for the present it must be regarded as something empirical added to and regulating the vector atom model."
I've found it in Atomic & Nuclear Physics by N. Subrahmanyam & Brij Lal.
My question is "What was the motivation behind Pauli's Exclusion Principle? Is that just an ad-hoc intuitive attempt of Pauli to explain the Zeeman effect or it was derived on the basis of mathematics? If mathematical then how? "
Initially Pauli exclusion principle imposed for fermions was purely phenomenological conception introduced in order to explain experimental facts. The reason is that there is no well-defined theoretical description of the relation between the spin and the statistics in non-relativistic quantum mechanics (within which Pauli has formulated his principle initially). But it arises in relativistic quantum theory. I'll try to explain these things below.
Qualitatively about the statistics
In the 3-dimensional space there are only two possibilities for the behavior of the wave-function $\psi(\mathbf x_{1},\mathbf x_{2})$ of the two identical particles under adiabatical changing of their positions due to which the particles becomes to be interchanged. Under this action, the two-particle wave-function can only be changed as $$ \tag 1 \psi(\mathbf x_{1},\mathbf x_{2}) \to \pm \psi(\mathbf x_{2},\mathbf x_{1}) $$ Note that this property follows from the topology of the relative phase space of two-particle state in the 3-dimensional space, and has nothing to do with any other arguments (in particular with experiment). This is purely theoretical argument.
Qualitatively about the spin
The spin is the quantity whose importance fundamentally emerges because of the Poincare symmetry of our world. Apart from its physical sense, it is the quantum number by which we mathematically characterize each (at least massive) particle due to its transformation properties under the Poincare group. Through the representations of the Poincare group, the description of the spin is related to the topology of the space-time.
The relation between the spin and statistics
So where these two conceptions, the spin and the statistics, meet each other? How the spin affects the statistics (and vice versa) and, in particular, how the Pauli exclusion principle follows from this? From the first point of view any relation between them is unnatural, at least from the point of view of topology. However, actually the relation exists.
The poincare invariance of the quantum theory requires that the hamiltonian density $\hat{H}(x)$ of the theory must commute with itself for space-like intervals: $$ \tag 2 [\hat{H}(x),\hat{H}(y)] = 0 \quad \text{for} \ (x-y)^2<0 $$ The hamiltonian is composed from the field operators $\hat{\psi}_{a}(x), \hat{\psi}_{b}(y)$ quantized in terms of creation-destruction operators $\hat{a}(\mathbf p, s),\hat{a}^{\dagger}(\mathbf p, s)$ of the particles with arbitrary spin $s$. From the relation $(1)$ we know that the creation-destruction operators must obey $$ \tag 3 [\hat{a}(\mathbf p, s),\hat{a}(\mathbf k, s)]_{\pm} =0,\quad [\hat{a}(\mathbf p, s),\hat{a}^{\dagger}(\mathbf k, s)]_{\pm} \sim \delta(\mathbf p -\mathbf k) $$ Note that this statement is purely theoretical and has nothing to do with the phenomenology. By taking into account both $(2),(3)$, we obtain that $$ \tag 4 [\hat{\psi}_{a}(x),\hat{\psi}_{b}(y)]_{\pm} = 0, \quad (x-y)^{2} < 0 $$ The structure of the field operators are completely fixed by their transformation properties and in particular by the spin value. The expression $(4)$ is the place where the spin and statistics meet each other. By treating it analytically, one obtains the Pauli exclusion principle.
Why there is no the relation between the spin and the statistics in non-relativistic quantum mechanics?
In the spirit of the statements written above it's not hard to understand why in non-relativistic physics the spin-statistics relation doesn't have theoretical base and can be only phenomenological. The reason is that in non-relativistic physics there is no requirement similar to $(2)$. Really, ortochronous Galilei group transformation, which ("plus" translations) represents the space-time symmetry in non-relativistic quantum mechanics, leaves unchanged chronological ordering in the S-operator. Contrary to the Galilei group, the Poincare group changes chronological ordering for space-like intervals. The latter is the underlying reason due to which we require $(2)$...
The Pauli exclusion principle can be derived from relativistic quantum field theory. Even though the Pauli exclusion principle is still important in the non-relativistic approximation, the reason for it is deeply rooted in relativistic QFT.
The theorem is called the spin-statistics theorem. The inputs to the theorem include
Lorentz symmetry,
Something called the spectrum condition, which says that the total energy must have a finite lower bound.
The result is easiest to derive for the case of a free (non-interacting) spin-1/2 particles. This is explained in more detail in Name YYY's answer. Here's a words-only version that includes a little extra context. In QFT, particles are phenomena that the theory predicts, rather than ingredients used to construct the theory. The theory is constructed in terms of quantum fields, which are represented as operators acting on a Hilbert space. For example, all electrons (and antielectrons) are manifestations of the electron field, just like all photons are manifestations of the electromagnetic field. For a free spin-1/2 field, Lorentz symmetry requires that the equation of motion be linear in the space- and time-derivatives, rather than being quadratic like it would be for a scalar or vector field. Because of this, the field operators need to anticommute with each other at spacelike separation (rather than commuting with each other as they would for a scalar or vector field) in order for the energy operator to have an energy spectrum with a finite lower bound. This anticommutativity is the Pauli exclusion principle. This derivation, in the special case of a free spin-1/2 field, is included in many QFT textbooks, such as section 3.5 in Peskin and Schroeder's An Introduction to Quantum Field Theory.
For general derivations:
The classic reference is the book by Streater and Wightman, PCT, Spin and Statistics, and All That (1980). They prove the theorem starting from a system of axioms for relativistic QFT called the Wightman axioms. They also prove another general theorem that is usually called the CPT theorem (they called it "PCT" instead of "CPT" — same thing), which is often presented together with the spin-statistics theorem because they share the same inputs.
It can also be derived from a different system of axioms in the context of algebraic QFT. The relevant theorem in this case is called the Doplicher-Roberts reconstruction theorem, which actually proves much more: it helps explain why field operatores are useful in the first place, in addition to explaining the Pauli exclusion principle. This theorem is reviewed on page 92 in "Algebraic Quantum Field Theory", https://arxiv.org/abs/math-ph/0602036.
From my research, the short answer to your question is Pauli came up with the principle to explain more on what was happening within an atom, but he himself could not explain where it comes from. The Pauli exclusion principle has no derivation. While Wolfgang Pauli was trying to answer the big quantum question at the time, why don't all the electrons go into the lowest energy state, he had come up with the Pauli exclusion principle. At the time many physicists were confused by this because they had no idea where it even came from, and why it even exists in the first place. In fact, Pauli himself was troubled that he couldn't explain his own principle from logic and that he couldn't derive it from any quantum equations. This is a common theme within quantum physics, where there are many laws and equations that have no proper explanation on why they happen, and scientists are still trying to find more answers.
My source: https://www.aps.org/publications/apsnews/200701/history.cfm
Fermions are described by antisymmetric wave functions. $$\Psi(x_{1},x_{2})=-\Psi(x_{2},x_{1})\tag{1}$$ We define: $$\Psi(x_{1},x_{2})=\frac{1}{\sqrt{2}}\left[\Psi_{1}(x_{1})\Psi_{2}(x_{2})-\Psi_{1}(x_{2})\Psi_{2}(x_{1})\right]\tag{2}$$ Using this definition we can see that: $$\Psi(x_{2},x_{1})=\frac{1}{\sqrt{2}}\left[\Psi_{1}(x_{2})\Psi_{2}(x_{1})-\Psi_{1}(x_{1})\Psi_{2}(x_{2})\right]\tag{3}$$ So equation (1) is satisfied. We can write equation (2) as a determinant: $$\Psi(x_{1},x_{2})=\frac{1}{\sqrt{2}}\begin{vmatrix} \Psi_{1}(x_{1}) & \Psi_{2}(x_{1})\\ \Psi_{1}(x_{2}) & \Psi_{2}(x_{2}) \end{vmatrix}$$ We can write this determinant for any number of fermions: $$\Psi(x_{1},x_{2},\dots,x_{N})=\frac{1}{\sqrt{N!}}\begin{vmatrix} \Psi_{1}(x_{1}) & \Psi_{2}(x_{1}) & \cdots & \Psi_{N}(x_{1})\\ \Psi_{1}(x_{2}) & \Psi_{2}(x_{2}) & \cdots & \Psi_{N}(x_{2})\\ \vdots & \vdots & \ddots & \vdots\\ \Psi_{1}(x_{N}) & \Psi_{2}(x_{N}) & \cdots & \Psi_{N}(x_{N})\\ \end{vmatrix}$$ This is known as a Slater determinant. You can see that if two or more fermions share the same quantum state, the determinant is zero.
At the Center of all that is the so-called "Spin-Statistics-Theorem" which says that particles with integer spin, aka. bosons, must follow Bose-Statistics and particles with half-integer spin, aka. fermions, follow Fermi-Statistics, which include the Pauli exclusion principle. That theorem follows from very general statements like Lorenz invariance, locality, unitarity, positive norm, finite energy, ... To make all that a bit more mathematical, consider the so called "construction of the Fock state", which is part of standard QFT courses. One starts with a single vacuum state, $|0>$, and acts on it with creation operators $a'_p$ (it should be the adjoint of a, but I really don't know how to write that here). For bosons, the operators $a'_p$ commute, while for fermions, they are often called $c'_p$ and do anticommute, so $c'_kc'_p=-c'_pc'_k$. That is because otherwise, there would be negative energy states etc. and the theory would not work out. What they do is create a particle with momentum p or k, but they may also contain other information about the state. So, a one-boson state with momentum p is $$a'_p|0>$$ and a two boson state where both bosons are of momentum p is proportional to $$a'_pa'_p|0>$$ Now, what happens if we want to create two fermions in the same state? That would be $$c'_kc'_k|0>=-c'_kc'_k|0>=0$$ because the creation operators anticommute. That is a simple (although admittably not very general) prove of the Pauli exclusion principle.
