Dependence of spin on classical vs non-classical physics? Textbook derivations often state that spin can be derived by adding relativity to quantum mechanics. The general argument comes in several steps : 


*

*Schrödinger first tried to describe quantum particles by quantizing $E^2 = p^2 + m^2$, which leads to the Klein-Gordon equation (KGE). The density associated to solutions of KGE can be negative, which is problematic if density is interpreted as a probability density. Ultimately, Schrödinger dropped the relativistic constraint and derived his equation using semi-classical approximations to the Hamilton-Jacobi equation.

*Pauli heuristically modified the Schrödinger equation to account for spin, leading to the so-called Pauli equation

*Dirac understood that treating time and space on the same footing is the only constraint one should impose on a field equation, even if this leads to using operators as coefficients. By linearizing KGE, he obtained the Dirac equation. Interpreting the equation solutions as quantum wave functions allows to explain spin, but also leads to negative energy states.

*Imposing anti-commutation relations on the solution to the Dirac equation, and hence promoting it into the field, solves the negative energies problem.


The conclusion of these steps is generally that spin is derived by going to a relativistic setting. 
However, we know since the work of Lévy-Leblond that this is not true. Indeed, if we interpret irreducible representations of the Galilei group as spin states, the wave equation obtained by imposing Galilei invariance is the Schrödinger equation. But this is not the only invariant equation one can derive ! Factoring the Schrödinger equation to have a linear equation with first order time and space derivatives leads to a non-relativistic equation describing spin 1/2 particles, which is exactly the Pauli equation, and the non-relativistic limit of the Dirac equation ! The Schrödinger equation can thus be seen as an exact evolution equation for spin 0 non-relativistic particles or as an approximation of the evolution of all non-relativistic quantum particles when neglecting spin.
Therefore, the conclusion is that spin (especially spin 1/2) comes from imposing a linear relation on time and space derivatives. 
My questions will start here : 


*

*Am I wrong when thinking that spin has nothing to do with either relativity nor quantum mechanics ? One can define linear equations for classical fields, leading to spinor solutions, which can be non-relativistic.

*We know that the heat equation is related to the Schrödinger one by a Wick rotation. What is the equivalent of a Wick-rotated Pauli or Dirac equation ? Of course, it is expected to be a classical field equation constraining the evolution of macroscopic variables of a thermodynamic system, but which one ?

*Is there a notion of spin attached to the classical field modes solving the previous equation ? How can we test a spin 1/2 property in such case ?

*Is there some deep relationship between linearizing an equation and the Wick rotation ? Both somehow correspond to taking the "square-root of a geometry", using respectively space derivative order reduction and complex time.


Thanks in advance for your help ! 
EDIT : Thanks for the answers, which helped me focus a bit more my questions into a single one : do you know of any paper about the statistical system obtained by Wick-rotating the Dirac or Pauli equation ? My goal is to check if spin 1/2 classical modes exist for such a system and how they could be detected.
 A: Regarding your first question, 

Am I wrong when thinking that spin has nothing to do with eitherrelativity nor quantum mechanics?

