\section{Introduction and summary}
According to their modern description and construction in chiral effective field
theory, nuclear forces are organized in a hierarchical way \cite{evgeni,hammer}. For 
generic few- and many-body observables the contributions arising from two-nucleon 
interactions are larger than those from three-nucleon forces, and the latter are 
again more important than possible  corrections due to four-body forces. After a
renormalization-group evolution to lower resolution scale such chiral low-momentum 
interactions exhibit desirable convergence properties in 
perturbative calculations of many-nucleon systems and infinite nuclear matter 
\cite{vlowkreview,achimnuc,3bodycalc,hebeler}. In particular, the inclusion of the 
leading-order chiral three-nucleon interaction (consisting of contact, 
$1\pi$-exchange, and $2\pi$-exchange components) is essential in order to achieve 
saturation of isospin-symmetric nuclear matter in this framework. The calculation of 
the sub-leading chiral three-nucleon interaction, built up by many pion-loop diagrams 
etc., has been completed in ref.\,\cite{3bodyn3lo}. By employing chiral two-, three-, 
and four-nucleon interactions up to order N$^3$LO the equation 
of state of neutron matter at zero temperature has been studied in detail by the 
Darmstadt group \cite{achim,darmstadt}. A good convergence pattern has been observed 
for the second and third order particle-particle diagrams with two-nucleon 
potentials only and when including additionally three-nucleon forces in the form of 
density-dependent two-body interactions. Moreover, in ref.\,\cite{darmstadt} the 
N$^3$LO three-body and four-body contributions to the energy per particle of 
isospin-symmetric nuclear matter have been estimated and large effects from the 
sub-leading chiral three-nucleon force have been found together with small 
corrections from the leading-order chiral four-nucleon interactions.

An analytical calculation of chiral four-body interactions in nuclear matter and pure
neutron matter has been performed recently in ref.\,\cite{4body}. The long-range 
four-nucleon interaction is generated by pion-exchanges that are combined by the 
chiral $4\pi$-vertex (off-shell $\pi\pi$-interaction) or the chiral 
$N\!N3\pi$-vertex. The contributions arising from the 2-ring and 1-ring diagrams 
related to this ''irreducible'' four-nucleon interaction turned out to be very small, 
with values below $0.6\,$MeV for densities $\rho<0.4\,$fm$^{-3}$. However, when 
including the $\Delta(1232)$-isobar as an explicit degree of freedom and counting 
the delta-nucleon mass splitting $\Delta = 293$\,MeV as a small scale (comparable to 
the pion mass $m_\pi$ or the Fermi momentum $k_f$) a new class of leading-order 
long-range four-nucleon interactions arises. These are mediated by two-fold  
$\Delta(1232)$-isobar excitation of nucleons and the exchange of three pions. The 
analytical calculation of the pertinent 3-ring, 2-ring, and 1-ring diagrams at 
four-loop order has lead to a moderately repulsive contribution of $\bar E(\rho_0)
=2.35\,$MeV at nuclear matter saturation density $\rho_0 = 0.16\,$fm$^{-3}$. However, 
the curve for $\bar E(\rho)$ rises strongly with the density, reaching about 10 
times that value at $\rho= 2\rho_0$. In the case of pure neutron matter the same 
$\Delta(1232)$-induced four-neutron interactions lead to a repulsive contribution 
that is about half as strong (see Fig.\,13 in ref.\,\cite{4body}). 

The complete set of leading-order four-nucleon interactions in chiral effective 
field theory has been constructed by Epelbaum in refs.\,\cite{evgeni4na,evgeni4nb}
using the method of unitary transformations. This method allows to project the 
dynamics of the interacting pion-nucleon system into the purely nucleonic subspace 
relevant for few-nucleon systems at energies below the pion-production threshold. 
In addition to the ''irreducible'' chiral four-nucleon interaction mentioned above, 
which follows (via the Feynman diagram technique) directly from the  vertices of the 
chiral Lagrangian ${\cal L}_{\pi N}^{(1)}+ {\cal L}_{\pi\pi}^{(2)}$, the method of unitary 
transformations gives rise to eight classes of ''reducible'' four-nucleon diagrams. 
In ref.\,\cite{evgeni4nb} it is shown that disconnected diagrams always lead to 
vanishing $4N$-forces and with the kinematically suppressed pion-energies at a 
Weinberg-Tomozawa  $\pi\pi N\!N$-vertex one obtains that actually three classes (III, 
VI, and VIII) vanish at leading order. The nonvanishing $4N$-interactions, grouped 
into the classes I, II, IV, V, and VII, involve multiple pion-exchanges and a 
spin-spin contact-interaction with coupling constant $C_T$. The associated isoscalar 
central $N\!N$-contact term proportional $C_S$ is absent because it drops out of the
commutators arising in the construction of the reducible $4N$-interaction. Note 
that the ``irreducible'' four-nucleon interaction, which constitutes a part of class 
II, is called $V^e+V^f$ in ref.\,\cite{evgeni4na}. Moreover, one finds that in an 
anti-symmetrized four-neutron state only the (diagonal) matrix-element of $V^a$ 
(class I) is nonvanishing and an explicit expression has been given for it in 
eq.(A14) of ref.\cite{darmstadt}. Analogous expressions for the four-nucleon 
matrix-elements of the other classes have so far not been published.

