What do the Pauli matrices mean?

Let me first remind you of (or perhaps introduce you to) a couple of aspects of quantum mechanics in general as a model for physical systems. It seems to me that many of your questions can be answered with a better understanding of these general aspects followed by an appeal to how spin systems emerge as a special case.

General remarks about quantum states and measurement.

The state of a quantum system is modeled as a unit-length element $|\psi\rangle$ of a complex Hilbert space $\mathcal H$, a special kind of vector space with an inner product. Every observable quantity (like momentum or spin) associated with such a system whose value one might want to measure is represented by a self-adjoint operator $O$ on that space. If one builds a device to measure such an observable, and if one uses that device to make a measurement of that observable on the system, then the machine will output an eigenvalue $\lambda$ of that observable. Moreover, if the system is in a state $|\psi\rangle$, then the probability that the result of measuring that quantity will be the eigenvalue of the observable is \begin{align} p(\lambda) = |\langle \lambda|\psi\rangle|^2 \end{align} where $|\lambda\rangle$ is the normalized eigenvector corresponding to the eigenvalue $\lambda$.

Specialization to spin systems.

Suppose, now, that the system we are considering consists of the spin of a particle. The Hilbert space that models the spin state of a system with spin $s$ is a $2s+1$ dimensional Hilbert space. Elements of this vector space are often called "spinors," but don't let this distract you, they are just like any other vector in a Hilbert space whose job it is to model the quantum state of the system.

The primary observables whose measurement one usually discusses for spin systems are the cartesian components of the spin of the system. In other words, there are three self-adjoint operators conventionally called $S_x, S_y, S_z$ whose eigenvalues are the possible values one might get if one measures one of these components of the system's spin. The spectrum (set of eigenvalues) of each of these operators is the same. For a system of spin $s$, each of their spectra consists of the following values: \begin{align} \sigma(S_i) = \{m_i\hbar\,|\, m_i=-s,-s+1,\dots, s-1,s\} \end{align} where in my notation $i=x,y,z$. So for example, if you build a machine to measure the $z$ component of the spin of a spin-$1$ system, then the machine will yield one of the values in the set $\{-\hbar, 0, \hbar\}$ every time. Corresponding to each of these eigenvalues, each spin component operator has a normalized eigenvector $|S_i, m_i\rangle$. As indicated by the general remarks above, if the state of the system is $|\psi\rangle$, and one wants to know the probability that the measurement of the spin component $S_i$ will yield a certain value $m_i\hbar$, then one simply computes \begin{align} |\langle S_i, m_i |\psi\rangle|^2. \end{align} For example, if the system has spin-$1$, and if one wants to know the probability that a measurement of $S_y$ will yield the eigenvalue $-\hbar$, then one computes \begin{align} |\langle S_y, -1|\psi\rangle|^2 \end{align}

Spinors.

In the above context, spinors are simply the matrix representations of states of a particular spin system in a certain ordered basis, and the Pauli spin matrices are, up to a normalization, the matrix representations of the spin component operators in that basis specifically for a system with spin-$1/2$. Matrix representations often facilitate computation and conceptual understanding which is why we use them.

More explicitly, suppose that one considers a spin-$1/2$ system, and one chooses to represent states and observables in the basis $B =(|S_z, -1/2\rangle, |S_z, 1/2\rangle)$ consisting of the normalized eigenvectors of the $z$ component of spin, then one would find the following matrix representations in that basis \begin{align} [S_x]_B &= \frac{\hbar}{2}\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} = \frac{\hbar}{2}\sigma_x\\ [S_y]_B &= \frac{\hbar}{2}\begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} = \frac{\hbar}{2}\sigma_y\\ [S_z]_B &= \frac{\hbar}{2}\begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} =\frac{\hbar}{2}\sigma_z\\ \end{align} Notice that these representations are precisely the Pauli matrices up to the extra $\hbar/2$ factor. Moreover, each state of the system would be represented by a $2\times 1$ matrix, or "spinor" \begin{align} [|\psi\rangle]_B = \begin{pmatrix} a \\ b\end{pmatrix}. \end{align} And one could use these representations to carry out the computations referred to above.

