\section{Introduction}
The first positive measurements of spin polarization of $\Lambda$ hyperons made lately by the STAR Collaboration \cite{STAR:2017ckg,Adam:2018ivw,Niida:2018hfw,Adam:2019srw} have revived the interest in the studies of the relation between the vorticity of the matter produced in relativistic heavy-ion collisions and the averange spin polarization of particles produced in these processes \cite{Becattini:2009wh,Becattini:2013fla,Montenegro:2017rbu,Montenegro:2017lvf,Becattini:2018duy,Boldizsar:2018akg,Prokhorov:2018bql,Yang:2018lew,Florkowski:2019voj,Weickgenannt:2019dks,Hattori:2019lfp,Ambrus:2019ayb,Sheng:2019kmk,Prokhorov:2019cik,Ivanov:2019wzg,Hattori:2019ahi,Xie:2019jun,Liu:2019krs,Wu:2019eyi,Becattini:2019ntv,Zhang:2019xya,Li:2019qkf,Florkowski:2019gio,Prokhorov:2019yft,Fukushima:2020qta,Liu:2020ymh,Bhadury:2020puc,Tabatabaee:2020efb,Liu:2020flb,Yang:2020hri,Deng:2020ygd,Taya:2020sej}; for a recent review see \cite{Becattini:2020ngo}. Recently, it has been shown that the thermal-based models~\cite{Becattini:2016gvu,Karpenko:2016jyx,Li:2017slc,Xie:2017upb} which correctly describe the global polarization unfortunately are not able to explain the differential observables \cite{Adam:2019srw}. These models are based on the assumption that the spin polarization of particles emitted at freeze-out is entirely determined by the quantity known as thermal vorticity~\cite{Becattini:2007sr,Becattini:2013fla} and do not include the possibility of its independent dynamical evolution, which may take place during the fluid expansion. In this work,  following ideas put forward in Refs.~\cite{Florkowski:2017ruc,Florkowski:2017dyn,Florkowski:2018myy,Florkowski:2018fap,Florkowski:2018ahw,Florkowski:2019qdp}, we study such possibility within the framework of relativistic perfect-fluid hydrodynamics with spin.
\section{Equilibrium distribution functions}
Relativistic fluid dynamics may be derived from the underlying kinetic theory assuming that the distribution function describing the equilibrium state of the system is known \cite{Florkowski:2017olj}. Herein, following works by Becattini et al.  \cite{Becattini:2013fla},  we assume that the local equilibrium state of the relativistic system of particles ($+$) and antiparticles ($-$) with spin $\onehalf$ and mass $m$ is described by the following phase-space distribution functions (spin density matrices)
\bel{fplusrsxp}
f^+_{rs}(x,p) =  {\bar u}_r(p) X^+ u_s(p), \qquad
f^-_{rs}(x,p) = - {\bar v}_s(p) X^- v_r(p),
\end{eqnarray}
where $x$ is the space-time position and $p$ is the four-momentum, and $u_r(p)$ and $v_r(p)$ are Dirac bispinors ($r,s = 1,2$) with the normalization ${\bar u}_r(p) u_s(p)=\,\delta_{rs}$ and ${\bar v}_r(p) v_s(p)=-\,\delta_{rs}$.

The matrices $X^{\pm}$ have the form of generalized relativistic Boltzmann distributions
\bel{XpmM}
X^{\pm} =  \exp\left[\pm \xi(x) - \beta_\mu(x) p^\mu  \pm \f{1}{2} \omega_{\mu\nu}(x)  {\Sigma}^{\mu\nu} \right], \nonumber
\end{eqnarray}
where $\beta^\mu \equiv U^\mu/T$ and $\xi \equiv \mu/T$, with $T$, $\mu$ and $U^\mu$ denoting the temperature, baryon chemical potential and four-velocity,  respectively. The quantity $\omega_{\mu\nu}$ is the spin polarization tensor satisfying $\omega^{\mu\nu}=-\omega^{\nu\mu}$ and ${\Sigma}^{\mu\nu}  \equiv  \f{i}{4} [\gamma^\mu,\gamma^\nu]$ is the spin operator.

Employing expressions derived in Ref.~\cite{DeGroot:1980dk} and using definitions \rfn{fplusrsxp} we can determine the corresponding equilibrium Wigner functions
\begin{eqnarray}
{\cal W}^\pm_{\rm eq}(x,k) &=& \frac{e^{\pm \xi}}{4 m}  \int dP
\,e^{-\beta \cdot p }\,\, \delta^{(4)}(k \mp p) \label{eq:Weqpxk2} \\
&& \times  \left[2m (m \pm \slashed{p}) \cosh(\zeta) \pm \f{\sinh(\zeta)}{2\zeta}  \, \omega_{\mu\nu} \,(\slashed{p} \pm m) {\Sigma}^{\mu\nu} (\slashed{p} \pm m) \right], \nonumber
\end{eqnarray}
where $k$ is the off-mass-shell four-momentum of particles,
$dP = d^3p/((2 \pi )^3 E_p)$  with $E_p = \sqrt{m^2 + {\boldsymbol p}^2}$ being the on-shell particle energy, and $\zeta =  \f{1}{2 \sqrt{2}} \sqrt{ \omega_{\mu\nu} \omega^{\mu\nu}}$.