The purpose present paper is to evaluate analytically the five classes (VII, V, IV, 
II, I) of reducible chiral four-nucleon interactions ($V^n, V^l, V^k, V^c, V^a$ in 
the notation of ref.\,\cite{evgeni4na}) in isospin-symmetric nuclear matter and 
pure neutron matter. In order to facilitate a uniform treatment we will first
introduce in section\,2 a general representation of $4N$-interactions in terms of a 
product of four (effective) single-nucleon vertices and three ''propagator''-factors.
The latter will be ordinary pion-propagators, their squares, or just constants. The 
first-order ''Hartree-Fock'' contribution to the energy density of nuclear matter is 
obtained from the connected four-body diagram by closing its four nucleon-lines in 
all possible ways. The  actually occurring chiral vertices have the convenient 
property that closing any nucleon-line to itself leads to a vanishing spin-trace or 
isospin-trace. Therefore, one needs to consider here only the closed 2-ring and 
1-ring diagrams generated by the reducible chiral $4N$-interactions. These four-loop 
diagrams are then evaluated individually for the five different classes (in the order 
VII, V, IV, II, I) in the subsections 2.1 up to 2.5. In most cases the occurring 
integral over the product of four Fermi spheres of radius $k_f$ can be reduced to an 
easily manageable one- or two-parameter integral, that depends on the dimensionless 
ratio $k_f/m_\pi$. The 1-ring diagrams with crossed pion-exchanges related to class 
II and I allow only for a partial analytical reduction, such that six-parameter 
integrals remain. The numerical results for the energy per particle $\bar E(\rho)$ 
as a function of the nucleon density $\rho=2k_f^3/3\pi^2$ are presented and discussed 
in section 3. We consider two opposite values of the spin-spin contact-coupling, 
$C_T=0.22\,$fm$^2$ and  $C_T=-0.45\,$fm$^2$, taken from table I in ref.\,\cite{
darmstadt}. One observes substantial cancelations between the contributions from
 different classes. In particular, the $C_T$-independent classes II and I provide 
the largest repulsive and attractive contribution, respectively. Altogether, the net 
attraction generated by the reducible chiral four-nucleon interaction stays above 
values of $-1.3\,$MeV for densities $\rho<2\rho_0=0.32$\,fm$^{-3}$. This feature 
holds for the whole range $|C_T|<0.5\,$fm$^2$ of the spin-spin contact-coupling 
$C_T$. In the case of pure neutron matter, where only class I contributes, the 
result for the energy per particle $\bar E_n(\rho_n)$ is determined entirely by a 
2-ring diagram, since the three types of 1-ring diagrams compensate each other almost 
completely. The corresponding attractive energy per particle goes approximately as 
$\rho_n^{7/3}$ and it reaches a value of about $-2.1$\,MeV at $\rho_n=0.4$\,fm$^{-3}$.  

In summary one finds that the chiral four-nucleon correlations studied in this work 
are at least one order of magnitude smaller than those provided by the strongly 
coupled $\pi N\Delta$-system \cite{4body} with its small mass-gap of $293$\,MeV.  
 



\section{Diagrammatic calculation of four-body interactions in nuclear 
matter}
In this section we derive analytical expressions for the first-order contributions 
to the energy per particle of nuclear matter as they arise from the reducible chiral 
four-nucleon interactions of ref.\,\cite{evgeni4nb}. In order to facilitate a 
uniform treatment of the five different classes (VII, V, IV, II, I) we introduce 
first a generic representation of $4N$-interactions in terms of a product of four 
(effective) single-nucleon vertices and three ''propagator''-factors. The latter 
will be ordinary pion-propagators, their squares, or just constants. Fig.\,1 shows 
the intended generic form of the reducible chiral four-nucleon interaction, 
represented as an on-shell scattering process $N_1N_2N_3N_4\to N'_1N'_2N'_3N'_4$. The 
labels $a,b,c$  denote isospin-indices of exchanged pions and $i,j$ denote 
components of (spin)-vectors. These vector-indices are introduced by the spin-spin 
contact-term proportional to $C_T$, and if present, the left and/or right 
''propagator''-factor is a constant. The momentum transfers $\vec q_a=\vec p_1\,\!\!'-
\vec p_1$, $\vec q_c=\vec p_4\,\!\!'-\vec p_4$ and $\vec q_b = \vec p_2\,\!\!'-
\vec p_2+\vec q_a= \vec p_3 -\vec p_3\,\!\!'-\vec q_c$ are given by differences of 
out- and in-going nucleon momenta.  
\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.6,clip]{4-bodydiagram.eps}
\end{center}
\vspace{-11.3cm}
\caption{Generic form of the ''reducible'' chiral 4N-interaction. The square-box and 
filled circle symbolize effective single-nucleon vertices. The labels $a,b,c$ denote 
isospin-indices of exchanged pions and $i,j$ denote components of (spin)-vectors.}
\end{figure}

