\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}