It is convenient to consider the Clifford-algebra expansion of the Wigner function \rfn{eq:Weqpxk2}
\begin{eqnarray}
{\cal W}^{\pm}_{\rm eq}(x,k) &=& \f{1}{4} \left[ {\cal F}^{\pm}_{\rm eq}(x,k) + i \gamma_5 {\cal P}^{\pm}_{\rm eq}(x,k) + \gamma^\mu {\cal V}^\pm_{{\rm eq}, \mu}(x,k) \right. \nonumber \\
&& \left.  \hspace{1cm} + \gamma_5 \gamma^\mu {\cal A}^\pm_{{\rm eq}, \mu}(x,k)
+ {\Sigma}^{\mu\nu} {\cal S}^\pm_{{\rm eq}, \mu \nu}(x,k) \right], \nonumber
\label{eq:wig_expansion}
\end{eqnarray}
where the coefficient functions $\mathcal{X} \in\left\{\mathcal{F}, \mathcal{P}, \mathcal{V}_{\mu}, \mathcal{A}_{\mu}, \mathcal{S}_{\nu \mu}\right\}$ can be extracted by calculating the trace of ${\cal W}^{\pm}_{\rm eq}(x,k)$ multiplied first by: $\left\{\mathbf{1},-i \gamma_{5}, \gamma_{\mu}, \gamma_{\mu} \gamma_{5}, 2 \Sigma_{\mu \nu}\right\}$.
\section{Kinetic equations}
General Wigner function satisfies the kinetic equation
\bel{eq:eqforWC}
\left(\gamma_\mu K^\mu - m \right) {\cal W}(x,k) = C[{\cal W}(x,k)],
\end{eqnarray}
where the differential operator reads $K^\mu = k^\mu + \frac{i \hbar}{2} \,\partial^\mu$.
In the case of global equilibrium, the Wigner function satisfies exactly \rf{eq:eqforWC} with $C[{\cal W}(x,k)] = 0$. The usual treatment of \rf{eq:eqforWC}  is to consider the semi-classical expansion of the coefficient functions
\bel{eq:semi}
\mathcal{X}=\mathcal{X}^{(0)}+\hbar \mathcal{X}^{(1)}+\hbar^{2} \mathcal{X}^{(2)}+\cdots. \nonumber
\end{eqnarray}
The analysis of \rf{eq:eqforWC} up to the next-to-leading order in $\hbar$ yields the following kinetic equations for the two independent coefficients $\mathcal{F}_{\mathrm{eq}}$ and $\mathcal{A}_{\mathrm{eq}}^{\nu}$,
\bel{eq:ke}
k^{\mu} \partial_{\mu} \mathcal{F}_{\mathrm{eq}}(x, k)=0, \quad k^{\mu} \partial_{\mu} \mathcal{A}_{\mathrm{eq}}^{\nu}(x, k)=0, \quad k_{\nu} \mathcal{A}_{\mathrm{eq}}^{\nu}(x, k)=0.
\end{eqnarray}
In global equilibrium  \rfm{eq:ke} are exactly fulfilled which results in the conditions that $\beta_{\mu}$ is a Killing vector, and $\xi$ and $\omega_{\mu \nu}$ are constant, however $\omega_{\mu \nu}$ does not have to be equal to  thermal vorticity $\varpi_{\mu \nu} = -\frac{1}{2} \left(\partial_\mu \beta_\nu - \partial_\nu \beta_\mu \right) = \hbox{const}$.
\section{Hydrodynamic equations}
In local equilibrium \rfm{eq:ke} are not satisfied exactly. In this case, we
adopt the standard treatment \cite{Denicol:2012cn}, namely,  by allowing for $x$ dependence of the $\beta$, $\xi$ and $\omega$, we require that only certain moments in momentum space of the kinetic equations \rfn{eq:ke} are satisfied. This method leads to equations expressing conservation laws for charge, energy, linear momentum and spin\cite{Florkowski:2018ahw}
\begin{eqnarray}
\quad\partial_\mu N^\mu = 0, \label{Ncons} \\
\partial_\mu T^{\mu\nu}_{\rm GLW} = 0, \label{Tcons}\\
\partial_\lambda  S_{\rm GLW}^{\lambda, \alpha \beta} =0,
\label{Scons}
\end{eqnarray}
where the baryon current, the energy-momentum tensor, and the spin tensor are given by the de Groot - van Leeuwen - van Weert (GLW)~\cite{DeGroot:1980dk}  expressions
\begin{eqnarray}
N^\alpha &=& n U^\alpha,\\
 T^{\alpha\beta}_{\rm GLW} &=& (\varepsilon + P ) U^\alpha U^\beta - P g^{\alpha\beta},\\
 S^{\alpha , \beta \gamma }_{\rm GLW}
&=& \cosh(\xi) \LS n_{(0)} U^\alpha \omega^{\beta\gamma}  +  {\cal A}_{(0)} \, U^\alpha U^\d U^{[\beta} \omega^{\gamma]}_{\HP\d}  \right. \\
&& +\left. \, {\cal B}_{(0)} \, \Big(
U^{[\beta} \Delta^{\alpha\d} \omega^{\gamma]}_{\HP\d}
+ U^\alpha \Delta^{\d[\beta} \omega^{\gamma]}_{\HP\d}
+ U^\d \Delta^{\alpha[\beta} \omega^{\gamma]}_{\HP\d}\Big) \right],
\end{eqnarray}
where $\Delta^{\mu\nu} =g^{\mu\nu} - U^\mu U^\nu$ is the projector on the space orthogonal to $U$.

In the leading order in the polarization tensor the energy density $\varepsilon$, the pressure $P$, and the baryon density $n$ are given by the formulas
\begin{eqnarray}
n &=& 4 \, \sinh(\xi)\, n_{(0)}(T), \label{n0small} \\
\varepsilon &=& 4 \, \cosh(\xi) \, \varepsilon_{(0)}(T),\label{e0small}\\
P &=& 4 \, \cosh(\xi) \, P_{(0)}(T), \label{P0small}
\end{eqnarray}
where we defined the auxiliary quantities describing thermodynamic properties of the system of spin-less and neutral massive Boltzmann particles\cite{Florkowski:2010zz}
\begin{eqnarray}
n_{(0)}(T) &=&  \f{ T^3}{2\pi^2}\,  \hat{m}^2 K_2\left( \hat{m}\right), \label{n0c}\\
\varepsilon_{(0)}(T) &=& \f{ T^4 }{2\pi^2}  \, \hat{m} ^2
 \Big[ 3 K_{2}\left( \hat{m} \right) + \hat{m}  K_{1} \left( \hat{m}  \right) \Big],  \label{e0c}\\
 %