From the generic four-body diagram in Fig.\,1 the first-order "Hartree-Fock" 
contribution to the energy density $\rho \bar E(\rho)$  of nuclear matter is obtained 
by closing the four nucleon-lines in all possible ways. Inspecting the factorized 
expressions for the five reducible chiral $4N$-interactions written in 
eqs.(1,5,7,10,19) below, one sees that closing any nucleon-line to itself leads to a 
vanishing spin-trace or isospin-trace. In the case of four neutrons one has either 
a vanishing spin-trace or a zero from the totally anti-symmetric $\epsilon^{abc}
$-tensor, since $\pi^0$-exchange fixed $a=b=3$. As a consequence of this property, 
closed 4-ring and 3-ring diagrams vanish trivially and one needs to consider here 
only the closed 2-ring and 1-ring diagrams generated by the reducible chiral 
$4N$-interactions.
\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.6,clip]{2-Ring.eps}
\end{center}
\vspace{-3.cm}
\caption{Closed two-ring diagrams generated by the reducible chiral 4N-interaction. 
A symmetry factor $1/2$ occurs, if the  square-box and filled-circle vertex are 
identical. The diagram at the bottom vanishes due to the zero-momentum 
$\vec q_b = \vec 0$ of the exchanged central pion.}
\end{figure}
\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.6,clip]{1-Ringuzx.eps}
\end{center}
\vspace{-2.cm}
\caption{Closed one-ring diagrams generated by the reducible chiral 4N-interaction. 
The diagrams in the first, second, and third row are referred to as U-type, Z-type, 
and X-type, respectively. A symmetry factor $1/2$ belongs to the Z-type diagrams, 
if the  square-box and filled-circle represent identical vertices.}
\end{figure}


The three possible 2-ring topologies are shown in Fig.\,2. One notices that the 
diagram at the bottom vanishes due to the zero-momentum $\vec q_b = \vec 0$ of the 
exchanged central pion. At least one of the adjoined vertices involves $\vec q_b$ 
linearly. Furthermore, the six possible 1-ring topologies are shown in Fig.\,3. 
According to the shape of the dashed (pion) line the diagrams in the first, second, 
and third row are referred to as U-type, Z-type, and X-type, respectively. Note 
also that if the square-box and filled-circle vertex represent identical vertices 
(this happens for class VII and I) a symmetry factor $1/2$ belongs to the 2-ring 
diagrams and to the Z-type 1-ring diagrams, and obviously identical copies of the 
U-type and X-type 1-ring diagrams need to be discarded. After these preparations we 
can now turn to the evaluation of the individual 1-ring and 2-ring diagrams 
generated by the five classes of reducible chiral $4N$-interactions.    
\subsection{Class VII}
We start (in reverse order) with class VII, which is quadratic in the spin-spin 
contact-coupling $C_T$. The corresponding factorized expression for this chiral
$4N$-interaction reads:  
\begin{equation}  V^{n} = -{C_T^2 g_A^2 \over f_\pi^2} \, \sigma_1^i \,
(\vec \sigma_2 \times \vec q_b)^i \, \tau_2^b \,{1\over(m_\pi^2+\vec q_b^{\,2})^2} 
\,(\vec \sigma_3\times \vec q_b)^j \, \tau_3^b \,\sigma_4^j\,.\end{equation} 
Note that there is an extra factor $2$ in comparison to eq.(3.60) in 
ref.\,\cite{evgeni4nb} due to the interchange symmetry $(12)\leftrightarrow(43)$ 
between pairs of nucleons. A straightforward  evaluation of the left 2-ring 
diagram in Fig.\,2 gives the following contribution from class VII to the energy per 
particle of isospin-symmetric nuclear matter:
\begin{eqnarray} \bar E(\rho) &=& {3C_T^2 g_A^2m_\pi^7 \over (2\pi^2 u)^3f_\pi^2}
\int_0^u\!\!dx\int_0^u\!\!dy \, x(u^2-y^2)(u-x)^2(2u+x) \bigg\{-8x y +6x \nonumber \\ 
&& \times \big[\arctan(2x+2y)- \arctan(2x-2y) \big] +(2x^2 -2y^2 -1)\ln{1+4(x+y)^2 
\over 1+4(x-y)^2}\bigg\}\,, \end{eqnarray}
with the dimensionless variable $u =k_f/m_\pi$. The nucleon density $\rho= 2k_f^3/3
\pi^2$ is related to the Fermi momentum $k_f$ in the usual way. In order to arrive 
at a double-integral representation we have used the reduction formula:
\begin{equation} \int_0^u\!\!dp\! \int_{-1}^1\!\!dz \, p^2 F\big(p z+\sqrt{u^2-p^2
(1-z^2)}\,\big) = \int_0^u\!\!dy \, (u^2-y^2) F(2y)\,, \end{equation} 
to eliminate a complicated angular integral, and a master formula for integrating $
F(|\vec p_1\!-\!\vec p_2|)$ over two Fermi spheres $|\vec p_{1,2}|<k_f$. The latter 
introduces the weighting-function $x^2(u-x)^2(2u+x)$. Actually, there exists an 
analytical solution of the double-integral in eq.(2), namely: 
\begin{eqnarray} \bar E(\rho) &=& {C_T^2g_A^2m_\pi^7 \over 140\pi^6f_\pi^2}
\Bigg\{ {u\over 160}-{31u^3\over 48}+{229u^5\over 20}-{179u^7\over 10} 
\nonumber \\ && +\bigg({3\over 16}+{343u^2\over 40}+20u^4-{142u^6\over 3}
\bigg)\arctan 2u\nonumber \\ && +\bigg({176u^6\over 3}-10u^4-{343u^2\over 
80}-{3 \over 32}\bigg) \arctan 4u \nonumber \\ && +\bigg(28u^5-{32u^7 
\over 5}+{35u^3 \over 8}-{7u\over 8}-{5 \over 128u} -{1 \over 3840 u^3}
\bigg) \ln(1+4u^2) \nonumber \\ && + \bigg({32u^7 \over 5}-28u^5-{35u^3 
\over 8}+{7u \over 32}+{5 \over 512u} +{1 \over 15360u^3} \bigg) 
\ln(1+16u^2) \Bigg\}\,. \end{eqnarray} 
Note that in the chiral limit $m_\pi \to 0$ only a term $u^7(128\ln 2-179)/10= 
-9.0277\,u^7$ remains. Let us also comment on the features of the other diagrams 
for class VII: The right 2-ring diagram in Fig\,2 leads to a vanishing isospin-trace 
in symmetric nuclear matter, whereas in pure neutron matter the left and right 
2-ring diagram cancel each other. The U-type and Z-type 1-ring diagrams in Fig.\,3 
cancel each other after taking the spin-trace, and the X-type 1-ring diagram leads 
to a vanishing spin-trace by itself.
\subsection{Class V}
Next we come to class V, which is induced by the Weinberg-Tomozawa $\pi\pi 
N\!N$-vertex in combination with the spin-spin contact-coupling $C_T$. The 
corresponding factorized expression for this chiral $4N$-interaction reads 
\cite{evgeni4nb}:  
\begin{equation}  V^l={C_T g_A^2\over 8f_\pi^4}\,\vec\sigma_1\cdot \vec q_a \,
\tau_1^a\,  {1\over m_\pi^2+\vec q_a^{\,2}} \, \epsilon^{abc} \, \tau_2^c \,
{1\over m_\pi^2+\vec q_b^{\,2}} \,(\vec \sigma_3\times \vec q_b)^j \, \tau_3^b
\,\sigma_4^j\,,\end{equation} 
and it obviously vanishes for four neutrons, where only $\pi^0$-exchange is 
possible ($\epsilon^{33c}=0$). The evaluation of the X-type 1-ring 
diagram in Fig.\,3 gives the following contribution from class V to the energy per 
particle of isospin-symmetric nuclear matter:
\begin{eqnarray} \bar E(\rho) &=& {3C_T g_A^2m_\pi^7 \over (2\pi)^6f_\pi^4u^3}
\int_0^u\!\!dx\int_0^u\!\!dy \, {x(u^2-y^2)\over 1+4x^2}(u-x)^2(2u+x) \bigg\{
-4x y \nonumber \\ && \times (1+4x^2+4y^2)+{1\over 4}\big[1+4(x+y)^2\big]
\big[1+4(x-y)^2\big] \ln{1+4(x+y)^2 \over 1+4(x-y)^2}\bigg\}\,. \end{eqnarray}
Note that the integral $\int_0^u\!dy\dots$ could be solved in terms of elementary 
functions ($\ln, \arctan$), but we refrain from presenting this rather lengthy 
expression. The same applies to the subsequent results given in eqs.(8,9,11,20,21).
Concerning the other diagrams for class V, one finds that both 2-ring diagrams in 
Fig.\,2 lead to a vanishing spin- and isospin-trace. Moreover, the U-type and Z-type 
1-ring diagrams in Fig.\,3  cancel each other after taking the spin-trace.

