The statement that the second term in your momentum flux in terms of the relative velocity corresponds to the viscous stresses is not correct. Instead it corresponds to the entire stress tensor and therefore includes the pressure term. This can be shown taking momenta of the Boltzmann equation assuming only a small perturbation from macroscopic equilibrium, given by a Maxwell-Boltzmann equilibrium distribution, with the use of perturbation theory, the so called Chapman-Enskog analysis. The entire momentum tensor in this form emerges then macroscopically from microscopic interactions of particles. This is a mathematically very involved procedure and therefore I will only sketch out the derivation below.
I will use my own nomenclature as I always try to keep it consistent throughout my answers so I can refer to them in future answers. My momentum flux tensor $\Pi_{ij}$ is the equivalent to your $J_{ij}$, the viscous stresses $\tau_{ij}$ replace your $\sigma_{ij}$ and I will use $\sigma_{ij}$ for the entire stress tensor consisting of pressure and viscous stresses. For the individual particle's velocity I will use $\vec \xi$ instead of $\vec v$, for the main flow $\vec u$ and thus for the relative velocity $\vec v := \vec \xi - \vec u$ instead of $\vec w$.
Momentum flux from a continuum mechanics approach
With this nomenclature the momentum flux tensor is given as
$$\Pi_{ij} := \rho u_i u_j + \underbrace{p \delta_{ij} - \tau_{ij}}_{- \sigma_{ij}}$$
where the viscous stresses $\tau_{ij}$ for a Newtonian fluid can be derived from a macroscopic continuum perspective
$$\sigma_{ij} = - p \delta_{ij} + \underbrace{2 \mu S_{ij} - \frac{2}{3} \mu \sum\limits_{k \in \mathcal{D}} S_{kk} \delta_{ij}}_{\tau_{ij}},$$
and the strain rate tensor
$$S_{ij} := \frac{1}{2} \left( \vec \nabla \otimes \vec u + (\vec \nabla \otimes \vec u)^T \right) = \frac{1}{2} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) $$
is the symmetric part of the velocity gradient.
It allows us to rewrite the momentum equation in fluid mechanics
$$\frac{\partial (\rho u_i )}{\partial t} + \sum\limits_{j \in \mathcal{D}}\frac{\partial (\rho u_i u_j )}{\partial x_j} = \sum\limits_{j \in \mathcal{D}} \frac{\partial \sigma_{ij}}{\partial x_j } + \rho g_i$$
that itself can be derived by a simple force balance on a continuum element in a pretty convenient way
$$\frac{\partial (\rho u_i )}{\partial t} + \sum\limits_{j \in \mathcal{D}}\frac{\partial \Pi_{ij}}{\partial x_j} = \rho g_i.$$
Kinetic theory of gases and the Boltzmann equation
In kinetic theory one tries to describe a dilute fluid as a collection of particles that generally only interact with collisions but in more sophisticated methods also in far-field interactions. There exists a range of models such as the simple $1/8$ model generally used to link pressure to changes of momentum and further to the kinetic energy of a gas that does not consider collisions but instead uses symmetries and equipartition to account for them.
A more sophisticated approach for this - based on multi-body mechanics and the Louiville equation - is finding an evolution equation that describes how the particle distribution function
$$ f := \frac{d N}{d \vec x \, d \vec \xi}$$
(where $N$ is the number of particles) evolves over time due to linear motion of and collisions between particles. This particle distribution can be seen as an extended concept of density that might still hold if the limiting values, that one assumes to exist in fluid mechanics
$$ \rho := \lim\limits_{\Delta V \rightarrow 0} \frac{\Delta m}{\Delta V} \qquad \sigma_{ij} := \lim\limits_{\Delta A_i \rightarrow 0} \frac{\Delta F_j}{\Delta A_i} \qquad g_i := \lim\limits_{\Delta m \rightarrow 0} \frac{\Delta G_i}{\Delta m}, $$
can't be found.
The macroscopic values can be found by summing up through integration, for example the density and the momentum can be found to
$$\rho = m_P \int f d \vec \xi, \qquad \rho \vec u_i = m_P \int \xi_i f d \vec \xi.$$
The evolution equation of this distribution is given by the Boltzmann equation
$$ \underbrace{ \frac{\partial f}{\partial t} + \vec \xi \boldsymbol{\cdot} \vec \nabla f + \frac{\vec F}{m} \boldsymbol{\cdot} \vec \nabla_{\vec \xi} f }_\text{Propagation} = \underbrace{ \Omega(f) }_\text{Collision} $$
where the left-hand side can be derived by basic chain rule (the third term is a force term resulting from external forces) while the collision term on the right-hand side can be derived analytically for elastic collisions (I have posted a possible derivation some time ago over here as well as the corresponding collision cross-section $A_c$), the so called "Stoßzahlansatz" to
$$ \Omega_{Stoß} = df |_{\text{collision}} = \frac{\Delta N_{\text{gain}} - \Delta N_{\text{loss}}}{\Delta V \Delta \vec \xi \Delta t} = \int\limits_{ \vec \xi_1 } \int\limits_{ A_c } |\vec g| (f_1' f' - f f_1 ) d A_c d \vec \xi_1. $$
Furthermore one can find a equilibrium distribution, the so called Maxwell-Boltzmann distribution (derived it already in another post),
$$ f^{(eq)}(\vec x, |\vec v|, t) = n \underbrace{ \frac{1}{(2 \pi R_m T )^{\frac{3}{2}}} \, e^{-\frac{|\vec v|^2}{2 R_m T}} }_\text{Gaussian distribution}. $$
(where $n$ is the particle density $n := \frac{N}{V} = \frac{\rho}{m_P}$) and proof that a system without external disturbances evolves towards such a distribution.
BGK collision operator
Now that we know that this system in the bulk evolves towards a Maxwell-Boltzmann distribution this might give us the idea to not use the exact collision term but instead only model every collision as a tiny step to the final equilibrium by a tiny relaxation time $\tau$, a model named after the three scientists that introduced it, the Bhatnagar–Gross–Krook (BGK) collision operator:
$$ \Omega_{BGK} := \frac{f^{(eq)} - f}{\tau}$$
The reason that I mention this is that you may find the following Chapman-Enskog expansion for both, the Boltzmann "Stoßzahlansatz" as well as the simplified BGK operator in the literature.
Perturbation analysis: Chapman-Enskog expansion
You might be telling yourself now: "Well, that's a fancy-ass equation in some obscure variable but in the end it tells me nothing about the macroscopic nature of the fluid." That is luckily not the case, you can gain some insight into the macroscopic behaviour of the fluid by considering moments of the particle distribution as these are linked to the macroscopic behaviour of the dilute gas, this is done by perturbation theory.
Perturbation theory is a mathematical method for finding approximate solutions to a variety of problems including transcendental and differential equations and traditionally is used in multi-body and celestial mechanics. Starting with the exact solution of a simple unperturbed problem we expand it using a power series in a small parameter $\epsilon$ that quantifies the deviation from the exactly solvable problem $x_0$.
$$ x = \sum_{n=0}^{\infty} \epsilon^n x_n = x_0 + \epsilon x_1 +\epsilon^2 x_2 + \cdots$$
Normally an approximate perturbation solution is obtained by truncating the series and usually keeping the low-order terms only. Now this can be applied to approximate a system where the exact solution to the unperturbed problem $x_0$ (in our case the Maxwell-Boltzmann distribution) is known using the fundamental theorem of perturbation theory:
If
$$ x_0 + \epsilon x_1 +\epsilon^2 x_2 + \cdots + \epsilon^n x_n + \mathcal{O}(\epsilon^{n+1})=0 $$
for $\epsilon \rightarrow 0$ and $x_0, x_1, \cdots$ independent of $\epsilon$, then
$$ x_0 = x_1 = x_2 = \cdots = x_n = 0 $$
For systems where this regular perturbation approach leads to secular terms, terms that can't be cancelled by choosing parameters accordingly and hence the solution grows without bound, there exist more sophisticated approaches that introduce additional scaling for variables. The Chapman-Enskog analysis is such a multiple-scales perturbation series.
In the Chapman-Enskog analysis this perturbation series takes the form
$$f(t_0, t_1, t_2, \ldots) = \underbrace{f^{(0)}(t_0)}_{f^{(eq)}} + \epsilon f^{(1)}(t_1) + \epsilon^2 f^{(2)}(t_2) + \mathcal{O}(\epsilon^3)$$
where
$$ t_n = \epsilon^n t_0$$
and $\epsilon$ is generally seen as the Knudsen number
$$ Kn := \frac{\lambda}{L} \phantom{spacespace} \frac{\text{mean free path}}{\text{representative physical length scale}}$$
as it characterises the deviation from a collision-dominated (right-hand side) macroscopic view in the dimensionless Boltzmann equation (here with BGK operator)
$$f^*=f \frac{c_0^3 L^6}{n}, \qquad t^*=t \frac{c_0}{L}, \qquad x_i^*=\frac{x_i}{L}, \qquad \xi_i^*=\frac{\xi_i}{c_0}, \qquad g^*=g \frac{L}{c_0^2}, \qquad \tau^*=\tau \frac{c_0}{\lambda},$$
$$\frac{\partial f^*}{\partial t^*} + \xi_j^* \frac{\partial f^*}{\partial x_j^*} + g_j^* \frac{\partial f^*}{\partial \xi_j^*} = \frac{1}{Kn} \frac{1}{\tau^*} \left( f^{(eq)*} - f^* \right).$$
Only the first three contributions are required to derive the equation system for a dense fluid, the Navier-Stokes equations while the meaning of higher-order contributions, such as the Burnett- and Super-Burnett equations, is not completely known. We apply this perturbation series to the Boltzmann equation and then evaluate the moments of it to obtain differential equations for mass, momentum and energy.
Applying this perturbation series only for equilibrium, for the first term $n=0$, we can find the Euler equation for inviscid flow (meaning the pressure emerges from $f^{(eq)}$ and does not vanish even in thermodynamic equilibrium) while considering the first three contributions (up to order $n=2$) we can find the Navier-Stokes equation with the following relationships between microscopic distributions and macroscopic variables:
$$ m_P \int (\xi_i - u_i) (\xi_j - u_j) f d \vec \xi = m_P \int v_i v_j f d \vec \xi = p \delta_{ij} - \tau_{ij} = - \sigma_{ij}, $$
$$ \frac{m_P}{3} \int |\vec \xi - \vec u|^2 f d \vec \xi = p, $$
$$ m_P \int \xi_i \xi_j f d \vec \xi = \rho u_i u_j + \underbrace{p \delta_{ij} - \tau_{ij}}_{-\sigma_{ij}} = \Pi_{ij}.$$
In the expression for the thermodynamic pressure we can see the Stokes' hypothesis, where it is assumed that the mechanical pressure
$$\overline{p} := - \frac{1}{dim(\mathcal{D})} \sum\limits_{j \in \mathcal{D}} \sigma_{jj} = - \frac{1}{3} \left( \sigma_{11} + \sigma_{22} + \sigma_{33} \right)$$
corresponds to the thermodynamical pressure
$$ p \approx \overline{p}.$$
The normal viscous stresses are measured relative to an average value, the pressure, and thus the pressure corresponds to the average fluctuations of the relative velocity in the normal directions.
It is also interesting to note that the pressure emerges from the equilibrium terms $f^{(0)} = f^{(eq)}$ while the dissipative parts, such as the viscous stresses, emerge from the higher-order non-equilibrium contributions ($f^{(1)}$ and $f^{(2)}$) in the Chapman-Enskog expansion. From that we can clearly see that even in equilibrium there are fluctuations in velocity but they conform to the Maxwell-Boltzmann distribution and are thus isotropic in space. In non-equilibrium these fluctuations are in dis-balance and result in viscous stresses.