P_{(0)}(T) &=& T \, n_{(0)}(T) . \label{P0c}
\end{eqnarray}
 %
 The quantities ${\cal B}_{(0)} $ and ${\cal A}_{(0)}$ are defined as follows
 %
\begin{eqnarray}
{\cal B}_{(0)} =-\frac{2}{\hat{m}^2} s_{(0)}(T) ,\qquad
{\cal A}_{(0)}  = -3{\cal B}_{(0)} +2 n_{(0)}(T)
\end{eqnarray}
with $s_{(0)} =   \LR\varepsilon_{(0)}+P_{(0)}\right) / T$ being the entropy density and $\hat{m}=m/T$.
\section{Bjorken expansion}
Similarly to the Faraday tensor the polarization tensor $\omega_{\mu\nu}$ may be decomposed into electric-like ($\kappa$) and magnetic-like ($\omega$) components
\begin{eqnarray}
\omega_{\mu\nu} &=& \kappa_\mu U_\nu - \kappa_\nu U_\mu + \epsilon_{\mu\nu\alpha\beta} U^\alpha \omega^{\beta},\label{spinpol1}
\end{eqnarray}
where the four-vectors $\kappa$ and $\omega$ satisfy the conditions
 %
\begin{eqnarray}
\kappa\cdot U = 0,\qquad \omega \cdot U = 0.
\label{orth}
\end{eqnarray}
In the case of transversely homogeneous systems undergoing boost-invariant expansion in the longitudinal direction, also known as the Bjorken flow \cite{Bjorken:1982qr}, it is convenient to introduce the following four-vector basis
\begin{eqnarray}
U^\alpha &=& \frac{1}{\tau}\LR t,0,0,z \right) = \LR \cosh(\eta), 0,0, \sinh(\eta) \right), \nonumber \\
X^\alpha &=& \LR 0, 1,0, 0 \right),\nonumber\\
Y^\alpha &=& \LR 0, 0,1, 0 \right), \nonumber\\
Z^\alpha &=& \frac{1}{\tau}\LR z,0,0,t \right) = \LR \sinh(\eta), 0,0, \cosh(\eta) \right),
\label{BIbasis}
\end{eqnarray}
where $\tau = \sqrt{t^2-z^2}$ is the longitudinal proper time and $\eta = \half \ln((t+z)/(t-z))$ is the space-time rapidity.

The basis \rfn{BIbasis} satisfies the conditions
\begin{eqnarray}
 &&U \cdot U = 1,\nonumber\\
X \cdot X &=&   Y \cdot Y \,\,=\,\, Z \cdot Z \,\,=\,\, -1,\nonumber \\ \label{orthXYZ}
X \cdot U  \,\,&=& Y \cdot U\,\, \,\,=\,\, Z \cdot U \,\,=\,\, 0,   \\
X \cdot Y  &=&  Y \cdot Z \,\,=\,\, Z \cdot X \,\,=\,\, 0.  \nonumber
\end{eqnarray}
Using \rfm{orth} and \rfm{orthXYZ}, one can decompose the vectors $\kappa^{\mu}$ and $\omega^{\mu}$ as follows
\begin{eqnarray}
\kappa^\alpha &=&  C_{\kappa X} X^\alpha + C_{\kappa Y} Y^\alpha + C_{\kappa Z} Z^\alpha,  \nonumber\\
\omega^\alpha &=&  C_{\omega X} X^\alpha + C_{\omega Y} Y^\alpha + C_{\omega Z} Z^\alpha, \label{decom}
\end{eqnarray}
where the coefficients ${C}_{\kappa X}$, ${C}_{\kappa Y}$, ${C}_{\kappa Z}$, ${C}_{\omega X}$, ${C}_{\omega Y}$, and ${C}_{\omega Z}$ are scalar functions of proper time solely.

Using \rfm{decom} in \rf{Scons} and projecting the latter on $U_\mu X_\nu$, $U_\mu Y_\nu$, $U_\mu Z_\nu$, $X_\mu Y_\nu$, $X_\mu Z_\nu$ and $Y_\mu Z_\nu$, we obtain the following six evolution equations
\begin{equation}
{\rm diag}\LR
\cal{L}, \cal{L}, \cal{L}, \cal{P}, \cal{P}, \cal{P}\right) \,\,
\Dot{{\boldsymbol C}} ={\rm diag}\LR
{\cal{Q}}_1, {\cal{Q}}_1, {\cal{Q}}_2, {\cal{R}}_1, {\cal{R}}_1, {\cal{R}}_2 \right)\,\,
{\boldsymbol C}, \label{cs}
\end{equation}
where ${\boldsymbol C} = \LR C_{\kappa X}, C_{\kappa Y}, C_{\kappa Z},  C_{\omega X}, C_{\omega Y}, C_{\omega Z} \right)$, $\dot{(\dots)} \equiv U \cdot \partial = \partial_\tau$ and
\begin{eqnarray}
{\cal L}(\tau)&=&{\cal A}_1-\frac{1}{2}{\cal A}_2-{\cal A}_3,\nonumber\\
{\cal P}(\tau)&=&{\cal A}_1,\nonumber\\
{\cal{Q}}_1(\tau)&=&-\left[\dot{{\cal L}}+\frac{1}{\tau}\left( {\cal L}+ \frac{1}{2}{\cal A}_3\right)\right],\nonumber\\
 {\cal{Q}}_2(\tau)&=&-\left(\dot{{\cal L}}+\frac{{\cal L}}{\tau}   \right),\nonumber\\
  {\cal{R}}_1(\tau)&=&-\left[\Dot{\cal P}+\frac{1}{\tau}\left({\cal P} -\frac{1}{2} {\cal A}_3 \right) \right],\nonumber\\
 {\cal{R}}_2(\tau)&=&-\left(\Dot{{\cal P}} +\frac{{\cal P}}{\tau}\right).\nonumber
 \label{LPQR}
 \end{eqnarray}
 %
 with
 \begin{eqnarray}
{\cal A}_1 &=& \cosh(\xi) \LR n_{(0)} -  {\cal B}_{(0)} \right) \label{A1} ,\nonumber\\
{\cal A}_2 &=& \cosh(\xi) \LR {\cal A}_{(0)} - 3{\cal B}_{(0)} \right) \nonumber \label{A2} , \\
{\cal A}_3 &=& \cosh(\xi)\, {\cal B}_{(0)}\label{A3},\nonumber
\end{eqnarray}