\subsection{Class IV}
Next we come to class VI, which is induced by the spin-spin contact-coupling $C_T$ in 
combination with a ''contraction'' of two ordinary $\pi N$-vertices. The 
corresponding factorized expression for this chiral $4N$-interaction reads 
\cite{evgeni4nb}: 
\begin{equation}  V^k = -{C_T g_A^4 \over 4f_\pi^4} \, \vec\sigma_1\cdot \vec q_a \,
\tau_1^a\,  {1\over m_\pi^2+\vec q_a^{\,2}} \, \big[\epsilon^{abc} \, \tau_2^c \,
\vec q_a\cdot\vec q_b+\delta^{ab}\,\vec\sigma_2\cdot(\vec q_a\times\vec q_b)\big]\,
{1\over (m_\pi^2+\vec q_b^{\,2})^2} \,(\vec \sigma_3\times \vec q_b)^j \, \tau_3^b
\,\sigma_4^j\,.\end{equation} 
Whereas the left 2-ring diagram has a vanishing isospin-trace, one obtains from the 
right 2-ring diagram the following contribution to the energy per particle 
of isospin-symmetric nuclear matter:
\begin{eqnarray} \bar E(\rho) &=& {3C_T g_A^4m_\pi^7 \over 16\pi^6f_\pi^4u^3}
\int_0^u\!\!dx\int_0^u\!\!dy \, {x(u^2-y^2)\over 1+4x^2}(u-x)^2(2u+x) \nonumber 
\\ && \times \bigg\{x y(1-12x^2+4y^2) +8x^3\big[\arctan(2x+2y)- \arctan(2x-2y)
\big] \nonumber \\ && +\Big(3x^4-{3x^2 +y^2 \over 2}-2x^2y^2-y^4-{1 \over 
16} \Big) \ln{1+4(x+y)^2 \over 1+4(x-y)^2} \bigg\} \,. \end{eqnarray}
By rewriting the spin- and isospin-factors as a product of two anti-commutators one 
finds that the U-type and Z-type 1-ring diagrams add to zero. On the other hand, the
X-type 1-ring diagrams in Fig.\,3 give a contribution to $\bar E(\rho)$ of the form:
\begin{eqnarray} \bar E(\rho) &=& {3C_T g_A^4m_\pi^7 \over 16\pi^6f_\pi^4u^3}
\int_0^u\!\!dx\int_0^u\!\!dy \, {x(u^2-y^2)\over 1+4x^2}(u-x)^2(2u+x) \nonumber 
\\ && \times \bigg\{ 2x y(1+4x^2+4y^2) +8x^3\big[\arctan(2x-2y)-\arctan(2x+2y)
\big] \nonumber \\ &&  +\Big(x^2-2x^4 +4x^2y^2-y^2-2y^4-{1\over 8}\Big)
\ln{1+4(x+y)^2 \over 1+4(x-y)^2} \bigg\}\,. \end{eqnarray}
In the case of pure neutron matter the left and right 2-ring diagram cancel each 
other, likewise do this the U-type and Z-type 1-ring diagrams,  and the X-type 1-ring 
diagrams vanish by themselves (because for four neutrons the $\epsilon^{abc}$-term 
in $V^k$ is absent).