let us take spin to mean 'intrinsic angular momentum'. Then we can ask, is there a consistent classical Galilean model of spin? Is there a consistent classical Lorentz-covariant theory of spin? The answer is yes, to both questions.
For the Galilean case, it is quite simple. You can do Hamiltonian mechanics on any manifold that has what is called a symplectic structure. Fortunately, the two-sphere $S^2$ admits precisely such a structure. Thus the dynamics can be defined by, e.g., the Hamiltonian function $$H = \frac{(\mathbf p - q\mathbf  A)^2}{2m} + q\phi + \mu \mathbf B(\mathbf x) \cdot \mathbf S $$
and the Poisson brackets $$\{x_i, p_j\} = \delta_{ij} \quad \{s_i, s_j \} = \epsilon_{ijk} s_k.$$
(See how the Poisson bracket for the spin is the same as the commutation relation for the Pauli matrices? It's of course not a coincidence - it's because $\mathfrak{su}(2) \sim \mathfrak{so}(3)$.) The equations of motion are $$\dot f = \{ f, H\}$$
as usual. You can work out that you get the usual Lorentz force plus a gradient force, and the spin precesses around $\mathbf B$.
Now what about the Lorentz case? Actually, such a theory precedes Pauli's introduction of his famous matrices. It was published in covariant form Frenkel in 1926 [1] and later that year Thomas derived his eponymous precession [2]. Thomas's equation of motion for the spin was later rediscovered by Bargmann, Michel, and Telegdi [3]. Thomas considered only homogeneous fields. R.H. Good [4] and I. Y. Tamm (unfortunately I don't have the original reference at the moment, but the paper is in Russian anyway) extended this to inhomogeneous fields. The Frenkel model and the Tamm-Good equations are reviewed in Refs. [5, 6].
Now the question is of course can you distinguish between classical and quantum models? Well, the Stern-Gerlach experiment [7] does that. The real story is a bit more complicated according to historians of science [8-10]. (I should add a disclaimer that I have not read these papers in detail as I've not had the time. But I think you will find them interesting.) On the other hand, Bell and GHZ experiments are rather convincing that we do live in a quantum world.
We can also look for Stern-Gerlach effects in the relativistic regime by using the Frenkel model/Tamm-Good equations and comparing to a calculation based on the Dirac equation. (To avoid having to use full QED, one should use the Foldy-Wouthuysen transformation [11-14].) This was done recently by Weng, Bauke, and Keitel [15], showing that strongly relativistic electrons in extremely intense laser fields have different motions according to the Frenkel and Foldy-Wouthuysen models. According to Weng, Bauke, and Keitel, the difference is plausibly within experimental limits at least with next-generation lasers.
References
[1] J. Frenkel, Nature 117, 653 (1926).
[2] L. H. Thomas, Nature 117, 514 (1926).
[3] V. Bargmann, L. Michel, and V. L. Telegdi, Phys. Rev. Lett. 2, 435 (1959).
[4] R. H. Good, Phys. Rev. 125, 2112 (1962).
[5] V. G. Bagrov and V. A. Bordovitsyn, Sov. Phys. J. 23, 128 (1980).
[6] I. M. Ternov and V. A. Bordovitsyn, Sov. Phys. Uspekhi 23, 679 (1980).
[7] B. Friedrich and D. Herschbach, Phys. Today 56, 53 (2003).
[8] F. Weinert, Stud. Hist. Philos. Mod. Phys. 26, 75 (1995).
[9] D. Giulini, Stud. Hist. Philos. Mod. Phys. 39, 557 (2008).
[10] M. Morrison, Stud. Hist. Philos. Mod. Phys. 38, 529 (2007).
[11] L. L. Foldy and S. A. Wouthuysen, Phys. Rev. 78, 29 (1950).
[12] J. P. Costella and B. H. J. McKellar, Am. J. Phys. 63, 1119 (1995).
[13] A. J. Silenko, Phys. Rev. A 77, 12116 (2008).
[14] D.-W. Chiou and T.-W. Chen, 1 (2015).
[15] M. Wen, H. Bauke, and C. H. Keitel, Sci. Rep. 6, 31624 (2016).
A: 
Am I wrong when thinking that spin has nothing to do with either
  relativity nor quantum mechanics?

Spin is part of angular momentum (the other part being orbital angular momentum). It comes from the part of the Lorentz group that deals with three-dimensional rotations. As such, you are correct, it does not have anything to do with special relativity, because it does not involve the boosts.
You are also correct that it does not need to involve quantum mechanics, because rotations are also relevant in classical theories. In fact the idea of the commutation relations associated with rotations were already proposed by Hamilton in 1843, long before the advent of quantum mechanics.

We know that the heat equation is related to the Schrödinger one by a
  Wick rotation. What is the equivalent of a Wick-rotated Pauli or Dirac
  equation?

The Wick rotation converts the time-dimension into an Euclidean dimension, as opposed to the Minkowski time-dimension found in special relativity. I don't know whether anybody has ever looked at the Dirac equation in this way, but one can try to think what this would mean or represent. [I'll try to add the math later.] One way to think of the Dirac equation, for instance, is to view it as the `square root' of the Klein-Gordon equation. From this perspective, the Dirac matrices can be derived by requiring that they reproduce the metric tensor when squaring the Dirac equation the recover the Klein-Gordon equation. In the Wick rotated version of the Klein-Gordon equation, one ends up with a 4-dimensional Poisson equation. This means that the Wick rotated version of the Dirac equation would be the square root of this 4-dimensional Poisson equation. The corresponding matrices (analogous to the Dirac matrices) must then reproduce the identity matrix upon squaring. The resulting equation would then represent a kind of spin in 4-dimensions.

Of course, it is expected to be a classical field equation
  constraining the evolution of macroscopic variables of a thermodynamic
  system, but which one?

Not sure I understand where this `of course' comes from, but the connection to thermodynamics is not clear to me. I don't think the analogy between the heat equation and the Schrödinger equation works in all cases with Wick rotations.

Is there a notion of spin attached to the classical field modes
  solving the previous equation? How can we test a spin 1/2 property in
  such case?

Naturally when one imposes a Wick rotation one changes the symmetry properties of the system. One moves from Minkowski space to a 4-dimensional Euclidean space. Hence, the symmetry group changes from the Lorentz group SO(1,3) to SO(4). The former contains 6 generators associated with the three boosts and the three rotations. The latter has also been well studied (see here). It does contain the well-known spin representations in terms of sub groups.

Is there some deep relationship between linearizing an equation and
  the Wick rotation?

I assume that the term linearizing refers to the process that I've called the square root, which for instance relates the Dirac equation to the Klein-Gordon equation. In this case I do not think that this is related to the Wick rotation. The latter changes the properties of the space in which the equations are defined. Neither do I think that the square root process relating different equations can be seen as a `square-root of a geometry,' unless geometry in this context refers to the spin representation. It does not change the nature of the space in which the equations are defined, but it does change the kind of fields that the equation deals with.
A: 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 ;-).
A: This is not a complete answer but my idea on that topic. Pauli-Matrizes (SU(2)) obey the same commutation relations as space rotations (SO(3)). This is because SU(2) is the universal cover of SO(3). When using irreducible representations of a given group SO(n) you essential throw away the specific group structure and restrict yourself to the commutator structure only. That is, these representations are the ones of the universal cover. This is ok for (quantum) physics, because our theories are local theories and the local structure of the former group is the same as the one of the universal cover.
Now comes the key point. SO(3) is both a subgroup of the Galilean group and the Poincaré group. So when using the projective representations (as pointed out by one of the comments), both of them will lead to SU(2), and thus spin, being "part" of the projective unitary group structure.
Spin is a generalization of space-rotations and as such "part" of the projective representations of the Galilean group and the Poincaré group.
I don't have any idea on that Wick-rotation part of your question, though.
A: LATE UPDATE (at Vote = -2), re Issam’s Question 1, and generally:
AFAIK, no-one’s referred directly, here, to linearising Schrödinger’s equation (as distinct from linearising the K-GE). 
It’s the subject of Section 4.2 of Chapter 4 “Pauli Spinors” in the English translation of Jean Hladik’s small gem of a text “Spinors in Physics” (Springer, 1999) , which provides some concise summaries, scanned into here. 
1)  His intro on p. 100, just prior to Section 4.2 In Ch. 4: 
It was Pauli's theory, and then Darwin's—which sought to introduce the electron's magnetism in a way conforming to the theory of relativity by defining four functions representing the components of a space-time vector—which inspired Dirac to invent his theory of the relativistic electron. Dirac, examining the previous relativistic theories, was led to a new hypothesis, that the equations controlling the evolution of the components ψi of the wave function must be first order with respect to the four variables z, y, z, t, although the relativistic equations generalizing Schrödinger's equation were second order in these variables. [This para included for context only]
The idea of linearizing Schrödinger's equation itself was developed next, inspiring Dirac's works, and allowed the Pauli equations to be recovered in 
which the spin was then introduced automatically. The existence of spin is consequently not a purely relativistic effect, but becomes a consequence of the linearization of the wave equations. We are going to establish these linearized equations from the paper by Levy-Leblond (1967). [My bold and italic]
2)  Section 4.2 concludes (sorry - I haven’t the time to scan and correct 5 pp. of maths)
Thus we obtain the Pauli equation in which the term (eħ/2mc)σ.B appears, and which represents the interaction energy of the magnetic field with the intrinsic magnetic moment of the electron. Although Pauli had added this term in Schrödinger's equation in such a way as to make the theoretical and experimental results agree, we see that the spin is here introduced automatically as a consequence of the postulate of linearization of the wave equation. Moreover, this latter theory gives the correct value of the intrinsic magnetic moment of the electron. [END OF QUOTES]
Just to be clear: Hladik deals with the Dirac equation much later, in Chapter 7, so the subject of Section 4.2 is not Dirac’s work as such. [END OF LATE UPDATE]
INITIAL ANSWER (following some earlier edits):
This is an attempt to answer the OP’s Question 4 [the first Answer, from when the Question was "Spin and equations linearization?"], in the light of: 
a) @ACuriousMind’s early Comment on the vagueness of Issam’s Question 4 (above)
b) Michael Atiyah’s famous (if somewhat gnomic) statement about spinors and “the square root of geometry” (see e.g. JamalS’s answer to physics.stackexchange.com/questions/141995/how-should-i-think-about-the-dirac-equation)
and
c) the absence, so far, of a clarification of his Q. 4 from Issam himself:
Maybe Issam is simply using natural language rather loosely. Let’s examine it more closely…
Might “linearization” in Q. 4 refer mainly to Dirac’s development of his equation by forcing it into first order w.r.t. space and time, resulting in spinor solutions, and perhaps more generally to ‘Dirac Operators’? If so, then fair enough.
For Wick Rotation, on the other hand, Q. 4 might be referring to the distance measure ds in Minkowski space as being “linearised” by taking the square root of (dτ^2+dx^2+dy^2+dz^2) when it’s converted to a Euclidean form by setting it  = τ  (multiplication by the complex i acts as a rotation, hence Wick Rotation). While one might describe that metric expression/quadratic form, very loosely, as “linearised” (though the individual component terms are still squared) nevertheless, the null cone in Minkowski space can be parametrised by spinors. 
So spinors appear to provide at least a suggestive link between Wick Rotation and Issam’s use of Linearization. 
And further to the earlier suggestion that more of use might be found via: physics.stackexchange.com/questions/21261/wick-rotation-and-spinors, @udrv has latterly suggested (Comment below) that Issam may find the answer sought, on Wick rotations of Dirac fields, via arxiv.org/abs/hep-th/9611043  and the review mentioned therein of older literature: arxiv.org/abs/hep-th/9608174 (referenced by Qmechanic in his 1st Comment in .../21261/...).
I hope this helps to stimulate further discussion on the fundamentals of spin, it's descriptions, and its origins. 