From \rfm{cs} we observe that in the case of Bjorken expansion the ${C}$ coefficients evolve independently. Due to the rotational symmetry in the transverse plane the functions ${C}_{\kappa X}$ and ${C}_{\kappa Y}$ (as well as ${C}_{\omega X}$ and ${C}_{\omega Y}$) obey the same differential equations.

Employing the Bjorken symmetry,  the conservation of the charge current \rfn{Ncons} can be written as
\begin{eqnarray}
\frac{dn}{d\tau}+\frac{n}{\tau}=0\label{ns}
\end{eqnarray}
while the conservation of energy and linear momentum \rfn{Tcons}  (projected on $U$) gives
\begin{eqnarray}
\frac{d\varepsilon}{d\tau}+\frac{(\varepsilon+P)}{\tau}=0.\label{ts}
\end{eqnarray}
\section{Spin polarization of particles at freeze-out}
The information about the space-time evolution of the spin polarization tensor may be used to determine the average spin polarization per particle which is defined as follows~\cite{Florkowski:2018ahw}
\begin{eqnarray}
\langle\pi_{\mu}\rangle=E_p\frac{d\Pi _{\mu }(p)}{d^3 p}/E_p\frac{d{\cal{N}}(p)}{d^3 p}, \label{pi}
\end{eqnarray}
where $E_p\frac{d\Pi _{\mu }(p)}{d^3 p}$ is the total value (integrated over freeze out hypersuface $\Delta \Sigma_{\lambda }$) of the Pauli-Luba\'nski vector,
\begin{eqnarray}
E_p\frac{d\Pi _{\mu }(p)}{d^3 p} = -\f{ \cosh(\xi)}{(2 \pi )^3 m}
\int
\Delta \Sigma _{\lambda } p^{\lambda } \,
e^{-\beta \cdot p} \,
\tilde{\omega }_{\mu \beta }p^{\beta },  \nonumber
\end{eqnarray}
and
\begin{eqnarray}
E_p\frac{d{\cal{N}}(p)}{d^3 p}&=&
\f{4 \cosh(\xi)}{(2 \pi )^3}
\int
\Delta \Sigma _{\lambda } p^{\lambda }
\,
e^{-\beta \cdot p}, \nonumber
\end{eqnarray}
is the momentum density of particles and antiparticles.

The polarization vector $\langle\pi^{\star}_{\mu}\rangle$ in the particle rest rame is obtained by performing the canonical boost \cite{Leader:2001gr}  of \rfn{pi}
\begin{eqnarray}
\langle\pi^{\star}_{\mu}\rangle=-\frac{1}{8m }\left[\begin{array}{c}
0 \\ \\
\left(\frac{p_x\sh y_p }{b}\right) a_1 +
  \LR\frac{ \chi \,p_x \ch y_p  }{b} \right) a_2\!+\!2 C_{\kappa Z} p_y  \!-\!\chi C_{\omega X}{m}_T \\ \\
 %
\left(\frac{p_y\sh y_p }{b}\right)a_1+   \LR\frac{ \chi \,p_y \ch y_p  }{b} \right) a_2\!-\!2 C_{\kappa Z} p_x \!-\!\chi C_{\omega Y}{m}_T \\ \\
 %
 -\left(\frac{m\ch y_p+m_T}{b}\right)a_1
-\LR\frac{\chi \,m\,\sh y_p}{b}\right) a_2 \nonumber\\
\end{array}
\right],\\
\label{PLVPPLPRF}
\end{eqnarray}
where $a_1=\chi\left(C_{\kappa X} p_y-C_{\kappa Y} p_x\right)+2 C_{\omega Z} m_T$, $a_2=C_{\omega X} p_x+C_{\omega Y} p_y$, $b=m_T \ch y_p+m$, and $\chi=\left( K_{0}\left( \hat{m}_T \right)+K_{2}\left( \hat{m}_T \right)\right)/K_{1}\left( \hat{m}_T \right)$ with  $\hat{m}_T=m_T/T$.  Above we used the following parametrization of the four-momentum  $p^\lambda = \left( m_T\ch y_p,p_x,p_y,m_T\sh y_p \right)$.

\begin{figure}
\centering
 \includegraphics[width=0.6\textwidth]{Tmu.pdf}
\caption{Proper-time dependence of the temperature scaled by its initial value (solid black line) and the ratio of baryon chemical potential over temperature rescaled by the initial ratio (dotted blue line).}
\label{fig:Tmu}
\end{figure}

\begin{figure}
\centering
 \includegraphics[width=0.6\textwidth]{C.pdf}
\caption{Proper-time dependence of the coefficients $C_{\kappa X}$ (solid black line), $C_{\kappa Z}$ (dashed-dotted blue line), $C_{\omega X}$ (dotted red line) and $C_{\omega Z}$ (dashed green line).}
\label{fig:C}
\end{figure}
\section{Results}
In this section we present the results obtained by solving numerically the differential equations \rfn{cs}, \rfn{ns}, and \rfn{ts}. We initialize the system at the proper time $\tau_0 = 1$ fm with initial temperature $T_0=T(\tau_0)=150$ MeV and the baryon chemical potential $\mu_0=\mu(\tau_0)=800$ MeV. We assume that the system consists of $\Lambda$ particles with mass $m=1116$ MeV. In Fig.~\ref{fig:Tmu}, we show the proper-time dependence of the (properly scaled) temperature and baryon chemical. We reproduce the well known results that the temperature of such a system decreases with proper-time while the ratio of chemical potential and temperature increases. In Fig.~\ref{fig:C}, we show the proper time dependence of the $C$ coefficients that describe the evolution of the spin polarization.