\subsection{Class II}
Next we come to class II, which is induced by the Weinberg-Tomozawa $\pi\pi 
N\!N$-vertex in combination with the ''contraction''of two ordinary $\pi 
N$-vertices. The corresponding factorized expression for this $3\pi$-exchange 
$4N$-interaction reads 
\cite{evgeni4nb}: 
\begin{eqnarray}  V^c &=& -{g_A^4 \over 32f_\pi^6} \, \vec\sigma_1\cdot \vec q_a \,
\tau_1^a\,  {1\over m_\pi^2+\vec q_a^{\,2}} \, \epsilon^{abd}\, \tau_2^d \, {1\over 
m_\pi^2+\vec q_b^{\,2}} \, \nonumber \\ && \times \big[\epsilon^{bce} \, \tau_3^e \,
\vec q_b\cdot\vec q_c +\delta^{bc}\,\vec\sigma_3\cdot(\vec q_b\times\vec q_c)
\big]\,{1\over m_\pi^2+\vec q_c^{\,2}} \,\vec \sigma_4\cdot \vec q_c \, \tau_4^c\,,
\end{eqnarray} 
and it obviously vanishes for four neutrons, since $\epsilon^{33d}=0$. The evaluation 
of the left 2-ring diagram in Fig.\,2 gives the following contribution to the energy 
per particle of isospin-symmetric nuclear matter:
\begin{eqnarray} \bar E(\rho) &=& {3g_A^4m_\pi^7 \over (2\pi f_\pi)^6u^3}
\int_0^u\!\!dx\int_0^u\!\!dy \, {x^3(u^2-y^2)\over (1+4x^2)^2}(u-x)^2(2u+x) 
\bigg\{ 4x y\nonumber \\ && \times (1+4x^2+4y^2)-{1\over 4}\big[1+4(x+y)^2
\big] \big[1+4(x-y)^2\big] \ln{1+4(x+y)^2 \over 1+4(x-y)^2} \bigg\}\,, 
\end{eqnarray}
whereas the right 2-ring diagram produces a vanishing spin-trace. In the case of
the U-type 1-ring diagrams in Fig.\,3, three (out of four) Fermi sphere integrals 
over the pion-propagators and momentum-dependent interactions factorize by making use 
of tensor-contractions. For this reason one can represent the corresponding 
contribution to  $\bar E(\rho)$  as a one-parameter integral of the form:
\begin{equation} \bar E(\rho) = -{3g_A^4m_\pi^7 \over 4(4\pi f_\pi)^6u^3}
\int_0^u\!\!dx \, {G_V^2(x) \over x} \Big[2G_S(x)+G_T(x)\Big]\,. \end{equation}
Here, we have introduced the auxiliary functions: 
\begin{equation} G_V(x)= u(1+u^2+x^2)-{1\over 4x}\big[1+(u+x)^2\big]\big[
1+(u-x)^2\big]  \ln{1+(u+x)^2 \over 1+(u-x)^2}\,, \end{equation}
\begin{eqnarray} G_S(x)&=& {4u x\over 3}(2u^2 -3)+4x \big[ 
\arctan(u+x) + \arctan(u-x)\big]\nonumber \\ &&  +(x^2-u^2-1) \ln{1+(u+x)^2
\over 1+(u-x)^2} \,, \end{eqnarray}
\begin{eqnarray}G_T(x)&=& {u x \over 6}(8u^2 +3x^2)-{u \over 
2x}(1+u^2)^2+{1 \over 8}\bigg[{(1+u^2)^3 \over x^2}\nonumber \\ && -x^4  
+(1-3u^2)(1+u^2-x^2)\bigg]\ln{1+(u+x)^2 \over 1+(u-x)^2} \,, \end{eqnarray}
defined by vectorial and tensorial Fermi sphere integrals over a pion-propagator
(see also eq.(A.1) in ref.\,\cite{4body}). In the case of the Z-type 1-ring diagrams
in Fig.\,3, the two Fermi sphere integrals over the pion-propagators associated to 
the short pion-lines factorize via tensors. For the remaining two Fermi sphere 
integrals the angular part can be carried out, and in this procedure the third 
pion-propagator introduces the logarithmic function: 
\begin{equation} {\bf L}= \ln{1+(x+y)^2 \over 1+(x-y)^2} \,. \end{equation}
Putting all pieces together, one obtains the following double-integral representation 
for the contribution of the Z-type 1-ring diagrams from class II to the energy per 
particle of nuclear matter:
\begin{eqnarray} \bar E(\rho) &=& {3g_A^4m_\pi^7 \over (4\pi f_\pi)^6u^3}
\int_0^u\!\!dx\int_0^u\!\!dy \, G_V(x) \bigg\{G_S(y) \bigg[2y+{{\bf L} 
\over 2x}(x^2-y^2-1)\bigg] \nonumber \\ &&  + {G_T(y) \over 4y} \bigg[y^2-
3x^2-3+{{\bf L}\over 4xy} \Big(3(1+x^2)^2+2y^2-2x^2y^2-y^4\Big)\bigg]\bigg\}
\,. \end{eqnarray}
In the case of the X-type 1-ring diagrams in Fig.\,3 the assignment of momenta to 
the virtual pions is now such that the previous factorization approach does not 
work anymore. Only the Fermi sphere integral associated to the nucleon-line between 
the square-box and filled-circle vertex can be performed analytically. Hence, the 
contribution of the X-type 1-ring diagrams from class II to the energy per particle 
of nuclear matter result in the form:
\begin{equation}  \bar E(\rho) = {9g_A^4m_\pi^7 \over (4\pi f_\pi)^6 u^3}
\int\limits_{|\vec p_j|<u}\!\!\!{d^3p_1\, d^3p_2\,d^3p_3 \over (2\pi)^3} \, { G_V(|\vec 
\eta\, |)\over (1+ \vec q_1^{\,2}) (1+ \vec q_2^{\,2}) \, \vec \eta^{\,2}}
\, [(\vec q_1- \vec q_2)\cdot\vec q_2] \,\, (\vec q_1\cdot\vec \eta\,)
\,,  \end{equation}
with $\vec q_1 = \vec p_1-\vec p_3$, $\vec q_2 = \vec p_2-\vec p_3$ and 
$\vec \eta=\vec p_1+\vec p_2-\vec p_3$, where all momentum-vectors have been divided 
by the pion mass $m_\pi$. Note that due to rotational invariance it is sufficient to 
parametrize the nine-dimensional integral in eq.(18) by three radii ($0<p_{1,2,3}<u$), 
two directional cosines ($-1<z_{1,2}<1$), and one azimuthal angle ($0<\varphi<2\pi$).
   