Groups are abstract mathematical structures, defined by their topology (in case of continual (Lie) groups) and the multiplication operation.

But it is almost impossible to talk about abstract groups. That is why usually elements of groups are mapped onto linear operators acting on some vector space $V$:

$$ g \in G \rightarrow \rho(g) \in \text{End}(V), $$

where G is the group, $\text{End}(V)$ stands for endomorphisms (linear operators) on $V$, and $\rho(g)$ is the mapping. In order for this mapping to be meaningful, we have to map the group multiplication properly:

$$ \rho(g_1 \circ g_2) = \rho(g1) \cdot \rho(g2). $$

The inverse is also mapped to

$$ \rho(g^{-1}) = \rho(g)^{-1} $$

and the group identity is just

$$ \rho(e) = \text{Id}_V. $$

This is called the representation of the group $G$. $V$ transforms under the representation $\rho$ of group $G$.

In your case, the group of interest is the group of rotations in 3 dimensions which is usually denoted as SO(3). Our goal is to find different objects which can be rotated, i.e. representations (and representation spaces) of SO(3).

One such representation is the defining representation (which is used to define SO(3)), or the vector representation. In this case $V$ is just $R^3$ and matrices from $\rho(\text{SO(3)})$ are orthogonal $3\times 3$ matrices with unit determinant:

$$ A^{T} A = 1;\quad \det A = 1 $$

So vectors can be rotated in 3 dimensions. The result of such rotation by $g \in \text{SO(3)}$ is determined by acting on the initial vector with the operator $\rho(g)$.

Another representation is the spinor representation. The vector space is now 2-dimensional and complex. The image of this representation consists of unitary $2\times 2$ with unit determinant:

$$ A^{\dagger} A = 1;\quad \det A = 1. $$

This representation is not as obvious as the previous one, since spinors are something which we don't usually see in everyday life. But it can be mathematically proven that these representations are isomorphic and therefore are two different representations of the same group (actually, they are homomorphic and spinor representation is the double cover of the vector representation).

Now to the Pauli matrices. There is a general principle: for each Lie group $G$ there exist a corresponding linear space (Lie algebra) with a Lie bracket (an anti-commutative operation satisfying the Jacobi identity) which maps uniquely onto some neighbourhood of the group unity of $G$. This mapping is called the exponential.

So you can write an arbitrary (close enough to unity for global topological problems to be avoided) $2\times 2$ complex matrix from the spinor representation in form

$$ A = \exp \left[ \frac{i}{2} \alpha^a \sigma_a \right], $$

where $\alpha^a$ are three numbers which parametrize the group element whose representation is $A$, and $\frac{i}{2} \sigma_a$ are the Lie algebra basis, with $\sigma_a$ - 3 $2\times 2$ Pauli matrices. This equation pretty much specifies how a spinor is transformed under an arbitrary rotation.

In the vector representation there is also a Lie algebra basis, which consists of 3 $3\times 3$ matrices.

There are two other interpretation of the Pauli matrices that you might find helpful, although only after you understand JoshPhysics's excellent physical description. The following can be taken more as "funky trivia" (at least I find them interesting) about the Pauli matrices rather than a physical interpretation.

1. As a Basis for $\mathfrak{su}(2)$

The first interpretation is variously seen as (i) they are unit quaternions, modulo a sign change and reordering of the mathematician's definition of these beasts, (ii) as a basis for the Lie algebra $\mathfrak{su}(2)$ of $SU(2)$ when we use the matrix exponential to recover the group $SU(2) = \exp(\mathfrak{su}(2))$ through (iii) a three dimensional generalisation of De Moivre's theorem.

A general, traceless, $2\times2$ skew Hermitian matrix $H$ can be uniquely decomposed as:

$$H = \alpha_x \sigma_x + \alpha_y \sigma_y + \alpha_z \sigma_z\tag{1}$$