The knowledge of the evolution of thermodynamic parameters and $C$ coefficients allows us to calculate the components of the particle-rest-frame mean polarization vector $\langle\pi^{\star}_{\mu}\rangle$ at freeze-out as functions of particle three-momentum, see Fig.~\ref{fig:polarization1}.
We observe that the component $\langle\pi^{\star}_{y}\rangle$ is negative, which reflects the initial spin polarization of the system.  Due to the Bjorken symmetry the longitudinal component ($\langle\pi^{\star}_{z}\rangle$) is vanishing which does not agree with the characteristic quadrupole structure of the longitudinal polarization observed in the experiment. One the other hand we observe $\langle\pi^{\star}_{x}\rangle$ exhibits quadrupole structure. Clearly, we observe that the Bjorken symmetry is too restrictive to address the experimental measurements correctly.

\begin{figure}
\centering
 \includegraphics[width=0.32\textwidth]{px.pdf}
\includegraphics[width=0.32\textwidth]{py.pdf}
 \includegraphics[width=0.32\textwidth]{pz.pdf}
\caption{Components of the particle-rest-frame mean polarization three-vector of $\Lambda$ particles obtained with the initial conditions $\mu_0=800$~MeV,
$T_0=155$~MeV, $C_{\kappa, 0}=(0,0,0)$ and $ C_{\omega, 0}=(0,0.1,0)$ for $y_p=0$.}
\label{fig:polarization1}
\end{figure}
\section{Summary}
In this work we briefly reviewed the basic  ingredients  of the recently formulated  approach of relativistic perfect-fluid hydrodynamics with spin. Using the kinetic theory definitions for the local equilibrium distribution functions we derived  the evolution equations governing the system's expansion. Assuming the Bjorken flow of the matter we studied numerically the spin polarization dynamics of the system.
We have shown that the coefficient functions characterizing the spin polarization evolve independently. We have used these results to determine the spin polarization of particles at the freeze-out. We have shown that  within the simple Bjorken setup the characteristic features observed in the experiment can not be properly reproduced.
\section*{Acknowledgments}
Supported in part by the Polish National Science Center Grants No.   2016/23/B/ST2/00717 and No. 2018/30/E/ST2/00432.


\section{Construction of 3D Spin-Orbit coupling}

The 3D SO coupling is constructed with 3D optical Raman lattices, as outlined in Fig.\ref{fig1}(a).
The spin is defined from the magnetic sublevels of $^{87}$Rb atoms $\mid\uparrow\rangle=\left|1,-1\right\rangle$ and $\mid\downarrow\rangle=\left|1,0\right\rangle$.
A bias magnetic field $\boldsymbol{B}=14.5$G provides the quantization in $y$-axis and the Zeeman splitting of $10.2\text{MHz}$.
The $|1,+1\rangle$ state is effectively excluded.
Three retro-reflective laser beams $\boldsymbol{E}_{i}(i=x,y,z)$ with wavelength $\lambda=787$nm ($k_{0}=2\pi/\lambda$ and recoil energy $E_{\rm{r}}=\hbar^{2}k_{0}^{2}/2m$) shine on the atoms in $i$-th direction, each having two orthogonal polarization components $E_{ij}(j=x,y,z)$.
$E_{xy}$, $E_{yx}$ have frequency $\omega_{1}$, $E_{xz}$, $E_{yz}$ have frequency $\omega_{2}$, and $E_{zy}$, $E_{zx}$ have frequency $\omega_{3}$.
We set $\omega_3-\omega_1\approx 2\pi\times 10.2\text{MHz}$, matching the Zeeman splitting of $\mid\uparrow\rangle$ and $\mid\downarrow\rangle$, while $\omega_1-\omega_2=2\pi\times 200\text{kHz}$.
Two $\lambda/4$ waveplates with optical axis aligned along $\hat{z}$ are placed in front of the retro-reflectors,
producing $\pi$ phase shift between $E_{xz(yz)}$ and $E_{xy(yx)}$.

Our realization is based on configuration-tunable lattices. $E_{xy}$ and $E_{yx}$ form a square lattice, together with $E_{xz}$ and $E_{yz}$, deforming into the checkerboard lattice in x-y plane $\mathcal{V}_{2D}(x,y)=V_{2D}\cos{k_{0}x}\cos{k_{0}y}$ upon adjusting relative phase between $\boldsymbol{E}_{x}$ and $\boldsymbol{E}_{y}$ to $\phi=0$ (See Methods).
The 2D lattice potential $\mathcal{V}_{2D}(x,y)$ can be rewritten as $\mathcal{V}_{2D}(u,v)=V_{2D}\left(\cos^{2}{k_1u}+\cos^{2}{k_1v}\right)$ along $\hat u$ and $\hat v$ directions, with $\hat{u}=(\hat{x}+\hat{y})/\sqrt{2}$, $\hat{v}=(\hat{x}-\hat{y})/\sqrt{2}$ and $k_1=k_{0}/\sqrt{2}$.
The density profile of $\mathcal{V}_{2D}$ is shown in Fig.\ref{fig1}(b).
The total 3D lattice potential $\mathcal{V}_{Latt}\left(u,v,z\right)=\mathcal{V}_{2D}(u,v)+\mathcal{V}_{z}(z)$, where
$\mathcal{V}_{z}(z)=V_{z}\cos^{2}{k_{0}z}$ along $z$ direction, with depth $V_{z}\propto\left|E_{zx}\right|^{2}+\left|E_{zy}\right|^{2}$.