\subsection{Class I}
Finally we come to class I, which is induced by a pair of ''contractions'' of two 
ordinary $\pi N$-vertices. The corresponding factorized expression for this 
$3\pi$-exchange $4N$-interaction reads: 
\begin{eqnarray}  V^a &=& {g_A^6 \over 16f_\pi^6} \, \vec\sigma_1\cdot \vec q_a \,
\tau_1^a\,  {1\over m_\pi^2+\vec q_a^{\,2}} \, \big[\epsilon^{abd} \, \tau_2^d \,
\vec q_a\cdot\vec q_b +\delta^{ab}\,\vec\sigma_2\cdot(\vec q_a\times\vec q_b)\big]\,
\nonumber \\ && \times {1\over (m_\pi^2+\vec q_b^{\,2})^2} \,\big[\epsilon^{bce} \, 
\tau_3^e \,\vec q_b\cdot\vec q_c +\delta^{bc}\,\vec\sigma_3\cdot(\vec q_b\times
\vec q_c)\big]\,{1\over m_\pi^2+\vec q_c^{\,2}} \,\vec \sigma_4\cdot \vec q_c \, 
\tau_4^c\,, \end{eqnarray}
Note that there is again an extra factor $2$ in comparison to eq.(3.33) in 
ref.\,\cite{evgeni4nb} due to the interchange symmetry $(12)\leftrightarrow(43)$ 
between pairs of nucleons. The evaluation of the left 2-ring diagram in Fig.\,2 
gives the following contributions to the energy per particle of isospin-symmetric 
nuclear matter:
\begin{eqnarray} \bar E(\rho) &=& {3g_A^6m_\pi^7 \over (2\pi f_\pi)^6u^3}
\int_0^u\!\!dx\int_0^u\!\!dy \,{x^3(u^2-y^2)\over (1+4x^2)^2}(u-x)^2 (2u+x)
\nonumber \\ && \times \bigg\{-2x y(1+20x^2+4y^2)+32x^3\big[\arctan(2x+2y)-
\arctan(2x-2y)\big]\nonumber \\ && +\Big(10x^4-5x^2+y^2-12x^2y^2+2y^4+
{1\over 8}\Big) \ln{1+4(x+y)^2\over 1+4(x-y)^2}\bigg\}\,, \end{eqnarray}
and to the energy per particle of pure neutron matter:
\begin{eqnarray} \bar E_n(\rho_n) &=& {g_A^6m_\pi^7 \over (2\pi f_\pi)^6u^3}
\int_0^u\!\!dx\int_0^u\!\!dy \,{x^3(u^2-y^2)\over (1+4x^2)^2}(u-x)^2 (2u+x)
\nonumber \\ && \times \bigg\{x y(1-12x^2+4y^2)+8x^3\big[\arctan(2x+2y)-
\arctan(2x-2y)\big]\nonumber \\ && +\Big(3x^4-2x^2y^2-y^4-{3x^2+y^2 
\over 2}-{1 \over 16}\Big) \ln{1+4(x+y)^2 \over 1+4(x-y)^2}\bigg\}\,. 
\end{eqnarray} 
The neutron density $\rho_n$ is related to the neutron Fermi momentum $k_n$ by 
$\rho_n=k_n^3/3\pi^2$. In eq.(21) and all following formulas for $\bar E_n(\rho_n)$ 
the dimensionless variable $u$ has the meaning $u  =k_n/m_\pi$. Next one notices 
that for class I the right 2-ring diagram in Fig.\,2 produces vanishing spin-trace. 

