In my study of quantum mechanics thus far, I have not yet encountered the Dirac equation, but to the best of my knowledge, the Dirac equation is the first place where you can show mathematically that an electron has a spin of $\pm\hbar/2$. If a relativistic version of QM (i.e. Dirac's equation) is the first place that you can determine a particle's spin, then why does the notion of spin arise when you are considering the eigenvalues of the non-relativistic operators $L_z$ and $L^2$ (of course, when you look at the eigenvalues of these operators, you rename them as $J_z$ and $J^2$ to take into account the spin of the system)? I have learned about spin, and know that it arises from the need to rotate the components of some spinor valued wavefunction, but it seems strange that the notion of half-integer eigenvalues of spin would arise in a study of quantum mechanics before the wave-function is even considered to be a spinor, or even before a relativistic formulation of quantum mechanics is considered. So then, why do the eigenvalues of spin show up before considering spinor wavefunctions and before the Dirac equation is considered?
[Physics] Why does spin arise in non-relativistic quantum mechanics
quantum mechanicsquantum-spinspinors
Related Solutions
How does spin arise in (relativistic or not) quantum mechanics?
What are particles in the first place? And what precisely is this property that we call "spin"?
A modern way of arriving at the notion of particles, that may be more transparent than the historical version you summarize in your post, is the way they are introduced in Weinberg's textbook on Quantum Field Theory.
Take the Hilbert space $\mathcal{H}$ of your quantum theory, and take the group $G$ of symmetries of your spacetime. To ensure that an experiment will give the same results if I move it around in spacetime, rotate it, or take it with me on a cruise, there should exist a unitary representation $U$ of $G$ on $\mathcal{H}$. Ie. for any $g \in G$, there should exist a unitary operator $U(g)$ on $\mathcal{H}$, with: $$ U(1) = \text{id}_{\mathcal{H}} \;\&\; U(g.h) = U(g) U(h). $$
We can then decompose this representation $\mathcal{H}, U$ into simpler representations. This means writing $\mathcal{H}$ as a direct sum of smaller Hilbert spaces: $$ \mathcal{H} = \bigoplus_k \mathcal{H}_k $$ with each $\mathcal{H}_k$ being stabilized by all $U(g)$ for $g \in G$, so that $\mathcal{H}_k, U_k := \left. U \right|_{\mathcal{H}_k}$ is itself a unitary representation of G.
A representation that does not contain any smaller representation is called an irreducible representation, or irrep, and the simplest irreps are the ones that hold quantum states with just one particle (refer to Weinberg for what "simplest" exactly means here). So each (simple) irrep of the symmetry group $G$ that can be found in our theory $\mathcal{H}$ is what we define as a particle species.
Integer spins
Good, so if we want to know which particle species are physically possible, we just need to know what are the "simplest" representations of our group of symmetries. Fortunately, a full classification of those is known for the Poincaré or Galilean group. The way it is constructed would be too long to be reproduced here, but again it can be found in great details in Weinberg for the Poincaré group (brief accounts of both the Poincaré and the Galilean case can be found on wikipedia). The bottom line is that, for physically-admissible massive particles, they have the form: $$ \mathcal{H}_k = \text{Span} \left\{ \left| \vec{p}, m \right\rangle \middle| \vec{p} \in \mathbb{R}^3, m \in \mathbb{Z}, -s_k \leq m \leq +s_k \right\} $$ with the non-negative integer $s_k$ being what is called the spin of this particle species $k$ and $\vec{p}$ being its impulsion. The impulsion determines how the particle transform under a spatial translation: $$ U_k(\text{translation by } \vec{\alpha}) \left| \vec{p}, m \right\rangle = e^{i \vec{\alpha}.\vec{p}} \left| \vec{p}, m \right\rangle $$ (just as in the good old momentum representation of QM). Its spin number determines how it transforms under a rotation: $$ U_k(\text{rotation by } \vec{\theta}) \left| \vec{p}, m \right\rangle = \sum_{m^\prime} R^{(s_k)}_{mm^\prime}(\vec{\theta}) \left| R(\vec{\theta}) \vec{p}, m^\prime \right\rangle $$ with $R(\vec{\theta})$ the usual $3\times 3$ rotation matrix acting on $\vec{p}$ and $R^{(s_k)}$ an irrep of the rotation group $\mathcal{SO}(3)$.
The intuition behind this form of $U_k$ is that the spin $s_k$ captures the way the particle may be affected by a rotation beyond the obvious rotation of its impulsion $\vec{p}$. The classical analogy here is that of a rigid body, which changes not only its position but also its orientation under a rotation.
The reason the spin is an integer is because the irreps of $\mathcal{SO}(3)$ are labeled by integers: spin-0 is the trivial representation $R^{(0)}(\vec{\theta}) = \text{id}, \forall \vec{\theta} \in \mathbb{R}^3$ on a 1-dimensional vector space, spin-1 is the usual representation by $3 \times 3$ matrices, and so on.
Half-integer spins
But now there is a twist (figuratively and mathematically...). As mentioned in an earlier comment by ACuriousMind, and as explained in great details in the linked thread, the overall phase of a quantum state is not physically measurable. This means that we can get away with less that a strict unitary representation of $G$ on $\mathcal{H}$, and still ensure that all experimental results are invariant under $G$! Namely, we can replace: $$ U(g\cdot h) = U(g) U(h) \text{ by } U(g\cdot h) = e^{i \varphi(g,h)} U(g) U(h) $$ with the phase factors $\varphi(g,h)$ satisfying suitable consistency relations. Such unitary representations "up to additional phase factors" are called projective representations.
If one goes through the math, we find that for the Poincaré/Galilean group this gives a few additional possible irreps, corresponding to particle species with half-integer spin. They correspond to projective representations of the rotation group in which a rotation of $2\pi$ has a non-trivial (albeit non-detectable) action on the quantum state: $$ U_k(\text{rotation by } 2\pi) \left| \vec{p}, m \right\rangle = - \left| \vec{p}, m \right\rangle $$
Observable signature?
But wait! If this extra minus sign is not physically detectable anyway, how do we know that some particles have half-integer spin?
This has to do with the properties of quantum measurments, which will reveal the spectrum (aka. eigenvalues) of the measured observable. We cannot directly observe that the quantum state vector transforms under a projective representation but we can determine it indirectly because it gets imprinted in the spectrum of the angular-momentum operator.
What about classical (non-quantum) mechanic?
Take a classical mechanical system, say, for concreteness, a system of rigid bodies possibly interacting via conservative forces. The phase space of such a system carries a non-projective representation $T$ of the Galilean group (we can check this by writing it explicitly). But this representation is not a linear representation (at best it may be an affine representation, since translations act, well, by translations). So spin in the above sense does not immediately make sense.
Instead, we can do classical statistical physics for this system: ie. write a field equation for a probability distribution $\rho$ on the phase space (which can be seen as the classical counterpart of a quantum mechanical wave-function). The space $\mathcal{P}$ of such probability distributions naturally carries a linear representation $U$ of $G$ defined by: $$ \forall g \in G,\; [U(g) \rho](x) = \rho\big(T(g^{-1})x\big) $$ which is, again, a non-projective representation (strictly speaking, admissible probability distributions are positive and normalized, but we can study their spin properties by working in the vector space they span: this is analogous to considering the whole Hilbert space in quantum mechanics, although actual quantum states should be normalized).
So, what would a "half-integer-spin mode" for such a system be? According to the previously explained definition of spin, that would be a half-integer-spin irrep $\mathcal{P}_k \subseteq \mathcal{P}, U_k := \left. U \right|_{\mathcal{P}_k}$ appearing in the decomposition of $\mathcal{P}, U$. Can such a $\mathcal{P}_k$ exist? No!
Indeed, if it would, we would have a distribution $\rho \in \mathcal{P}_k \setminus \{0\}$ such that $$ U(\text{rotation by } 2\pi) \rho = U_k (\text{rotation by } 2\pi) \rho = -\rho, $$ but, since $U$ is a non-projective representation, we already know that $U(\text{rotation by } 2\pi) \rho = \rho$.
A similar argument can be applied for example to the classical electromagnetic field: the space of solutions of Maxwell's equations carries a non-projective linear representation of the Poincaré group (one could say: by historical definition of the latter).
What about a thermodynamical system?
Suppose I take a large number of mechanical bodies interacting via conservative forces (say molecules) and take the thermodynamical limit to derive effective equations for some macroscopic variable (eg. their density). Could such an equation exhibit half-integer-modes? Ie. could its space of solutions carry a projective representation of G? Let us do some though experiment:
Take two rigorously identical boxes containing this thermodynamical system and perform the exact same experiment on them, except that the second one is first subjected to a full $2\pi$-rotation (very slowly, so as to not perturb any (local) thermodynamical equilibrium). Because the underlying microscopic theory carries a non-projective representation of G, the two experiments should give the exact same result!
Note that in arguments of this kind, one has to be very careful. The thermodynamic limit can do funny things to the symmetries of a systems. This is known as symmetry-breaking: while the space of solutions of the underlying microscopic theory may be invariant under a certain group $G$, a given thermodynamical phase may have less symmetry because it fails to explore the full solution space (keyword: ergodicity, or more precisely lack thereof).
But, such a mechanism cannot turn a non-projective representation into a projective one: since a $2\pi$-rotation brings me back on the exact same microscopic configuration from which I started, I am guaranteed not to land in a different thermodynamical phase.
Can we cheat?
Suppose I come up with a mathematical description of some physically valid classical system in which, for technical reasons, I choose to introduce some auxiliary, non-measurable quantity (eg. a complex phase). Since the auxiliary quantity is non-measurable, I can let it transform in whatever way is mathematically convenient. In this way, I may arrive at a description of a classical system which exhibits a projective representation.
But still the original physical system will not exhibit any observable half-integer-spin behavior. As the truly measurable quantities have to be invariant under a full $2\pi$-rotation, there should exist a basic description of the same system, that refrains from introducing any auxiliary quantities, and carries a non-projective representation. Computing experimental predictions using this basic description, no half-integer-spins should show up.
TL;DR: This is crucially differently from the above discussed quantum mechanical case, in which you can hide a non-projective representation, so as to preserve $2\pi$-rotation-invariance, while nevertheless retaining some observable signature.
Bonus: Does Wick-rotating a quantum equation give a thermodynamical equation?
I do not think Wick rotation should be thought as some kind of magic transformation to turn a QM equation into a thermodynamical one.
There is a connection between quantum field theory on 3+1d Minkowski spacetime and statistical field theory in 4d Euclidean space. But statistical physics (the study of the probability distribution over (field) configurations) is not quite the same as thermodynamics (the derivation of effective equations for macroscopic variables in the large-number-of-particles limit).
I suspect the appearance of the heat equation as the complex-time Schrödinger equation is more a coincidence, coming from the fact that, well, there are only so many linear PDEs you can write with a certain order in space and time derivatives.
If you would like to investigate the Wick rotated Dirac equation anyway, I guess a good place to start would be the Wick-rotated gamma matrices. You will get a field equation carrying a projective representation of the 4d Euclidean group, sure. But Wick-rotating a physically valid quantum equation does not a priori guarantee any particular physical relevance for the resulting equation: in fact, such an equation cannot describe any actual physical system, if only because, as pointed out by flippiefanus, we do not live in 4d Euclidean space ;-).
There is no such thing as spin $1/4$ in 4-dimensional spacetime.
Spin has to do with representations of the Lorentz Lie algebra. So let me shed some light on how these are classified.
First, the complex-valued Lorentz algebra $so(1,3)$ is equivalent to $so(4)$. The $so(4)$ algebra is equal to the direct sum of two copies of $su(2)$, which means that irreducible representations of $so(1,3)$ are labeled by ordered pairs of the irreducibles of $su(2)$.
The $su(2)$ representation theory can be found in any textbook on Lie groups or even in some QFT textbooks. One of the crucial facts is that irreducibles of $su(2)$ are labeled by nonnegative half-integers called spins:
$$ j = 0, \; 1/2, \; 1, \; 3/2, \; 2, \; 5/2, \; \dots $$
This is enough to start constructing irreducibles of the Lorentz algebra. The basic blocks of the representation theory are the two fundamental representations $(1/2, 0)$ and $(0, 1/2)$ which are called the left and right Weyl spinors respectively. Both are two-dimensional.
Dirac spinors actually belong to the 4-dimensional reducible representation $$ (1/2, 0) \oplus (0, 1/2). $$
Another 4-dimensional representation is the irreducible $(1/2, 1/2)$ to which the 4-vectors belong.
As you can see, this is all just representation theory and there is no spin-$1/4$ representation in 4 spacetime dimensions.
However, in $2$ spacetime dimensions this is no longer valid and representations with fractional spins exist.
Best Answer
In quantum mechanics, the $\overrightarrow L$ operator and the associated observables $L^2$ and $L_z$ already existed in nonrelativistic QM. The eigenfunctions of $L^2$ and $L_z$ are the spherical harmonics and the eigenvalues are $\hbar^2l(l+1)$ and $\hbar m$ respectively. The values of $l$ can be integer (corresponding to the spherical harmonics) and half integer (corresponding to a divergent solution to the differential equation). Given a value of $l$, the allowed values of $m$ went from $-l$ to $l$ in integer steps, totalling $2l+1$ degeneracy.
You can use the fact that experimental observables are related to the angular momentum operators. In particular, the energy of a dipole in a magnetic field is proportional to $L_z$. When experiments were performed on hydrogen, they observed a number of effects that hinted at electrons having spin. Firstly, electrons occupied hydrogen orbitals (and other atom's orbitals) in pairs. The Pauli exclusion principle forbids electrons from occupying the same state, but if you add a new quantum number with two states, you can explain the behavior of the periodic table. Experiments with Hydrogen in magnetic fields showed that the energy levels diverged. There was the expected behavior of moving charge creating a magnetic dipole. This produced odd splittings in the energy levels (as was expected for moving electrons in spherical harmonic states). But there were also splitting of the energy states which produced an even effect. A last effect worth mentioning and the one that really demonstrates that electrons have intrinsic spin $\frac{1}{2}\hbar$ is that if you send a cathode ray beam through a magnetic field, it will split into two beams in the direction of the magnetic field, indicating that the electron was a magnetic dipole with only two spin states.
These and a few other observations lead physicists to the conclusion that electrons have intrinsic angular momentum of magnitude $\frac{1}{2}\hbar$.