with $\alpha_x,\,\alpha_y,\,\alpha_z\in \mathbb{R}$. This matrix fulfils the characteristic equation $H^2 = -\frac{\theta^2}{4}\,\mathrm{id}$, where $\mathrm{id}$ is the $2\times2$ identity and $\frac{\theta}{2} = \sqrt{\alpha_x^2+\alpha_y^2+\alpha_z^2}$. So, if we deploy the universally convergent matrix exponential Taylor series, and then reduce all powers of $H$ higher than the linear term with the characteristic equation, we get:

$$\exp\left(H\right) = \cos\left(\frac{\theta}{2}\right)\mathrm{id} + \hat{H}\sin\left(\frac{\theta}{2}\right)\tag{2}$$

which is seen to be a generalisation of De Moivre's formula for the "pure imaginary" unit

$$\hat{H} = \frac{\alpha_x \sigma_x + \alpha_y \sigma_y + \alpha_z \sigma_z}{\sqrt{\alpha_x^2+\alpha_y^2+\alpha_z^2}}\tag{3}$$

and all members of $SU(2)$ can be realised by an exponential such as in (2) (but be aware that the exponential of a Lie algebra, although the whole of $SU(2)$ in this case, is not always the whole Lie group unless the latter is (i) connected and (ii) compact). Thus every member of $SU(2)$ can be decomposed as a "unit length superpositon of the Pauli matrices and the identity matrix.

The reason for the factor 2 in the definition $\theta/2$ is so far mysterious: witness that for the purposes of the above, we might just as easily have replaced $\theta/2$ by $\theta$. The reason is related to the relationship between the Pauli matrices and the Celestial sphere, which I discuss later on. Quaternions represent rotations through a spinor map (BUT, as Joshphysics advises, don't be distracted too much by this word); if a vector in 3-space is represented by a purely imaginary quaternion of the form $x\,\sigma_x+y\,\sigma_y+z\,\sigma_z$, then its image under a rotation of angle $\theta$ about an axis with direction cosines $\gamma_x,\,\gamma_y,\,\gamma_z$ is given by:

$$x\,\sigma_x+y\,\sigma_y+z\,\sigma_z \mapsto U\,(x\,\sigma_x+y\,\sigma_y+z\,\sigma_z)\,U^\dagger;\quad U=\exp\left(\frac{\theta}{2}(\gamma_x\,\sigma_x+\gamma_y\,\sigma_y+\gamma_z\,\sigma_z)\right) \tag{4}$$

This spinor map is an example of the group $SU(2)$ acting on its own Lie algebra through the adjoint representation. It can be intuitively understood in terms of a triangle rule to work out the compositions of two rotations, as sketched in my diagram below. The arcs on the unit sphere represent a rotation through an angle twice that given by the angle subtended by the arc at the origin.

I explain this in detail in Example 1.4 "$2\times2$ Unitary Group $SU(2)$" on my web page "Some Examples of Connected Lie Groups" here.

There is also my interactive Mathematica demonstration "The $SU(2)$ Spinor Map: Rotation Composition by Graphical Quaternion Triangles" on the Wolfram Demonstrations site.

2. The Celestial Sphere

By expanding the 3 dimensional linear space of superpositions of Pauli matrices (which is the same as the linear space of traceless $2\times2$ skew-Hermitian matrices) to the 4 dimensional space spanned by the Pauli matrices and the identity matrices, then any transformation from the group $SL(2,\,\mathbb{C})$ acts on vectors of the form $t\,\mathrm{id}+x\,\sigma_x + y\, \sigma_y + z\,\sigma_z$ by the same spinor map as in (4). If we restrict ourselves to projective rays in this space, the group $SL(2,\,\mathbb{C})$, isomorphic to the Moebius group of Möbius transformations acts on this space of rays in exactly the same way as Möbius (fractional linear) transformations act on the Riemann sphere. $SL(2,\,\mathbb{C})$ is a double cover of the Lorentz group, and you can calculate how the view of a spacefarer changes as they undergo Lorentz transformations. See the section "Lorentz Transformations" on the Wikipedia "Möbius Transformation" page for further details.