Furthermore, one obtains from the U-type 1-ring diagram in Fig.\,3 the following 
contributions:
\begin{equation} \bar E(\rho) = {g_A^6m_\pi^7 \over 4(4\pi f_\pi)^6u^3}
\int_0^u\!\!dx\,\bigg\{ \Big[8G_S^2(x)+G_T^2(x)\Big]K_S(x)+ \Big[8G_S(x)
+G_T(x)\Big]G_T(x) K_T(x) \bigg\} \,, \end{equation}
\begin{equation} \bar E_n(\rho_n) = {g_A^6m_\pi^7 \over 12(4\pi f_\pi)^6u^3}
\int_0^u\!\!dx \,\bigg\{ \Big[2G_S^2(x)+G_T^2(x)\Big]K_S(x)+ \Big[G_T(x)-4G_S(x)
\Big]G_T(x) K_T(x) \bigg\} \,, \end{equation}
where we have introduced two new auxiliary functions:
\begin{equation} K_S(x)= 2u-3\arctan(u+x)-3\arctan(u-x) + 
 {2+u^2-x^2\over 2x}  \ln{1+(u+x)^2 \over 1+(u-x)^2} \,,\end{equation} 
\begin{equation} K_T(x)= {u\over 4x^2}(3+3u^2-x^2)+{1\over 16 x^3}\Big[
x^2(x^2+2u^2-2)-3(1+u^2)^2\Big]\ln{1+(u+x)^2 \over 1+(u-x)^2}\,, \end{equation}
defined by a tensorial Fermi sphere integral over a squared pion-propagator.

On the other hand the Z-type 1-ring diagrams give the following contributions:
\begin{eqnarray} \bar E(\rho) &=& {g_A^6m_\pi^7 \over (4\pi f_\pi)^6u^3}
\int_0^u\!\!dx\int_0^u\!\!dy\, \bigg\{2 G_S(x)G_S(y)\,[{\bf D}-{\bf L}]+G_S(x) G_T(y) 
\nonumber \\ && \times \bigg[2{\bf D}-{6x \over y}+{{\bf L}\over 2y^2}(3+3x^2-y^2)
\bigg]+ G_T(x) G_T(y) \nonumber \\ && \times \bigg[{{\bf D}\over 2} +{3y\over 8x}
+{27\over 16x y} -{{\bf L}\over 64 x^2y^2}(27+60y^2+2x^2y^2+6y^4)\bigg]\bigg\} \,,
\end{eqnarray} 
\begin{eqnarray} \bar E_n(\rho_n) &=& {g_A^6m_\pi^7 \over 6(4\pi f_\pi)^6u^3}
\int_0^u\!\!dx\int_0^u\!\!dy \,\bigg\{G_S(x)G_S(y)\, [{\bf D}-{\bf L}] +G_S(x) G_T(y) 
\nonumber \\ && \times \bigg[{6x \over y}-2{\bf D} +{{\bf L}\over 2y^2}(y^2-3x^2-3)
\bigg]+G_T(x)G_T(y) \nonumber \\ && \times\bigg[{\bf D}+{3y\over 4x}+{27\over 8x y} 
-{{\bf L}\over 32 x^2y^2}(27+60y^2+2x^2y^2+6y^4)\bigg]\bigg\} \,,\end{eqnarray} 
with the rational function:
\begin{equation}{\bf D} = {1 \over 1+(x-y)^2}-{1 \over 1+(x+y)^2}\,, \end{equation}
arising from an angular integral $\int_{-1}^1\!dz\dots$ over a squared 
pion-propagator.

Finally, the X-type 1-ring diagram leads to contributions of the form:
\begin{eqnarray}  \bar E(\rho) &=& {3g_A^6m_\pi^7 \over (4\pi f_\pi)^6 u^3}
\int\limits_{|\vec p_j|<u}\!\!\!{d^3p_1\, d^3p_2\,d^3p_3\over(2\pi)^3}\,{1\over(1+ 
\vec q_1^{\,2}) (1+ \vec q_2^{\,2})}\nonumber \\ && \times\bigg\{ 
3\vec q_1^{\,2}\,\vec q_2^{\,2} \Big[K_S(|\vec \eta\,|) -2 K_T(|\vec \eta\,|)\Big] +
(\vec q_1\cdot\vec q_2)^2 \Big[7K_T(|\vec \eta\,|)-K_S(|\vec \eta\,|)\Big]
\nonumber \\ && +\,{3 \over \vec \eta^{\,2}} K_T(|\vec \eta\,|) \Big[  6\,\vec q_1^{\,2}
\,(\vec \eta\cdot\vec q_2)^2- 7(\vec q_1\cdot\vec q_2)\,(\vec \eta\cdot \vec q_1)
\, (\vec \eta\cdot\vec q_2)\Big] \bigg\}\,, \end{eqnarray}
\begin{eqnarray}  \bar E_n(\rho_n) &=& {g_A^6m_\pi^7 \over (4\pi f_\pi)^6 u^3}
\int\limits_{|\vec p_j|<u}\!\!\!{d^3p_1\, d^3p_2\,d^3p_3\over(2\pi)^3}\,{1\over(1+\vec 
q_1^{\,2}) (1+ \vec q_2^{\,2})}\nonumber \\ && \times \bigg\{ 
3(\vec q_1\cdot \vec q_2)^2\Big[K_S(|\vec \eta\,|) +K_T(|\vec \eta\,|)\Big] 
- \vec q_1^{\,2}\,\vec q_2^{\,2} \Big[K_S(|\vec \eta\,|)+2K_T(|\vec \eta\,|)\Big]
\nonumber \\ && +\,{3 \over \vec \eta^{\,2}} K_T(|\vec \eta\,|) \Big[  
2\,\vec q_1^{\,2}\,(\vec \eta\cdot \vec q_2)^2- 3(\vec q_1\cdot\vec q_2)\,(\vec 
\eta\cdot\vec q_1) \, (\vec \eta\cdot\vec q_2)\Big] \bigg\}\,, \end{eqnarray}
with $\vec q_1 = \vec p_1-\vec p_3$, $\vec q_2 = \vec p_2-\vec p_3$ and $\vec \eta= 
\vec p_1+\vec p_2-\vec p_3$. 