Raman couplings are generated by beam pairs $\left(E_{zx},E_{xy}\right)$ and $\left(E_{zy},E_{yx}\right)$,
with $\Omega_{1}=\Omega_{zx}\cos{k_{0}z}\sin{k_{0}x}$ and $\Omega_{2}=\Omega_{zy}\cos{k_{0}z}\sin{k_{0}y}$, where Rabi frequencies $\Omega_{zx}\propto E^*_{zx}E_{xy}$ and $\Omega_{zy}\propto E^*_{zy}E_{yx}$ are set as $|\Omega_{zx}|=|\Omega_{zy}|=\Omega_{0}$.
The large frequency difference $\omega_2-\omega_1$ prevents the couplings between $(E_{zy}, E_{xz})$ and $(E_{zy}, E_{yz})$.
For the circular polarized $\boldsymbol{E}_{z}$, $\Omega_{1}$ and $\Omega_{2}$ have $\pi/2$ phase shift~\cite{Theory_Liu}.
Then the Raman couplings, having 2D forms in $(x,y,z)$ coordinates, take 3D forms in the deformed checkerboard lattice bases and $\Omega_{R}\left(u,v,z\right)=\Omega_{u}\left(u,v,z\right)\sigma_{x}+\Omega_{v}\left(u,v,z\right)\sigma_{y}$, where $\Omega_{u}\left(u,v,z\right)=\sqrt{2}\Omega_{0}\cos{k_{0}z}\sin{k_1u}\cos{k_1v}$ and $\Omega_{v}=\Omega_{u}\left(u\leftrightarrow v,z\right)$.
As shown in Fig.\ref{fig1}(b), $\Omega_{u}$ ($\Omega_{v}$) is antisymmetric along $\hat{u}$ ($\hat{v}$) and symmetric along $\hat{z}$, leading to spin-flipped hopping in $\hat{u}$ ($\hat{v}$) direction.
We finally reach the Hamiltonian as
\begin{equation}\label{3DSOC}
    H=\frac{\hbar^2\boldsymbol{k}^{2}}{2m}+\mathcal{V}_{Latt}(u,v,z)+\Omega_{u}(u,v,z)\sigma_{x}+\Omega_{v}(u,v,z)\sigma_{y}+\frac{\delta}{2}\sigma_{z},
\end{equation}
where $\hbar\boldsymbol{k}$ is momentum, $m$ is the mass of an atom and $\delta$ is the two-photon Raman detuning.

With the tight-binding (TB) approximation, Eq.(\ref{3DSOC}) in Bloch momentum space reads \cite{Theory_Liu,WangPRA2016}
\begin{equation}\label{TBH}
  \begin{split}
    H_{\text{TB}}(\boldsymbol{q})&=\boldsymbol{h}(\boldsymbol{q})\cdot\boldsymbol{\sigma}\\
    &=2t_\text{so}(\sin{q_{u}}\sigma_{x}+\sin{q_{v}}\sigma_{y})+\left[m_{z}-2t_{z}\cos{q_{z}}-2t_{1}\left(\cos{q_{u}}+\cos{q_{v}}\right)\right]\sigma_{z},
  \end{split}
\end{equation}
where $\boldsymbol{q}=(q_u, q_v, q_z)$ is dimensionless quasi-momenta, $t_{z,1}$ ($t_\text{so}$) denote the spin-conserved (spin-flipped) hopping coefficients, and $m_{z}=\delta/2$.
For each fixed $q_z$, $H_{\text{TB}}$ gives to a 2D QAH model in $u$-$v$ plane~\cite{Theory_Liu,Wu_Science,Baozong,Sun_2018A}, with topology modulated by both $m_Z$ and $q_z$. The Weyl points emerge when the Chern number of QAH models changes versus $q_z$.

The experiment starts from an optical trapped $^{87}\textrm{Rb}$ Bose-Einstein condensate (BEC) with $2.0\times10^{5}$ atoms in $\mid\uparrow\rangle=\left|1,-1\right\rangle$ state.
As the 3D Raman lattice beams are adiabatically ramped up in 100ms, the BEC is loaded into the ground state of Eq.(\ref{3DSOC}).
The probe beams are shining along $\hat{z}$ and $\hat{y}$ (shown in Fig.\ref{fig1}(a)) after the atoms are freely released for 22ms,
to take the spin-resolved time-of-flight (ToF) 2D images on CCD-Z and CCD-Y, as shown in Fig.\ref{fig1}(c).
The $\mid\uparrow\rangle$ clouds are from the diffraction of the 3D optical lattices.
On CCD-Z, the major portion of the BEC is at momentum $\left(k_{x},k_{y}\right)=(0,0)$ and four small dots at momenta $\left(k_{x},k_{y}\right)=(\pm k_0,\pm k_0)$.
On CCD-Y, there are two diffractions are at $(k_{x}, k_{z})=(0, \pm 2k_0)$ and the overlaps of the four diffractions at $(k_x=k_0, k_y=\pm k_0, k_z=0)$, $(k_x=-k_0, k_y=\pm k_0, k_z=0)$.
The $\mid\downarrow\rangle$ clouds are from the Raman couplings.
On CCD-Z, they are at $\left(k_{x},k_{y}\right)=(\pm k_0,0)$ and $\left(k_{x},k_{y}\right)=(0,\pm k_0)$. On CCD-Y, they are at $\left(k_{x},k_{z}\right)=(k_0,\pm k_0)$, $\left(k_{x},k_{z}\right)=(0,\pm k_0)$ and $\left(k_{x},k_{z}\right)=(-k_0,\pm k_0)$.
The 3D momentum distributions of the $\mid\downarrow\rangle$ clouds, as illustrated in Fig.\ref{fig1}(d), confirm the successful realization of 3D SO coupling.

\section{Topological phase diagram and Weyl node measurement}

The Weyl node momentum $\boldsymbol{q}^{\text{W}}$ is determined via $\left|\boldsymbol{h}(\boldsymbol{q}^{\text{W}})\right|=0$.
Around these points, the linear expansion $H(\boldsymbol{q}^\text{W}+\delta\boldsymbol{q})=2t_\text{so}\delta q_{u}\sigma_{x}+2t_\text{so}\delta q_{v}\sigma_{y}+2t_{z}\sin{q_z^{\text{W}}}\delta q_{z}\sigma_{z}$ renders the Weyl Hamiltonian around the $\boldsymbol{q}^{\text{W}}$ node.
Fig.\ref{fig2}(a) is the calculated topological phase diagram of Eq.(\ref{3DSOC}).
The regions with number of Weyl nodes from 0 to 8 are presented by colors.
Importantly, there are phases with only two Weyl nodes, corresponding to the IWSM phase.
The insets of Fig.\ref{fig2}(a) show the Weyl nodes in the 3D Brillouin zone, which marked by $\oplus$ or $\ominus$ chirality.

