The aim of this article devotes to establishing the (L^1)-decay of cubic order spatial derivatives of solutions to the Navier–Stokes equations, which is a long-time challenging problem. To solve this problem, new tools have to be found to overcome these main difficulties: (L^1-L^1) estimate fails for the Stokes flow; the projection operator (P:,L^1(mathbb {R}^n_+)rightarrow L^1_sigma (mathbb {R}^n_+)) becomes unbounded; the steady Stokes’s estimates does not work any more in (L^1(mathbb {R}^n_+)). We first give the asymptotic behavior with weights of negative exponent for the Stokes flow and Navier–Stokes equations in (L^1(mathbb {R}^n_+)), and these are also independent of interest by themselves. Secondly, we decompose the convection term into two parts, and translate the unboundedness of projection operator into studying an (L^1)-estimate for an elliptic problem with homogeneous Neumann boundary conditions, which is established by using the weighted estimates of the Gaussian kernel’s convolution. Finally, a crucial new formula is given for the fundamental solution of the Laplace operator, which is employed for overcoming the strong singularity in studying the cubic order spatial derivatives in (L^1(mathbb {R}^n_+)).