This completes our presentation of semi-analytical results for the energy per 
particle of isospin-symmetric nuclear matter $\bar E(\rho)$ and pure neutron matter 
 $\bar E(\rho)$ as derived from the reducible chiral $4N$-interactions. The given 
formulas allow for an easy and accurate numerical evaluation  \cite{code} of these 
leading-order chiral four-body correlations.  

\section{Results and discussion}
We are now in the position to present numerical results. The physical parameters 
related to pion-exchange are: $g_A = 1.29$ (nucleon axial-vector coupling constant), 
$f_\pi = 92.4$\,MeV (pion decay constant), and $m_\pi = 135\,$MeV (neutral pion mass).
For the spin-spin contact-coupling $C_T$ we choose two opposite values, $C_T=
0.22\,$fm$^2$ and $C_T=-0.45\,$fm$^2$, taken from table I in ref.\cite{darmstadt}.

The contributions to the energy per particle $\bar E(\rho)$ of isospin-symmetric 
nuclear matter are presented in Fig.\,4 for the small $C_T$-value, $C_T=0.22\,
$fm$^2$, in the density region $\rho_0/4<\rho <2 \rho_0$, with $\rho_0=0.16\,
$fm$^{-3}$ the empirical saturation density. The dashed lines (with appropriate 
labels) show separately the five classes (VII, V, IV, II, I) and the full line gives 
their total sum. One observes that class VII (proportional to $C_T^2$) is smallest 
in magnitude and that other classes (V+IV and II+I) cancel each other to a large 
extent. The largest repulsive and attractive contributions are provided by the 
$C_T$-independent classes II and I, respectively. As a consequence of these 
compensations the net attraction stays above $-1.3\,$MeV for all densities 
$\rho<0.36\,$fm$^{-3}$. The analogous results for the larger negative $C_T$-value, 
$C_T=-0.45\,$fm$^2$, are shown in Fig.\,5. One notices that class VII ($\sim C_T^2$) 
is still smallest in magnitude and that the classes V and IV have changed their 
position in comparison to the sequence in Fig.\,4. The net attraction reaches now a 
value of $-1.7\,$MeV at $\rho=0.36\,$fm$^{-3}$. Exploring the whole range $|C_T|
<0.5\,$fm$^2$ of the spin-spin contact-coupling $C_T$, one finds that the 
$4N$-attraction in nuclear matter does not exceed values of $-1.3\,$MeV for 
densities $\rho <2 \rho_0$.


\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.5,clip]{4Nreducib22.eps}
\end{center}
\vspace{-.6cm}
\caption{Reducible chiral $4N$-interactions in isospin-symmetric nuclear matter. 
The dashed lines show the individual contributions from the five classes for 
$C_T=0.22$\,fm$^2$, and the full line gives their total sum.}
\end{figure}


\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.5,clip]{4Nreducib45.eps}
\end{center}
\vspace{-.6cm}
\caption{Reducible chiral $4N$-interactions in isospin-symmetric nuclear matter. 
The dashed lines show the individual contributions from the five classes for 
$C_T=-0.45$\,fm$^2$, and the full line gives their total sum.}
\end{figure}

\begin{figure}[h!]
\begin{center}
\includegraphics[scale=0.5,clip]{4Nredneutron.eps}
\end{center}
\vspace{-.6cm}
\caption{Reducible chiral $4N$-interactions in pure neutron matter. The
conributions from 2-ring and 1-ring diagrams are shown separately. The U-type, 
Z-type, and X-type diagrams cancel each other almost completely.}
\end{figure}

In the case of pure neutron matter only class I (see subsection 2.5) contributes and 
the result for the energy per particle $\bar E_n(\rho_n)$ comes out parameterfree.
Fig.\,6 shows by the dashed lines the contributions from the 2-ring diagram and the 
sum of the three 1-ring diagrams in the density region $0.04\,$fm$^{-3}<\rho_n<0.4
\,$fm$^{-3}$. One recognizes that the latter part is extremely small. This feature 
has its origin in the almost complete compensation between the contributions from 
the U-type, Z-type and X-type 1-ring diagrams, which are shown by the three 
dashed-dotted lines in Fig.\,6. Interestingly, the net result for  $\bar E_n(\rho_n)$ 
is thus entirely determined by the 2-ring part written in eq.(21). The corresponding 
attractive energy per particle goes approximately as $\rho_n^{7/3}$ and it reaches 
the value of about $-2.1$\,MeV at the (relatively high) neutron density of $\rho_n=
0.4$\,fm$^{-3}$. 

In passing we note that our numerical results for the chiral four-body correlations 
in nuclear and neutron matter agree perfectly \cite{thomas} with those of the 
Darmstadt group  \cite{darmstadt} if the (unessential) regulator-function $f_R$ is 
set to $f_R=1$.  

In summary one can conclude that the chiral four-nucleon correlations studied in 
this work are at least one order of magnitude smaller than those provided by the 
strongly coupled $\pi N\Delta$-system \cite{4body} with its small mass-gap of  
$293$\,MeV\,$\simeq 1.1 k_{f0}$.  

\section*{Acknowledgements}
We thank E. Epelbaum, T. Kr\"uger and A. Schwenk for informative discussions. 