To determine the positions of Weyl nodes in experiment, one can characterize the Chern number of 2D band structures for given $q_{z}$ plane.
In fact, Hamiltonian Eq.(\ref{3DSOC}) delineates a set of 2D QAH bands in $u$-$v$ planes, stacked along $q_{z}$.
The locations of Weyl nodes are characterized as the points where the topology of 2D QAH bands changes via scanning $q_{z}$~\cite{Theory_Liu}.

In the experiment, about $2.0\times10^{5}$ atoms with temperature around 150nK are adiabatically loaded into the lowest bands of Eq.(\ref{3DSOC}).
After a spin-resolved ToF imaging along $\hat{z}$ direction, we obtain the 2D momentum distribution ($q_u$ and $q_v$) of the atom in $\mid\uparrow\rangle$ ($n_{\uparrow}(\boldsymbol{q})$) and $\mid\downarrow\rangle$ ($n_{\downarrow}(\boldsymbol{q})$) on CCD-Z, with $q_{z}$ being integrated out.
Spin polarization at $\boldsymbol{q}$ is calculated by $P(\boldsymbol{q})=\left[n_{\uparrow}(\boldsymbol{q})-n_{\downarrow}(\boldsymbol{q})\right]/\left[n_{\uparrow}(\boldsymbol{q})+n_{\downarrow}(\boldsymbol{q})\right]$.
As shown in Ref.~\cite{Theory_Liu}, the Weyl band has an emergent magnetic group symmetry, for which we can reconstruct the 3D topological bands from a series of $q_z$-integrated 2D spin textures.
In particular, for a given Raman detuning $\delta_0$ in Eq.(\ref{3DSOC}), the observed spin texture in $q_u$-$q_v$ plane on CCD-Z, with $q_z$ being integrated out, is identical to its 2D spin texture on the $q_z=\pi/2$ plane of the 3D structure of spin distribution.
The spin textures of other $q_z$ planes with $\delta=\delta_0$ are equivalent to the 2D $q_z$-integrated spin texture with $\delta=\delta_0+\delta'$.
Therefore, by scanning $\delta'$, we can achieve a series of 2D spin textures to reconstruct the 3D spin distribution of given $\delta_0$. This is called Virtual Slicing Image~\cite{Gyu-Boong_Jo,Theory_Liu} (For details, see Methods).
Owing to our stable bias magnetic field~\cite{Magnetic}, $\delta'$ can be tuned precisely.
Fig.\ref{fig2}(b) shows the typical observed 2D spin textures at $\delta_0=-0.5E_{\rm{r}}$, where $q_{z}=0, 0.2\pi, 0.65\pi$ are equivalent to $\delta'=0.396E_{r},0.315E_r,0.176E_r$, respectively.
The spin texture for $q_{z}=0, 0.2\pi$ exhibit band inversion ring patterns~\cite{LZhang2018} corresponding to spin polarizations $P(\boldsymbol{q})=0$,~\cite{Theory_Liu, Wu_Science, Sun_2018A} marked by black circles.
Chern number for each of these 2D bands with fixed $q_{z}$ is determined by the product $\Theta$ of signs of spin polarizations $P(\Lambda_j,q_z)$ at four high symmetric momenta $\left\{\Lambda_{j}\right\}=\left\{\Gamma(0,0),X_{1}(0,\pi),X_{2}(\pi,0),M(\pi,\pi)\right\}$ \cite{C4_Liu,Wu_Science}, $\Theta(q_z)=\prod_j\text{sgn}(P(\Lambda_j,q_z))$.
Thereby scanning the $q_{z}$ planes through the whole 3D Brillouin zone, the topology of the 3D Weyl semimetal band can be fully determined.
Moreover, Weyl nodes emerge on the planes where Chern numbers change by integers.

The stack 2D spin textures in Bloch momentum space with $\delta_0=-0.5E_{\rm{r}}$ is shown in Fig.\ref{fig2}(c).
It is obvious that only two Weyl nodes are here, marked as $\oplus$ and $\ominus$. It corresponds to the IWSM regime.
Black dashed lines sketch the virtual slices of the fusiform band inversion surfaces.
When $q_{z}$ close to 0, spin polarizations near $\Gamma$ point $(0,0,q_{z})$ are negative and outside are positive.
Band inversion ring (black circles in Fig.\ref{fig2}(b)) encloses $\Gamma(0,0,q_{z})$, leads to Chern number $\mathcal{C}=1$ on this 2D plane~\cite{LZhang2018,Sun-2018B}.
As $q_{z}$ goes away from 0 plane, the band inversion rings gradually shrink. Eventually, they vanish at Weyl nodes where $\mathcal{C}$ jumps from 1 to 0.
When $q_{z}$ is close to $\pm\pi$, spin polarizations over the whole layer are positive, and $\mathcal{C}=0$.
From the stack of spin textures in Fig.\ref{fig2}(c), $P(\Lambda_j,q_z)$ are extracted and plotted in Fig.\ref{fig2}(d) as functions of $q_{z}$. $\Theta(q_z)$ is positive for $q_{z}\sim [-\pi,-0.3\pi]$ and $[0.3\pi,\pi]$, while negative for $q_{z}\sim[-0.3\pi,0.3\pi]$. Hence, $\mathcal{C}$ jumps from 0 to 1 at $q_{z}=\pm(0.3\pm0.03)\pi$, indicating the positions of two Weyl points labelled by diamonds.
After removing the high-band thermal effects~\cite{Sun_2018A}, we obtain the corrected locations of Weyl points at $\boldsymbol{q^{\text{W}}}=(0,0,\pm(0.54\pm0.02)E_{\rm{r}})$, agree well with numerical calculations (For details, see Methods).

By varying the detuning term in Eq.(\ref{3DSOC}), we reconstruct the 3D Weyl semimetal bands for $\delta_0=0$ with the same protocol.
To precisely characterize the topology, the high-band effects are taking into account for correction.
The 3D spin structure with for original data (left) and high-band correction (right) are presented in Fig.\ref{fig3}(a).
Three typical 2D spin textures are shown in the insets, with band inversion rings marked.
For $q_z=0.7\pi$, the ring encloses the $\Gamma$ point, indicating $\mathcal{C}=1$. For $q_z=0.3\pi$, the ring encloses the $\text{M}$ point, indicating $\mathcal{C}=-1$. For $q_z=0.5\pi$, the ring touches $\text{X}_1$ and $\text{X}_2$ points, indicating $\mathcal{C}$ jumping by 2 across the plane.
The corresponding four Weyl points locate at $\boldsymbol{q^{\text{W}}}=(0,\pi,\pm(0.52\pm0.03)\pi)$ and $(\pi,0,\pm(0.52\pm0.03)\pi)$, which is also confirmed by the $P(\Lambda_j,q_z)$ and $\Theta(q_z)$ calculation, in Fig.\ref{fig3}(b).

\section{Measuring the Weyl nodes with quench dynamics}

The cold atom system enables the quench studies to characterize the topology with high controllability. We then elucidate the Weyl nodes by quench dynamics. Atom clouds with temperature of 200nK are initially prepared in $\mid\uparrow\rangle$. Raman couplings are effectively excluded with initial detuning $\delta_i=-200E_{r}$. Thereupon, $\delta$ is suddenly switched to $-0.5E_{r}$ in 1$\mu\text{s}$.
Hence, the Hamiltonian is quenched from topological trivial region to IWSM region.
CCD-Y is applied to record the time evolution of spin after the quench, with spin-resolved ToF imaging.

Fig.\ref{fig4}(a) presents the dynamic evolution of spin polarizations $P(q_x,q_z,t)$ in the Brillouin zone, from 0 to 2ms (with $q_y$ being integrated out).
We further get $P(q_z,t)$ by integrating out $q_x$ to emphasise the dynamics along $q_z$ direction, as plotted aside.
For typical $q_z$, $P(q_z,t)$ is plotted as function of evolution time $t$ in Fig.\ref{fig4}(b).
When $q_z$ closed to $\pm\pi$, the evolution of $P(q_z,t)$ is dominated by fast $s$-$p$ band oscillations.
When $q_z$ closed to $0$, for example $q_z=0$ and $\pm0.3\pi$, $P(q_z,t)$ exhibits both fast and slow oscillations.
We fit those evolutions with a double-frequency damped oscillator (For details, see Methods), and extract the low frequency component $f_{\text{low}}$ as functions of $q_z$, as plotted in Fig.\ref{fig4}(c). Around $q_z=-0.64\pi$ and $q_z=+0.62\pi$, $f_{\text{low}}$ reaches kink minima which we explain below.

The critical points $q_z=-0.64\pi$ and $0.62\pi$ are identified as positions of the two Weyl nodes, close to the result obtained from equilibrium approach.
In topological region with $|q_z|<|q_z^w|$, the 2D bands of $\mid\uparrow\rangle$ and $\mid\downarrow\rangle$ are inverted, with band inversion ring in the 2D planes~\cite{LZhang2018}.
The dynamics is dominated by on-resonance Raman-Rabi oscillations between the two inverted bands~\cite{Sun-2018B, Yi}, corresponding to the low frequency oscillations with large amplitude.
In trivial region, the bands of $\mid\uparrow\rangle$ and $\mid\downarrow\rangle$ are separated, therefore the on-resonance Raman-Rabi oscillations vanish.
At the Weyl points, the two bands touch together at a single point, make $f_{\text{low}}$ reaches minima~\cite{Theory_Liu}, while not zero due to background noises.
The numerical results of low frequency components using experimental parameters are also presented in Fig.\ref{fig4}(c) (see Methods for details), well supporting our observation.

\section{Conclusion and Outlook}

We have realized the 3D SO coupling for ultracold atoms and the ideal Weyl semimetal band (IWSM) which consists of only two Weyl cones and is hard to engineer in solid materials. The IWSM band is experimentally identified by resolving the Weyl points, which are detected by virtual slicing imagining technique and further measured in quench dynamics. We also demonstrate the high tunability by engineering the semimetal band with more Weyl points for comparison. The experimental observations are consistent with the numerical results.

The realization of the IWSM opens a broad avenue to explore exotic quantum phenomena with this optimal Weyl semimetal platform. For example, for the chiral anomaly~\cite{Nielsen1983,Vishwanath2014} the indirect signature was studied in condensed matter physics by measuring negative magnetic resistance~\cite{Huang2015}. The clean and direct probe of this phenomenon may be achieved in the IWSM band with high controllability of the ultracold atoms. Another related but more exotic phenomenon is the chiral magnetic effect~\cite{Fukushima2008}, which states that when an energy shift is introduced between the Weyl points by breaking reflection symmetry, a nonzero chiral current may be generated by a magnetic field, even without applying external electric field~\cite{Franz2013}. This effect is highly debating in condensed matter physics, but may be resolved unambiguously based on current IWSM bands and precise control of atom distributions.

The present protocol for 3D SO coupling and IWSM band is generic, and can be immediately applied to fermion systems, in which case, the various correlated phases shall be accessed by tuning strong interactions. In particular the highly-sought-after topological superfluids~\cite{ZhangC2008,Sato2009} could be achieved for the 3D SO coupled Fermi gases, where the mean-field theory captures essential physics~\cite{Chan2017}, with higher reliability than similar attempts for the 1D or 2D SO coupled systems.

\begin{addendum}
    \item[Online content] Any methods, extended data, supplementary information, acknowledgements, details of author contributions and competing interests are available online.
\end{addendum}

\bibliographystyle{naturemag}