\section{Introduction}

The most important problem of comtemporary fundamental theoretical physics is the unification of quantum theory
and general relativity. A quantum theoretical setting of general relativity as usual quantum field theory
in analogy to the description of the other fundamental interactions is not possible, since such a description
leads to divergencies, which cannot be cured. This fact seems to indicate that the problem of incorporating
general relativity to a quantum theoretical description of all fundamental interactions is directly connected
to the question, whether the basic postulates of quantum field theory in its usual setting can be transferred
to such a more fundamental description of nature. These basic postulates are related to the idea of the basic
structure of space-time. General relativity is a background independent theory, which in contrast to the
usual setting of quantum field theory, on which the other fundamental interactions of the standard model of
particle physics are based, does not presuppose a nondynamical metric structure of space-time. Quantum theory
on the other hand contains a constitutive element of nonlocality, which becomes manifest in the famous
Einstein-Podolsky-Rosen thought experiment for example. The basic postulates of general quantum theory itself,
this means the general Hilbert space formalism in the setting of Dirac without relation to additional field
theoretic elements, does not contain special assumptions about a space-time structure.
Therefore it seems to be possible that a quantum theoretical description of general relativity as well as a
unification of all fundamental interactions appearing in nature could indeed require a decisive modification
of the postulates usually presupposed in quantum field theory. This holds especially for the basic assumption
of local commutativity or microscopic causality, which implies that the commutator or anticommutator
respectively of field operators at two separated space-time points vanishes.

In accordance with the above explanations in the present paper is suggested a generalization of the quantization
principle as it is presupposed in quantum field theory and as it is transferred to canonical quantum gravity.
Usually, the commutator or anticommutator respectively between the field and the corresponding canonical conjugated
momentum is assumed to be proportional to the delta function. In the present paper this quantization principle is
extended in such a way that the delta function is replaced by a general function, which depends on the distance
between the space-time points the corresponding operators refer to. This implies that the commutator between
the field operator and the corresponding operators describing the canonical conjugated quantity on two separated
space-time points does not vanish anymore. Such a quantization principle establishes a nonlocal structure of the
quantum description of a theory, violates the posulate of local commutativity and thus can be called nonlocal
quantization principle. A natural assumption for the function replacing the delta function is a Gaussian function,
which as limiting case merges the delta function. It will be shown that under this assumption one obtains a
propagator, which corresponds to the propagator of usual quantum field theory multiplied with an exponential
damping factor depending on the squared three momentum. This implies that this generalized quantum field theory
contains no ultra violet divergences. The obtained propagator corresponds to the propagator within the coherent
state approach to noncommutative geometry, which has been developed in \cite{Smailagic:2003rp},\cite{Smailagic:2003yb},\cite{Smailagic:2004yy}, extended in
\cite{Huang:2008yy},\cite{Kober:2010um} and transferred to the gravitational field in
\cite{Kober:2011am}. In \cite{Modesto:2009qc} also appears a propagator of this structure.
It is also possible to formulate this generalized quantization principle within the path integral formalism.
Such a formulation is obtained by replacing the usual inner product between field eigenstates and eigenstates
of the canonical conjugated variables by the generalized one according to the nonlocal quantization principle.
This results in a generalized path integral expression. If the nonlocal quantization principle is considered
as a fundamental description of nature, then one has of course to assume that it is also valid for the
quantization of the gravitational field. Accordingly it has to be postulated a nonlocal generalization of the commutator
between the three metric and its canonical conjugated momentum in quantum geometrodynamics. With respect to the
quantization of the gravitational field, the intricacy arises that the distance between two space-time points
depends on the metric field. Therefore the commutator also depends on the three metric itself. A Generalized
quantization principle in canonical quantum gravity and quantum cosmology has already been considered in
\cite{Kober:2011uj}. Other considerations of a generalized quantization principle with respect to the
variables of quantum cosmology and black holes can be found in \cite{Majumder:2011ad},\cite{Majumder:2011bv}.
The considerations in \cite{Kober:2011uj},\cite{Majumder:2011ad},\cite{Majumder:2011bv} have transferred the
generalized uncertainty principle in quantum mechanics, see \cite{Maggiore:1993kv},\cite{Kempf:1994su},\cite{Hinrichsen:1995mf} for example, which leads to a minimal length, to the variables describing the gravitational field. Another
generalization of the quantization principle in canonical quantum gravity, which is based on quaternions
and has also been transferred to $N=1$ supergravity, has been developed in \cite{Kober:2015bkv}.
However, these generalizations of the quantization principle in canonical quantum gravity in contrast to the
approach of the present paper remain local. Nonlocality has been considered in the context of quantum field theory
\cite{Yukawa:1950eq},\cite{Yukawa:1950er},\cite{Alebastrov:1973vw},\cite{Marnelius:1974rq},\cite{Wreszinski:1975yz},\cite{Dubnickova:1976jb},\cite{Bandyopadhyay:1976sp},\cite{Buchholz:1978fv},\cite{Efimov:1980qp},\cite{Buchholz:1985xs},\cite{Mannheim:1986ad},\cite{Marino:1987tk},\cite{Mironov:1989hm},\cite{Marino:1989ve},\cite{Bernard:1990ys},\cite{Cornish:1991kg},\cite{Cornish:1991ja},\cite{Barci:1995ad},\cite{Solovev:1998gm},\cite{Esposito:1999hp},
\cite{Amorim:1999mr},\cite{Shirokov:2001zh},\cite{Denk:2003jj},\cite{Jain:2003xs},\cite{Denk:2004pk},\cite{Gorbar:2005sx},\cite{Novello:1994ym},\cite{Novello:1994zi},\cite{Wetterich:1997bz},\cite{Dobado:1997yh},\cite{Ahluwalia:1999yc},
general relativity
\cite{Chen:2001gx},\cite{Mukohyama:2001jv},\cite{Bergman:2001rw},\cite{Barvinsky:2003kg},\cite{Hehl:2008eu},\cite{Koshelev:2008ie},\cite{Bronnikov:2009az},\cite{Calcagni:2010ab},\cite{Blome:2010xn},\cite{Chicone:2011me},\cite{Barvinsky:2011rk},\cite{Elizalde:2011su}
and quantum gravity
\cite{Rayner:1989wd},\cite{Cornish:1991iz},\cite{Kim:1993iua},\cite{Barvinsky:1994hw},\cite{Esposito:1998hp},\cite{Shojai:1999jj},\cite{Dominguez:1999sq},\cite{Moffat:2002ce},\cite{Reuter:2002kd},\cite{vonBorzeszkowski:2002ri},\cite{Reuter:2003yb},\cite{Giddings:2006be},\cite{Barvinsky:2011hd}
and also of black holes \cite{Nicolini:2012eu}. But the concept of a nonlocal qusantization principle
has not been studied yet.

The paper is structured as follows:
The generalized formalism of quantum field theory based on this nonlocal quantization principle for a scalar field
is considered at the beginning. Subsequently, the corresponding generalized expression for the propagator is calculated.
Especially the mentioned special case is regarded, where the function the commutator depends on is assumed to be a
Gaussian function leading to a propagator containing an exponential damping function, which depends on the squared
three momentum with negative sign of the exponent as additional factor. This means that the corresponding quantum
field theory contains no ultra violet divergences. After this the nonlocal quantization principle is considered in
the context of canonical quantum gravity and the generalized representations of the operators as well as the
generalized constraints are formulated. At the end the nonlocal quantization principle is transferred to the
path integral formula.

\section{Nonlocal Quantization Principle}

In this section is presented the basic idea of the nonlocal quantization principle in a general way.
According to the explanations given in the introduction a generalization of the fundamental quantization
principle with respect to the quantization of fields is suggested. This quantization principle omits the
postulate of microscopic causality. The generalized nonlocal quantization principle is obtained by the
following transition:

\begin{equation}
\left[\hat \Phi(x),\hat \Pi(y)\right]=i\delta^3(x-y)\quad \longrightarrow \quad
\left[\hat \Phi(x),\hat \Pi(y)\right]=i\alpha^3(x-y),
\label{generalization_commutator}
\end{equation}
where $\hat \Phi(x)$ denotes a field operator, $\hat \Pi(x)$ denotes the corresponding canonical conjugated
operator, $x$ and $y$ denote three-vectors, $\delta^3(x-y)$ denotes the three-dimensional delta function
referring to a spacelike hypersurface and $\alpha^3(x-y)$ denotes a general function on the spacelike
hypersurface, which is symmetric in $x$ and $y$ and thus just depends on the distance of the points
$x$ and $y$. This assumption is natural, if homogenity and isotropy of space shall be maintained.
The field operator $\hat \Phi(x)$ and the corresponding canonical conjugated operator $\hat \Pi(x)$
are generalizations of the usual operators, which shall be denoted with $\hat \phi(x)$ and $\hat \pi(x)$
and obey the usual commutation relations:

\begin{equation}
\left[\hat \phi(x),\hat \pi(y)\right]=i\delta^3(x-y).
\end{equation}
The generalized operators can be represented in a specific way by using the usual operators.
One possibility is the representation, in which the generalized field operator $\Phi(x)$ is equal
to $\phi(x)$ and the generalized canonical conjugated operator is modified,

\begin{equation}
\hat \Phi(x)=\hat \phi(x),\quad \hat \Pi(x)=\int d^3 z\ \alpha^3(x-z)\ \hat \pi(z),
\end{equation}
and another possibility is the representation, where the generalized canonical conjugated operator
is equal to $\pi(x)$ and the generalized field operator is modified,

\begin{equation}
\hat \Phi(x)=\int d^3 z\ \alpha^3(x-z)\ \hat \phi(z),\quad \hat \Pi(x)=\hat \pi(x).
\end{equation}
The momentum eigenstates in the representation, where the momentum operator is generalized and
where the usual operators are represented with respect to the field, looks as following:

\begin{equation}
|\Pi[\Phi(x)]\rangle=\exp\left[i\int d^3 x\ d^3 y\ {\bar \alpha}^3(x-y)\Phi(x)\Pi(y)\right],
\label{eigenstate_Pi}
\end{equation}
where $\bar \alpha^3(x-y)$ is defined by the following relation:

\begin{equation}
\int d^3 z\ \alpha^3(x-z){\bar \alpha}^3(z-y)=\delta(x-y),
\end{equation}
and $\Pi(x)$ denotes the momentumn eigenvalue corresponding to the momentum eigenstate.
This can be seen as follows:

\begin{eqnarray}
\hat \Pi(x)|\Pi\rangle&=&-i\int d^3 y\ \alpha^3(x-y)\frac{\delta}{\delta \Phi(y)}
\exp\left[i\int d^3 w\ d^3 z\ {\bar \alpha}^3(w-z)\Phi(w)\Pi(z)\right]\nonumber\\
&=&\int d^3 y\ d^3 w\ d^3 z\ \alpha^3(x-y){\bar \alpha}^3(w-z)\delta(w-y)\Pi(z)
\exp\left[i\int d^3 w\ d^3 z\ {\bar \alpha}^3(w-z)\Phi(w)\Pi(z)\right]\nonumber\\
&=&\int d^3 y\ d^3 z\ \alpha^3(x-y){\bar \alpha}^3(y-z)\Pi(z)
\exp\left[i\int d^3 w\ d^3 z\ {\bar \alpha}^3(w-z)\Phi(w)\Pi(z)\right]\nonumber\\
&=&\int\ d^3 z\ \delta(x-z)\Pi(z)\exp\left[i\int d^3 w\ d^3 z\ {\bar \alpha}^3(w-z)\Phi(w)\Pi(z)\right]\nonumber\\
&=&\Pi(x)\exp\left[i\int d^3 w\ d^3 z\ {\bar \alpha}^3(w-z)\Phi(w)\Pi(z)\right]\nonumber\\
&=&\Pi(x)|\Pi\rangle.
\label{eigenwert_equation_Pi}
\end{eqnarray}
The corresponding position eigenstate is given by

\begin{equation}
|\Phi[\Phi(x)]\rangle=\delta\left(\Phi[\Phi(x)]-\Phi^{\prime}[\Phi(x)]\right),
\label{eigenstate_Phi}
\end{equation}
where $\Phi(x)$ denotes the eigenvalue to the position eigenstate. The inner product for two states
referring to the field can be defined as

\begin{eqnarray}
\langle \Psi[\Phi(x)]|\Omega[\Phi(x)]\rangle=\int d \Phi(x)\ d^3 x\ \Psi[\Phi(x)]\ \Omega[\Phi(x)].
\end{eqnarray}
This implies that the inner product of a momentum eigenstate with a position eigenstate is given by

\begin{equation}
\langle \Phi|\Pi \rangle=\exp\left[i\int d^3 w\ d^3 z\ {\bar \alpha}^3(w-z)\Phi(w)\Pi(z)\right].
\label{innerproduct_PhiPi}
\end{equation}

\section{Derivation of the Propagator in Quantum Field Theory with Nonlocal Quantization Principle}

The general concept of the nonlocal quantization principle can now be studied in the context of a scalar
field theory. In this section the corresponding propagator is calculated in its general shape meaning that
no special assumption for the special choice of the quantization function is made so far. To determine the
propagator, it is necessary to formulate the generalized expression for the plane wave expansion of the
scalar field. The following expression for the quantum field

\begin{equation}
\hat \Phi(x,t)=\int \frac{d^3 p}{\sqrt{(2\pi)^3 2\omega_p}}\ \left(a(p)\int d^3 z\ \alpha^3(x-z) e^{i{\bf pz}-i\omega_p t}
+a^{\dagger}(p)e^{-i{\bf px}+i\omega_p t}\right)
\label{generalized_field}
\end{equation}
and its canoniocal conjugated variable

\begin{equation}
\hat \Pi(x,t)=\int \frac{d^3 p\ i\omega_p}{\sqrt{(2\pi)^3 2\omega_p}}\ \left(a(p)
\int d^3 z\ \alpha^3(x-z) e^{i{\bf pz}-\omega_p t}-a^{\dagger}(p)e^{-i{\bf px}+i\omega_p t}\right),
\end{equation}
where $a(p)$ and $a^{\dagger}(p)$ are the usual creation and annihilation operators
fulfilling the commutation relations

\begin{equation}
\left[a(p),a^{\dagger}(p^{\prime})\right]=\delta(p-p^{\prime}),
\end{equation}
obey the generalized commutation relations given in ($\ref{generalization_commutator}$). This can be shown
as follows:

\begin{eqnarray}
\left[\hat \Phi(x,t),\hat \Pi(y,t)\right]&=&\frac{i}{(2\pi)^3}\int \frac{d^3 p}{\sqrt{2\omega_p}}
\ \int \frac{d^3 p^{\prime} \omega_{p^{\prime}}}{\sqrt{2\omega_p^{\prime}}} \int d^3 z
\left[\left[a(p),a^{\dagger}(p^{\prime})\right]\alpha^3(x-z)e^{i{\bf py}-i{\bf p^{\prime}z}}
+\left[a(p^{\prime}),a^{\dagger}(p)\right]\alpha^3(y-z)e^{i{\bf p^{\prime}x}-i{\bf pz}}\right]\nonumber\\
&=&\frac{i}{2(2\pi)^3}\int \frac{d^3 p}{\sqrt{\omega_p}}\ \int \frac{d^3 p^{\prime} \omega_p^{\prime}}
{\sqrt{\omega_p^{\prime}}} \int d^3 z \left[\delta^3(p-p^{\prime})\alpha^3(x-z)e^{i{\bf py}-i{\bf p^{\prime}z}}
+\delta^3(p^{\prime}-p)\alpha^3(y-z)e^{i{\bf p^{\prime}x}-i{\bf pz}}\right]\nonumber\\
&=&\frac{i}{2(2\pi)^3}\int d^3 p\ \int d^3 z \left[\alpha^3(x-z)e^{i{\bf py}-i{\bf pz}}
+\alpha^3(y-z)e^{i{\bf px}-i{\bf pz}}\right]\nonumber\\
&=&\frac{i}{2}\int d^3 z \left[\alpha^3(x-z)\delta^3(y-z)+\alpha^3(y-z)\delta^3(x-z)\right]\nonumber\\
&=&\frac{i}{2}\left[\alpha^3(x-y)+\alpha^3(y-x)\right]\nonumber\\
&=&i\alpha^3(x-y).
\end{eqnarray}
The corresponding propagator of the generalized quantum field ($\ref{generalized_field}$) is obtained by
inserting it to the genereal expression for the propagator in terms of the quantum field,
which is given by

\begin{equation}
G(x-y)=\langle 0|T\left[\hat \Phi(x)\hat \Phi(y)\right]|0 \rangle,
\end{equation}
where $T$ denotes the time ordering operator $|0 \rangle$ denotes the vacuum state.
This leads to the following expression:

\begin{eqnarray}
G(x-y)&=&\int \frac{d^3 p}{(2\pi)^3 2p_0}\ \int d^4 z
\left[\theta(x_0-y_0)\alpha^4(x-z)e^{-ipz+ipx}+\theta(y_0-x_0)\alpha^4(y-z)e^{-ipz+ipy}\right],
\end{eqnarray}
where the function $\alpha^4(x-y)$ has been defined as follows:

\begin{equation}
\alpha^4(x-y)=\alpha^3(x-y)\delta(x_0-y_0).
\end{equation}
By following the usual procedure of expresssing the $\theta$-function by an integral

\begin{equation}
\theta(t)=\lim_{\epsilon \to 0} -\frac{1}{2\pi i}\int dE \frac{e^{-iEt}}{E+i\epsilon},
\end{equation}
the expression of the propagator can be transformed to

\begin{eqnarray}
G(x-y)=i\int \frac{d^3 p\ dE}{(2\pi)^4 2\omega_p}\frac{1}{E+i\epsilon}
\int d^4 z\ \left[\alpha^4(x-z)e^{i{\bf p}({\bf z}-{\bf y})-i(\omega_p+E)(z_0-y_0)}
+\alpha^4(y-z)e^{i{\bf p}({\bf z}-{\bf x})-i(\omega_p+E)(z_0-x_0)}\right].\nonumber\\
\end{eqnarray}
Changing of the variables according to,

\begin{equation}
E^{\prime}=E+\omega_p,
\end{equation}
leads to

\begin{eqnarray}
G(x-y)&=&i\int \frac{d^3 p\ dE^{\prime}}{(2\pi)^4 2\omega_p}\frac{1}{E^{\prime}-\omega_p+i\epsilon}
\int d^4 z\ \left[\alpha^4(x-z) e^{i{\bf p}({\bf z}-{\bf y})-iE^{\prime}(z_0-y_0)}
+\alpha^4(y-z) e^{i{\bf p}({\bf z}-{\bf x})-iE^{\prime}(z_0-x_0)}\right]\nonumber\\
&=&i\int \frac{d^3 p\ dE^{\prime}\ d^4 z}{2(2\pi)^4}\left[\frac{e^{i{\bf p}({\bf z}-{\bf y})
-iE^{\prime}(z_0-y_0)}}{\omega_p E^{\prime}-\omega^2_p+i\epsilon}\alpha^4(x-z)+x \longleftrightarrow y\right]\nonumber\\
&=&i\int \frac{d^3 p\ dE^{\prime}\ d^4 z}{2(2\pi)^4}\left[\frac{e^{i{\bf p}({\bf z}-{\bf y})-iE^{\prime}(z_0-y_0)}}
{\sqrt{{\bf p}^2+m^2}E^{\prime}-{\bf p}^2-m^2+i\epsilon}
\alpha^4(x-z)+x \longleftrightarrow y\right].
\end{eqnarray}
By setting

\begin{equation}
p_0=E^{\prime}
\end{equation}
one obtains the following generalized expression for the propagator in position space:

\begin{eqnarray}
G(x-y)&=&i\int \frac{d^3 p\ dp_0\ d^4 z}{2(2\pi)^4}
\left[\frac{e^{i{\bf p}({\bf z}-{\bf y})-ip_0 (z_0-y_0)}}{\sqrt{{\bf p}^2+m^2}p_0-{\bf p}^2-m^2+i\epsilon}
\alpha^4(x-z)+x \longleftrightarrow y\right]\nonumber\\
&=&i\int \frac{d^4 p\ d^4 z}{2(2\pi)^4}\left[\frac{e^{ip(z-y)}}{p^2-m^2+i\epsilon}
\alpha^4(x-z)+x \longleftrightarrow y\right].
\label{generalized_propagator}
\end{eqnarray}

\section{Special Model for the Nonlocal Quantization Principle}

The generalized nonlocal quantum field theory shall now be considered for a special choice of the quantization
function $\alpha^3(x-y)$ appearing in ($\ref{generalization_commutator}$). Since this function should be
symmetric in $x$ and $y$, should become smaller, if the distance becomes larger, should vanish at infinity
and should merge the delta function in the limiting case, a natural choice seems to be a Gaussian function,

\begin{equation}
\alpha^3(x-y)=\frac{e^{-\frac{1}{2}\frac{({\bf x}-{\bf y})^2}{\Delta}}}{\Delta\sqrt{2\pi}},
\label{Gaussian_function_3}
\end{equation}
where $\Delta$ describes the width of the Gaussian function. Accordingly $\alpha^4(x-y)$ reads

\begin{equation}
\alpha^4(x-y)=\delta(x_0-y_0)\frac{e^{-\frac{1}{2}\frac{({\bf x}-{\bf y})^2}{\Delta}}}{\Delta\sqrt{2\pi}}.
\label{Gaussian_function_4}
\end{equation}
If the expression for $\alpha^4(x-y)$ given in ($\ref{Gaussian_function_4}$) is inserted into
($\ref{generalized_propagator}$), one obtains the expression for the corresponding special
propagator

\begin{equation}
G(x-y)=i\int \frac{d^4 p\ d^4 z}{2(2\pi)^4}\left[\frac{e^{ip(z-y)}}{p^2-m^2+i\epsilon}\delta(x_0-z_0)
\frac{e^{-\frac{1}{2}\frac{({\bf x}-{\bf z})^2}{\Delta}}}{\Delta\sqrt{2\pi}}+x \longleftrightarrow y\right].
\end{equation}
By using the substitution

\begin{equation}
w=x-z,
\end{equation}
this expression can be transformed to

\begin{eqnarray}
G(x-y)&=&i\int \frac{d^4 p\ d^4 w}{2(2\pi)^4}\left[\frac{e^{ip(x-w-y)}}{p^2-m^2+i\epsilon}\delta(w_0)
\frac{e^{-\frac{1}{2}\frac{{\bf w}^2}{\Delta}}}{\Delta\sqrt{2\pi}}+x \longleftrightarrow y\right]\nonumber\\
&=&i\int \frac{d^4 p\ d^3 w\ dw_0}{2(2\pi)^4}\frac{1}{\Delta\sqrt{2\pi}}\left[\frac{e^{ip(x-y)}}{p^2-m^2+i\epsilon}\delta(w_0)
e^{-\frac{1}{2}\frac{{\bf w}^2}{\Delta}+i{\bf p}{\bf w}-ip_0 w_0}+x \longleftrightarrow y\right]\nonumber\\
&=&i\int \frac{d^4 p\ d^3 w}{2(2\pi)^4}\frac{1}{\Delta\sqrt{2\pi}}\left[\frac{e^{ip(x-y)}}{p^2-m^2+i\epsilon}
e^{-\frac{1}{2}\frac{{\bf w}^2}{\Delta}+i{\bf p}{\bf w}}+x \longleftrightarrow y\right].
\end{eqnarray}
By substituting {\bf w} through {\bf a} according to

\begin{equation}
{\bf a}=\frac{{\bf w}}{\sqrt{2\Delta}}-i\sqrt{\frac{\Delta}{2}}{\bf p},
\end{equation}
the corresponding integral can be solved and one obtains the final expression for the
propagator represented in position space of the quantum field theory with the nonlocal
quantization principle under consideration of the special choice ($\ref{Gaussian_function_3}$)
for the quantization function

\begin{eqnarray}
G(x-y)&=&i\int \frac{d^4 p\ d^3 a}{(2\pi)^4}\frac{\Delta}{\sqrt{\pi}}
e^{-\frac{\Delta}{2}{\bf p}^2}\left[\frac{e^{ip(x-y)}}{p^2-m^2+i\epsilon}
e^{-{\bf a}^2}+x \longleftrightarrow y\right]\nonumber\\
&=&i\int \frac{d^4 p}{(2\pi)^4}\pi\Delta
e^{-\frac{\Delta}{2}{\bf p}^2}\left[\frac{e^{ip(x-y)}}{p^2-m^2+i\epsilon}
+x \longleftrightarrow y\right]\nonumber\\
&=&i\pi \Delta \int \frac{d^4 p}{(2\pi)^4} \left[\frac{e^{ip(x-y)}+e^{ip(y-x)}}
{p^2-m^2+i\epsilon}\right]e^{-\frac{\Delta}{2}{\bf p}^2}\nonumber\\
&=&i\Delta \int \frac{d^4 p}{(2\pi)^3}\frac{e^{ip(x-y)}
e^{-\frac{\Delta}{2}{\bf p}^2}}{p^2-m^2+i\epsilon}.
\end{eqnarray}
The corresponding expression for the propagator represented in momentum space is accordingly
given by

\begin{equation}
G(p)=\frac{i2\pi \Delta e^{-\frac{\Delta}{2}{\bf p}^2}}{p^2-m^2+i\epsilon}.
\end{equation}
This means that this quantum field theory in analogy to the situation in noncommutative geometry contains
no ultra violet divergences, since the propagator contains the damping factor $e^{-\frac{\Delta}{2}{\bf p}^2}$.
The special shape of the generalized propagator corresponds to the one calculated based on the coherent state
approach to noncommutative geometry \cite{Smailagic:2003rp},\cite{Smailagic:2003yb},\cite{Smailagic:2004yy}.
In \cite{Kober:2010um}, where the additional assumption of noncommuting momenta appears, has been found a
generalization of this propagator even containing no infra red divergences. Accordingly the presented generalized
quantum field theory, which is based on the concept of a nonlocal quantization principle represents at least
under certain consitions a promising candidate for an alternative theory to noncommutative geometry.

\section{Canonical Quantum Gravity with Nonlocal Quantization Principle}

If the nonlocal quantization principle is assumed to be a fundamental principle of nature, then of course
it has also be transferred to the quantum description of the gravitational field. In this section the
corresponding generalization of quantum geometrodynamics is considered. In quantum geometrodynamics
space-time is foliated into a time coordinate and a spacelike hypersurface. Accordingly the metric
field can be expressed in the following way:

\begin{equation}
g_{\mu\nu}=\left(\begin{matrix}N_a N^a-N^2 & N_b\\ N_c & h_{ab}\end{matrix}\right),
\end{equation}
where $h_{ab}$ denotes the three metric on the spacelike hypersurface, $N$ denotes the lapse function
and $N^a$ is the shift vector. The corresponding canonical conjugated momentum $\pi^{ab}$ reads

\begin{equation}
\pi^{ab}=\frac{\partial \mathcal{L}}{\partial \dot h_{ab}}=\frac{\sqrt{h}}{16 \pi G}\left(K^{ab}-K h^{ab}\right),
\label{canonical_conjugated_momentum}
\end{equation}
with $K^{ab}$ describing the extrinsic curvature

\begin{equation}
K_{ab}=\frac{1}{2N}\left(\dot h_{ab}-D_a N_b-D_b N_a\right).
\end{equation}
Since the generalization of the delta function should because of homogenity and isotropy of space just depend on
the distance of the two space points of the field operatopr and its canonical conjugated operator and with respect
to the canonical quantum description of general relativity the distance of two space points depends on the three
metric itself, the commutator depends on the three metric as well. This situation is similar to the generalized
quantization principle which has been developed in \cite{Kober:2011uj}. The distance $D$ between two points,
$x$ and $y$, on the spacelike hypersurface on which the three metric $h_{ab}$ is defined, is given by

\begin{equation}
D(x,y,h)=\int_x^{y} d\lambda\ \sqrt{h_{ab}(\lambda)\ \frac{dx^a}{d\lambda}\frac{dx^b}{d\lambda}}.
\end{equation}
In case of the variables of quantum geometrodynamics, this leads to the following generalized commutation
relation between the components of the three metric and its corresponding canonical conjugated quantity

\begin{equation}
\left[h_{ab}(x),p^{cd}(y)\right]=\frac{i}{2}\left(\delta_a^c \delta_b^d
-\delta_a^d \delta_b^c\right)\alpha^3 \left[D(x,y,h)\right].
\end{equation}
The corresponding nonlocal operators can be expressed as follows:

\begin{eqnarray}
\hat h_{ab}(x)|\Psi[h]\rangle=h_{ab}(x)|\Psi[h]\rangle,\quad
\hat \pi^{ab}(x)|\Psi[h]\rangle=-i\int d^3 y\ \alpha^3\left[D(x,y,h)\right]
\frac{\delta}{\delta h_{ab}(y)}|\Psi[h]\rangle.
\label{nonlocal_operators_gravity}
\end{eqnarray}
And the corresponding quantum constraints are obtained by inserting the operators ($\ref{nonlocal_operators_gravity}$)
to the corresponding classical constraints

\begin{eqnarray}
\hat {\mathcal{H}}(x)|\Psi[h]\rangle&=&\left[16\pi G G_{abcd}[\hat h(x)]\hat \pi^{ab}(x)\hat \pi^{cd}(x)
-\frac{\sqrt{\hat h(x)}}{16\pi G}\left(R[\hat h(x)]-2\Lambda \right)\right]|\Psi[h]\rangle=0,\nonumber\\
&=&\left[16\pi G G_{abcd}[h(x)]\int d^3 y\ d^3 z\ \alpha^3\left[D(x,y,h)\right]\alpha^3\left[D(x,z,h)\right]
\frac{\delta^2}{\delta h_{ab}(y)\delta h_{cd}(z)}
\right.\nonumber\\&&\left.
+\frac{\sqrt{h(x)}}{16\pi G}\left(R[h(x)]-2\Lambda \right)\right]|\Psi[h]\rangle=0,\nonumber\\
\hat {\mathcal{H}}_a(x)|\Psi[h]\rangle&=&-2 D_b \hat \pi_a^b(x)|\Psi[h]\rangle=
-2 D_b h_{ac}(x)\int d^3 z\ \alpha^3 \left[D(x,y,h)\right]\frac{\delta}{\delta h_{bc}(y)}|\Psi[h]\rangle=0,
\end{eqnarray}
where $G_{abcd}$ is defined according to

\begin{equation}
G_{abcd}=\frac{1}{2\sqrt{h}}\left(h_{ac}h_{bd}+h_{ad}h_{bc}-h_{ab}h_{cd}\right).
\end{equation}
If the special case of the assumption of a Gaussian function for the quantization function is considered,
this function is in case of quantum geometrodynamics of the following shape:

\begin{equation}
\alpha^3[D(x,y,h)]=\frac{e^{-\frac{1}{2}\frac{[D(x,y)]^2}{\Delta}}}{\Delta \sqrt{2\pi}},
\end{equation}
where the squared difference between the space-time points has to be replaced by the distance
function $D(x,y,h)$ depending on the three metric $h_{ab}$, too.

\section{Path Integral Quantization}

Of course it has to be possible to formulate the generalized quantization principle, which is given in
($\ref{generalization_commutator}$) by generalization of the usual commutation relations belonging
to canonical quantization, within the path integral formalism as well. To obtain such a formulation in
the framework of path integrals, first the general expression of the transition amplitude between
two states separated by an infinitesimal time interval has to be modified according to the
generalization of the canonical quantization condition. The transition amplitude between a state
$|\Phi(t)\rangle$ at time $t$ and a state $|\Phi^{\prime}(t+dt)\rangle$ at time $t+dt$ is given by

\begin{eqnarray}
\langle \Phi^{\prime}(t+dt)|\Phi(t)\rangle
&=&\int d \Pi\ \langle \Phi^{\prime}|\exp\left[-i\ dt\int d^3 x\ \hat {\mathcal{H}}\left(\hat \Phi(x),\hat \Pi(x)\right)\right]
|\Pi \rangle \langle \Pi|\Phi\rangle\nonumber\\
&=&\int d \Pi\ \langle \Phi^{\prime}|\exp\left[-i\ dt\int d^3 x\ {\mathcal{H}}\left(\Phi(x),\Pi(x)\right)\right]
|\Pi \rangle \langle \Pi|\Phi\rangle.
\label{transition_amplitude_A}
\end{eqnarray}
Usually the inner product between an arbitrary eigenstate of the field operator and an arbitrary eigenstate of the
momentum operator reads $\rangle \langle \Pi|\Phi\rangle=\exp\left[i\int d^3 x\ \Phi(x)\Pi(x)\right]$. In case of
the nonlocal quantization principle the general inner product is given by ($\ref{innerproduct_PhiPi}$). Inserting
of ($\ref{innerproduct_PhiPi}$) to ($\ref{transition_amplitude_A}$) leads to

\begin{eqnarray}
\langle \Phi^{\prime}(t+dt)|\Phi(t)\rangle&=&\int d \Pi\ \exp\left[-i\ dt\int d^3 x\ {\mathcal{H}}\left(\Phi(x),\Pi(x)\right)
+i\int d^3 w\ d^3 z\ {\bar \alpha}^3(w-z)\left(\Phi^{\prime}(w)-\Phi(w))\right)\Pi(z)\right].\nonumber\\
\label{transition_amplitude_B}
\end{eqnarray}
By using the expression of the transition amplitude for an infinitesimal time interval in case of the nonlocal quantization
principle, the corresponding transition amplitude for a non infinitesimal time interval can be determined as usual in
the following way:

\begin{eqnarray}
\langle \Phi(t_b)|\Phi(t_a)\rangle&=&\lim_{N \to \infty}\int d \Phi_1...d \Phi_N
\langle \Phi(t_b)|\Phi_N \rangle...\langle \Phi_1 |\Phi(t_a)\rangle\nonumber\\
&=&\lim_{N \to \infty}\int \prod_{k=1}^N d \Phi_k \prod_{k=0}^N d \Pi_k
\exp\left[i\sum_{k=1}^{N+1}\int d^3 x\left(\int d^3 z\ \bar \alpha^3(x-z)\left(\Phi_k(x)-\Phi_{k-1}(x)\right)\Pi_{k-1}(z)
\right.\right.\nonumber\\&&\left.\left.
\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad
\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad
-\mathcal{H}\left(\Phi_k(x),\Pi_{k-1}(x)\right)\right)\right]\nonumber\\
&=&\int \mathcal{D}\left[\Phi(x)\right]\mathcal{D}\left[\Pi(x)\right]\exp \left[i\int_{t_a}^{t_b} dt \int d^3 x
\left(\int d^3 z\ \bar \alpha^3(x-z)\Pi(x)\dot \Phi(z)-\mathcal{H}\left(\Phi(x),\Pi(x)\right)\right)\right].
\end{eqnarray}
This represents the analogue generalization from the usual quantization principle to the nonlocal
quantizastion principle in the path integral formalism.

\section{Summary and Discussion}

In this paper has been suggested a generalization of the concept of quantization within quantum field theory.
The generalization consists in the assumption that the commutator between the operator describing the
field, which is quantized, and the operator describing its corresponding canonical conjugated momentum
does not vanish, if the two operators refer to different space-time points, as it is presupposed in the
usual setting of quantum field theory. Accordingly the delta function appearing on the right hand side
of the quantization condition usually, has to be replaced by another function, which does not vanish,
if the space-time points the two operators within the commutator refer to are not equal. Because of
the postulate of homogenity and isotropy of space, this function should be assumed just to depend on
the distance between the two space coordinates. According to this nonlocal quantization principle
generalized field operators have to be defined fulfilling this generalized quantization conditions.
Based on the generalized expression of a scalar field operator, the corresponding generalized expression
for the propagator has been calculated, which of course directly depends of the generalization of the
delta function appearing in the generalized quantization condition. If the natural assumption is made
that the quantization function is described by a Gaussian function, then one obtains a special shape
for the propagator with an exponential containg the negative square of the three momentum as damping
factor. Thus the propagator corresponds to the propagator found in the coherent state approach
of noncommutative geometry and contains no ultra violet divegences. This is a very important hint that
this generalization of quantum field theory omitting the postulate of locality with respect to the
basic quantization condition could indeed represent a promising candidate for a generalization of
quantum field theory similar to quantum field theory on noncommutative space-time. Since such a
quantization principle should be interpreted as fundamental property of nature, it has analogously
to be transferred to the quantum description of all field theories including general relativity.
This implies the necessity to generalize the canonical quantum description of general relativity
as well. In this paper has been given a consideration with respect to quantum geometrodynamics, which
has not been elaborated completely yet. But the quantization function has in each case to be assumed
to depend on the gravitational field, since the distance between to space points depends on the
gravitational field. This leads to a similarity to a generalized quantization principle within
canonical quantum gravity \cite{Kober:2011uj} in accordance with a generalized uncertainty in quantum
mechanics, where the right hand side of the quantization condition depends on the operators themselves.
It has also to be possible to incorporate the nonlocal generalization of the quantization principle to
path integrals. This has been done by replacing the usual inner product between an arbitrary eigenstate
of the field operator and an arbitrary eigenstate of its canonical conjugated variable by the
corresponding generalized one. Based on this generalization of path integrals it is in principle
possible to incorporate the nonlocal quantization principle to gravity by using the path integral
formalism. In general it could be very interesting with respect to further research projects, to
formulate a full quantum description of general relativity based on this concept. Perhaps because of
the ultraviolet behaviour of the commutator it can even contribute to circumvent the problem of
nonrenormalizability of the quantrum description of gravity in the framework of usual quantum field
theory. An additional important task would be to determine an upper limit for the parameter $\Delta$
within the Gaussian function, since the generalization of the quantization principle must not get in
conflict with the usual description of quantum field theory at low energies.


\section{Introduction}
It is well known that the $n$-th Bernoulli number $B_{n}$ is defined by the series
\[\frac{x}{{{e}^{x}}-1}=\sum\limits_{n=0}^{\infty }{\frac{{{B}_{n}}}{n!}{{x}^{n}}}.\]
For example, $B_{0}=1, B_{1}=-\frac{1}{2}, B_{2}=\frac{1}{6}$ and $B_{2n+1}=0$ for all $n \ge 1$. There are many fascinating congruences related to some special sums and Bernoulli numbers. For instance, using partial
sum of multiple zeta series, Zhao \cite{Zhao1} proved that for any prime $p \ge 3$,
\begin{equation}\label{zhao1}
\sum\limits_{\begin{smallmatrix}
 i+j+k=p \\
 i,j,k>0
\end{smallmatrix}}^{{}}{\frac{1}{ijk}\equiv -2{{B}_{p-3}} \pmod {p}.}
\end{equation}
By using some combinatorial identities, Ji \cite{Ji} gave a simple proof of this congruence. Zhou and Cai \cite{zhouxia} gave a generalization by establishing a congruence involving arbitrary number of variables. Namely, they showed that
\begin{equation}\label{zhou}
\sum\limits_{\begin{smallmatrix}
 {{l}_{1}}+\cdots +{{l}_{n}}=p, \\
      {{l}_{1}},\cdots ,{{l}_{n}}>0
\end{smallmatrix}}{\frac{1}{{{l}_{1}}{{l}_{2}}\cdots {{l}_{n}}}}\equiv \left\{ \begin{array}{ll}
   -(n-1)!{{B}_{p-n}} \pmod{p} & \textrm{if $2 \nmid n$;} \\
 -\frac{n}{2(n+1)}n!{{B}_{p-n-1}}p  \pmod{{{p}^{2}}} & \textrm{if $2 | n$,}
\end{array} \right.
\end{equation}
where $p \ge 5$ is a prime and $n \le p-2$ is a positive integer.

In 2014, by replacing prime $p$ to prime power, the author and Cai \cite{WangCai} gave a new generalization of (\ref{zhao1}).
Let ${\mathcal{P}_{n}}$ denote the set of positive integers which are prime to  $n$. We proved that for any prime $p \ge 3$,
\begin{equation}\label{Wang}
\sum\limits_{\begin{smallmatrix}
 i+j+k={{p}^{r }} \\
 i,j,k\in {\mathcal{P}_{p}}
\end{smallmatrix}}^{{}}{\frac{1}{ijk}\equiv -2{{p}^{r-1}}{{B}_{p-3}} \pmod{p^r}}.
\end{equation}

In view of (\ref{zhou}), it would be attractive to find some congruences similar to (\ref{Wang}) by increasing the number of variables. For convenience, we define for $1 \le m < n$ that
\[S_{n}^{(m)}({{p}^{r}})=\sum\limits_{\begin{smallmatrix}
 {{l}_{1}}+\cdots +{{l}_{n}}=m{{p}^{r}} \\
 l_{i} < p^{r}, l_{i} \in {\mathcal{P}_{p}}, 1\le i \le n
\end{smallmatrix}}{\frac{1}{{{l}_{1}}{{l}_{2}}\cdots {{l}_{n}}}}.\]
In particular, we also use $S_{n}(p^r)$ to represent $S_{n}^{(1)}(p^r)$ for convention.

Zhao \cite{Zhaoarxiv} found an analogous congruence for $S_{4}(p^r)$. He proved that for any prime $p\ge 5$ and integer $r\ge 2$,
\begin{equation}\label{4variable}
\sum\limits_{\begin{smallmatrix}
 {{l}_{1}}+\cdots +{{l}_{4}}={{p}^{r}} \\
 {{l}_{1}},\cdots ,{{l}_{4}}\in {\mathcal{P}_{p}}
\end{smallmatrix}}{\frac{1}{{{l}_{1}}{{l}_{2}}{{l}_{3}}{{l}_{4}}}}\equiv -\frac{4!}{5}{{p}^{r}}{{B}_{p-5}} \pmod{{p}^{r+1}}.
\end{equation}

Recently, the author \cite{Wang2015} solved the case when there are five variables. We proved that for any prime $p > 5$ and integer $r\ge 2$,
\begin{equation}\label{5variable}
\sum\limits_{\begin{smallmatrix}
 {{l}_{1}}+\cdots +{{l}_{5}}={{p}^{r}} \\
 {{l}_{1}},\cdots ,{{l}_{5}}\in {\mathcal{P}_{p}}
\end{smallmatrix}}{\frac{1}{{{l}_{1}}{{l}_{2}}{{l}_{3}}{{l}_{4}}{{l}_{5}}}}\equiv -\frac{5!}{6}{{p}^{r-1}}{{B}_{p-5}} \pmod{{{p}^{r}}}.
\end{equation}
For other related works, see \cite{CaiShenJia}--\cite{WangJcomb} and \cite{Xia}.

Following their steps, the goal of this paper is to find some analogous result for $S_{6}(p^r)$. Firstly, we establish some congruences analogues to (\ref{zhou}). These congruences are of independent interests themselves  and will also be applied to give the modulo $p^r$ determination of $S_{6}(p^r)$.
\begin{thm}\label{Rn2}
Let $n$ be an even integer and $p>n+2$ be a prime. Then
\begin{displaymath}
\sum\limits_{\begin{smallmatrix}l_1+l_2+\cdots +l_n=2p \\ l_1,l_2,\cdots,l_n \in \mathcal{P}_{p} \end{smallmatrix}}{\frac{1}{l_1l_2\cdots l_n}} \equiv \frac{n!}{2}\sum\limits_{\begin{smallmatrix} a=2 \\ a \,\mathrm{odd}\end{smallmatrix}}^{n-2}\frac{B_{p-a}B_{p-n+a}}{a(n-a)} \pmod{p}.
\end{displaymath}
\end{thm}

\begin{thm}\label{Rn3}
Let $n$ be an even integer and $p>n+2$ be a prime. Then
\begin{displaymath}
\sum\limits_{\begin{smallmatrix}l_1+l_2+\cdots +l_n=3p \\ l_1,l_2,\cdots,l_n \in \mathcal{P}_{p} \end{smallmatrix}}{\frac{1}{l_1l_2\cdots l_n}} \equiv \frac{n!}{6}\sum\limits_{\begin{smallmatrix} a=2 \\ a \,\mathrm{odd}\end{smallmatrix}}^{n-2}{(2n-a+3)\frac{B_{p-a}B_{p-n+a}}{a(n-a)}} \pmod{p}.
\end{displaymath}
\end{thm}

Secondly, utilizing these two congruences, we are able to solve the case of 6 variables. It is worthy mention that by using integer relation detecting tool PSLQ, Zhao \cite{Zhaoarxiv} found that if there exists a constant $c\in \mathbb{Q}$ such that for any $r \ge 2$,
\begin{equation}\label{Zhaoconj}
S_{6}(p^r) \equiv cp^rB_{p-7} \pmod{p^{r+1}},
\end{equation}
then both the numerator and the denominator of $c$ must have at least 60 digits. Our result below shows that the congruence satisfied by $S_{6}(p^r)$ is different from (\ref{Zhaoconj}).
\begin{thm}\label{thm1}
Let $p \ge 11$ be a prime and $r\ge 2$ be an integer. We have
\[\sum\limits_{\begin{smallmatrix}
 {{l}_{1}}+\cdots +{{l}_{6}}={{p}^{r}} \\
 {{l}_{1}},\cdots ,{{l}_{6}}\in {\mathcal{P}_{p}}
\end{smallmatrix}}{\frac{1}{l_1l_2l_3l_4l_5l_6}} \equiv -\frac{5!}{18}{{p}^{r-1}}{{B}_{p-3}^{2}} \pmod{{{p}^{r}}}.\]
\end{thm}
As a by-product, we have
\begin{thm}\label{thm2}
Let $p \ge 11$ be a prime and $n$ be a positive integer. Suppose  $p^{r}|n$ but $p^{r+1} \nmid n$ for some positive integer $r$.\\
(i) If $r=1$, then
\begin{displaymath}
\sum\limits_{\begin{smallmatrix}
 {{l}_{1}}+\cdots +{{l}_{6}}=n \\
 {{l}_{1}},\cdots ,{{l}_{6}}\in {\mathcal{P}_{p}}
\end{smallmatrix}}{\frac{1}{l_1l_2l_3l_4l_5l_6}}\equiv \frac{5!}{18} \cdot \Big(\big(\frac{n}{p}\big)^{5}-\big(\frac{n}{p}\big)^3\Big){{B}_{p-3}^{2}} \pmod{p}.
\end{displaymath}
(ii) If $r\ge 2$, then
\begin{displaymath}
\sum\limits_{\begin{smallmatrix}
 {{l}_{1}}+\cdots +{{l}_{6}}=n \\
 {{l}_{1}},\cdots ,{{l}_{6}}\in {\mathcal{P}_{p}}
\end{smallmatrix}}{\frac{1}{l_1l_2l_3l_4l_5l_6}}\equiv -\frac{5!}{18} \cdot \frac{n}{p}{{B}_{p-3}^{2}} \pmod{{{p}^{r}}}.
\end{displaymath}
\end{thm}
In particular, if $n=p^r$ ($r\ge 2$), then part (ii) of Theorem \ref{thm2} becomes Theorem \ref{thm1}.

The paper is organized as follows. In Sec. 2 we collect some ingredients which will be used to prove the theorems above. In Sec. 3 we will prove Theorems \ref{Rn2} and \ref{Rn3} by using some properties of the $p$-restricted multiple harmonic sums defined by
\begin{displaymath}
H_{N}(\alpha_1,\alpha_2,\cdots, \alpha_n) = \sum\limits_{\begin{smallmatrix}0<k_1<k_2<\cdots <k_n<N \\k_1,k_2,\cdots ,k_n \in \mathcal{P}_{p}\end{smallmatrix}}{\frac{1}{k_{1}^{\alpha_1}k_2^{\alpha_2}\cdots k_n^{\alpha_{n}}}}.
\end{displaymath}
From these two theorems, we can determine $S_{6}^{(k)}(p)$ for $1 \le k \le 5$ modulo $p$. Then in Sec. 4 we establish the relation between $S_{6}^{(k)}(p^{r+1})$ and $S_{6}^{(k)}(p^r)$, from which Theorem \ref{thm1} follows by induction on $r$. Finally, in Sec. 5 we give a rough predication of the general congruences about $R_{n}(p)$ and $S_{n}(p^r)$ for any positive integer $n$. We also point out the difficulties to find the general congruences using the current method.

\section{Preliminaries}
In this section, we introduce some notations and lemmas which will be fundamental to our proofs.

We define for $m\ge 1$ that
\begin{displaymath}
R_n^{(m)}({p^r}) = \sum\limits_{\begin{smallmatrix}{l_1} + {l_2} + \cdots + {l_n} = m{p^r}\\
{l_i} \in {\mathcal{P}_p}, 1 \le i\le n\end{smallmatrix}} {\frac{1}{{{l_1}{l_2}\cdots {l_n}}}}.
\end{displaymath}
The relation between $R_{n}^{(m)}(p)$ and $S_{n}^{(m)}(p)$ is given by the following lemma.
\begin{lem}\label{RSrelation}
For any integers $m$ and $n$, we have
\[S_{n}^{(m)}(p) \equiv \sum\limits_{k=0}^{m-1}{n \choose k}(-1)^kR_{n}^{(m-k)}(p) \pmod{p}.\]
\end{lem}
\begin{proof}
By Inclusion-Exclusion Principle, we have
\begin{displaymath}
\begin{split}
S_{n}^{(m)}(p)&=\sum\limits_{\begin{smallmatrix}l_1+l_2+\cdots +l_n=mp \\l_1,l_2,\cdots, l_n\in \mathcal{P}_{p} \end{smallmatrix}}{\frac{1}{l_1l_2\cdots l_{n}}}+\sum\limits_{k=1}^{m-1}(-1)^{k}\sum\limits_{1\le a_1<\cdots <a_{k}\le n}\sum\limits_{\begin{smallmatrix}l_1+l_2+\cdots +l_n=mp \\l_1,l_2,\cdots, l_n\in \mathcal{P}_{p} \\l_{a_1},\cdots ,l_{a_{k}}>p \end{smallmatrix}}\frac{1}{l_1l_2\cdots l_{n}}\\
&= \sum\limits_{k=0}^{m-1}(-1)^{k}{n \choose k}\sum\limits_{\begin{smallmatrix}l_1+l_2+\cdots +l_n=(m-k)p \\l_1,l_2,\cdots, l_n\in \mathcal{P}_{p}  \end{smallmatrix}}\frac{1}{(l_1+p)(l_2+p)\cdots (l_{k}+p)l_{k+1}\cdots l_{n}}\\
&\equiv \sum\limits_{k=0}^{m-1}{n \choose k}(-1)^k\sum\limits_{\begin{smallmatrix}l_1+l_2+\cdots +l_n=(m-k)p \\l_1,l_2,\cdots, l_n\in \mathcal{P}_{p}  \end{smallmatrix}}\frac{1}{l_1l_2\cdots l_{n}} \pmod{p}.
\end{split}
\end{displaymath}
This completes our proof.
\end{proof}

We also define
\begin{displaymath}
U_{N}(\alpha_1,\alpha_2,\cdots,\alpha_n)=\sum\limits_{\begin{smallmatrix}0<l_1,\cdots, l_n <N \\ l_i \ne l_j, \forall i \ne j, l_i \in \mathcal{P}_{p}\end{smallmatrix}}\frac{1}{l_1^{\alpha_1}\cdots l_{n}^{\alpha_n}}.
\end{displaymath}
There are some relations between $H_{N}$ and $U_{N}$. For example, let $S_{n}$ denote the symmetric group of $\{1,2, \cdots, n\}$. It is easy to see that
\begin{equation}\label{HU}
\sum\limits_{\sigma \in S_{n}}{H_{N}(\alpha_{\sigma(1)}, \alpha_{\sigma(2)}, \cdots, \alpha_{\sigma(n)})} = U_{N}(\alpha_1,\alpha_2,\cdots,\alpha_n).
\end{equation}

\begin{lem}\label{U}
Let $\alpha_1, \cdots, \alpha_n$ be positive integers, $r=\alpha_1+\cdots+ \alpha_n \le p-3$, where $p$ is a prime. For any positive integer $b$, we have
\begin{displaymath}
U_{bp}(\alpha_1, \cdots, \alpha_n) \equiv \left\{ \begin{array}{ll}
   (-1)^{n}(n-1)!\frac{b^2r(r+1)}{2(r+2)}B_{p-r-2}p^2  \pmod{p^3} & \textrm{if $2 \nmid r$;} \\
 (-1)^{n-1}(n-1)!\frac{br}{r+1}B_{p-r-1}p  \pmod{{{p}^{2}}} & \textrm{if $2 | r$.}
\end{array} \right.
\end{displaymath}
\end{lem}
\begin{proof}
When $b=1$, this has already been proved, see \cite[Lemma 3]{zhouxia}. In particular, if $b=n=1$, we obtain
\begin{equation}\label{Uspecial}
\sum\limits_{1 \le l < p}{\frac{1}{l^{\alpha}}}=U_{p}(\alpha) \equiv   \left\{ \begin{array}{ll}
   -\frac{\alpha(\alpha+1)}{2(\alpha+2)}B_{p-\alpha-2}p^2  \pmod{p^3} & \textrm{if $2 \nmid \alpha$;} \\
 \frac{\alpha}{\alpha+1}B_{p-\alpha-1}p  \pmod{{{p}^{2}}} & \textrm{if $2 | \alpha$.}
\end{array} \right.
\end{equation}
For all $k\ge 1$, we have
\begin{displaymath}
\begin{split}
\sum\limits_{kp<l<(k+1)p}{\frac{1}{l^{\alpha}}}&=\sum\limits_{l=1}^{p-1}{\frac{1}{(l+kp)^\alpha}} \\
&\equiv \sum\limits_{l=1}^{p-1}{\Bigg(1-\frac{\alpha kp}{l}+\frac{\alpha(\alpha+1)}{2l^2}k^2p^2\Bigg)\frac{1}{l^\alpha}}\\
&\equiv U_{p}(\alpha)-\alpha kpU_{p}(\alpha+1) +\frac{\alpha(\alpha+1)}{2l^2}k^2p^2U_{p}(\alpha+2) \pmod{p^3}.
\end{split}
\end{displaymath}
By (\ref{Uspecial}) we see that
\begin{displaymath}
\sum\limits_{kp < l< (k+1)p}{\frac{1}{l^{\alpha}}} \equiv   \left\{ \begin{array}{ll}
   -\frac{\alpha(\alpha+1)}{(\alpha+2)}\Big(\frac{1}{2}+k\Big)B_{p-\alpha-2}p^2  \pmod{p^3} & \textrm{if $2 \nmid \alpha$;} \\
 \frac{\alpha}{\alpha+1}B_{p-\alpha-1}p  \pmod{{{p}^{2}}} & \textrm{if $2 | \alpha$.}
\end{array} \right.
\end{displaymath}
Hence we have
\begin{displaymath}
\sum\limits_{1 \le l < bp, p\nmid l}{\frac{1}{l^{\alpha}}} \equiv   \left\{ \begin{array}{ll}
   -\frac{b^{2}\alpha(\alpha+1)}{2(\alpha+2)}B_{p-\alpha-2}p^2  \pmod{p^3} & \textrm{if $2 \nmid \alpha$;} \\
 \frac{b\alpha}{\alpha+1}B_{p-\alpha-1}p  \pmod{{{p}^{2}}} & \textrm{if $2 | \alpha$.}
\end{array} \right.
\end{displaymath}
This proves the lemma for $n=1$. Now assume the lemma is true when the number of variables is less than $n$. We have
\begin{displaymath}
\begin{split}
U_{bp}(\alpha_1,\cdots, \alpha_n)&=\sum\limits_{\begin{smallmatrix} 1\le l_1, \cdots, l_{n-1} <bp \\ l_i \ne l_j, l_i \in \mathcal{P}_{p}\end{smallmatrix}}{\frac{1}{l_1^{\alpha_1}\cdots l_{n-1}^{\alpha_{n-1}}}\Bigg(\sum\limits_{1\le l_{n}< bp, l_n \in \mathcal{P}_{p}}{\frac{1}{l_n^{\alpha_n}}}-\sum\limits_{i=1}^{n-1}{\frac{1}{l_i^{\alpha_n}}}\Bigg)} \\
&\equiv U_{bp}(\alpha_1, \cdots, \alpha_{n-1})\Bigg(\sum\limits_{1\le l_{n}< bp, l_n \in \mathcal{P}_{p}}{\frac{1}{l_n^{\alpha_n}}}\Bigg) \\
&\quad -\sum\limits_{i=1}^{n-1}{U_{bp}(\alpha_1, \cdots, \alpha_{i-1}, \alpha_{i}+\alpha_{n}, \alpha_{i+1}, \cdots, \alpha_{n-1})}.
\end{split}
\end{displaymath}
From the assumption above, we have
\begin{displaymath}
 U_{bp}(\alpha_1, \cdots, \alpha_{n-1})\Bigg(\sum\limits_{1\le l_{n}< bp, l_n \in \mathcal{P}_{p}}{\frac{1}{l_n^{\alpha_n}}}\Bigg) \equiv \left\{\begin{array}{ll}
 0 \pmod{p^3}, &\textrm{if $2 \nmid r$;} \\
0 \pmod{p^2}, &\textrm{if $2|r$.}
\end{array} \right.
\end{displaymath}
If $r$ is odd, then
\begin{displaymath}
\begin{split}
U_{bp}(\alpha_1,\cdots, \alpha_n)& \equiv -(n-1)U_{bp}(\alpha_1, \cdots, \alpha_{n-2}, \alpha_{n-1}+\alpha_{n})\\
&\equiv -(n-1)(-1)^{n-1}(n-2)!\frac{b^2r(r+1)}{2(r+2)}p^2B_{p-r-2} \\
&\equiv (-1)^n(n-1)!\frac{b^2r(r+1)}{2(r+2)}p^2B_{p-r-2} \pmod{p^3}.
\end{split}
\end{displaymath}
If $r$ is even, similarly we can deduce
\begin{displaymath}
U_{bp}(\alpha_1,\cdots, \alpha_n) \equiv (-1)^{n-1}(n-1)!\frac{br}{r+1}pB_{p-r-1} \pmod{p^2}.
\end{displaymath}
The proof of Lemma \ref{U} is complete by induction on $n$.
\end{proof}

For fixed prime $p$, let $C^{(m)}_{a}$ ($1 \le a, m \le 5$) denote the number of integer solutions $(x_1,\dots,x_6)$ of the equation
$$x_1 + \dots + x_6 =mp-a, \quad 0 \le x_i < p, \quad \forall i=1,\cdots ,6.$$
We need the following facts.
\begin{lem}\label{Cam}
For any prime $p$, we have \\
$(i)$ $C_{a}^{(m)} \equiv 0 \pmod{p}$;\\
$(ii)$ $C_{1}^{(1)}+C_{5}^{(1)} \equiv \frac{2}{5}p \pmod{p^2}$, $C_{2}^{(1)} +C_{4}^{(1)}\equiv -\frac{1}{10}p \pmod{p^2}$, $C_{3}^{(1)} \equiv \frac{1}{30}p \pmod{p^2}$; \\
$(iii)$ $C_{1}^{(2)}+C_{5}^{(2)} \equiv -\frac{8}{5}p \pmod{p^2}$, $C_{2}^{(2)}+C_{4}^{(2)} \equiv \frac{2}{5}p \pmod{p^2}$, $C_{3}^{(2)}\equiv -\frac{2}{15}p \pmod{p^2};$ \\
$(iv)$ $C_{1}^{(3)}+C_{5}^{(3)} \equiv \frac{12}{5}p \pmod{p^2}$, $C_{2}^{(3)}+C_{4}^{(3)} \equiv -\frac{3}{5}p \pmod{p^2}$, $C_{3}^{(3)}\equiv \frac{1}{5}p \pmod{p^2}$.
\end{lem}
\begin{proof}
By Inclusion-Exclusion Principle, we have
\begin{eqnarray*}
{{C}_{a}^{(m)}} &=&  \sum\limits_{k=1}^{m}{(-1)^{k-1}\#\Big\{(x_{1},\cdots,x_{6})\in \mathbb{N}^{6}\mid{{{x}_{1}}+\cdots +{{x}_{6}}=mp-a,\exists x_{i_{1}},\cdots,x_{i_{k-1}}\ge p}\Big\}} \\
 & =& \sum\limits_{k=1}^{m}{(-1)^{k-1}\#\Big\{(x_{1},\cdots,x_{6})\in \mathbb{N}^{6}\left|{{{x}_{1}}+\cdots +{{x}_{6}}=(m-k+1)p-a}\right.\Big\}}\binom{6}{k-1} \\
& =&\sum\limits_{k=1}^{m}{(-1)^{k-1}\binom{(m-k+1)p-a+5}{5}}\binom{6}{k-1}\\
\end{eqnarray*}

Let $X_{a}(u)=\binom{up-a+5}{5}$. It follows that
\begin{equation}\label{Camgen}
C_{a}^{(m)}=\sum\limits_{k=1}^{m}{(-1)^{k-1}X_{a}(m-k+1)\binom{6}{k-1}}
\end{equation}
In particular, we have
\begin{equation}\label{Camexpress}
\begin{split}
C_{a}^{(1)}&=X_{a}(1), \quad C_{a}^{(2)}=X_{a}(2)-6X_{a}(1), \\
C_{a}^{(3)}&=X_{a}(3)-6X_{a}(2)+15X_{a}(1).
\end{split}
\end{equation}

Since $1\le a \le 5$, we have
\begin{equation}\label{Xau}
\begin{split}
X_{a}(u) &=\frac{(up-a+5)(up-a+4)(up-a+3)(up-a+2)(up-a+1)}{5!} \\
& \equiv (-1)^{(a-1)}\frac{(a-1)!(5-a)!}{5!}up \pmod{p^2}.
\end{split}
\end{equation}
Hence from (\ref{Camgen}) we see that $C_{a}^{(m)} \equiv 0$ (mod $p$) and (i) is proved.

Furthermore, from (\ref{Xau}) we deduce that
\begin{equation}\label{Xauadd}
\begin{split}
X_{1}(u)+X_{5}(u)&\equiv \frac{2}{5}up \pmod{p^2}, \quad X_3(u) \equiv \frac{1}{30}up \pmod{p^2} \quad \textrm{and} \\
X_2(u)+X_4(u)&\equiv -\frac{1}{10}up \pmod{p^2}.
\end{split}
\end{equation}
From (\ref{Camexpress}) and (\ref{Xauadd}), (ii)--(iv) follows by simple calculations.
\end{proof}

\section{Proofs of Theorems \ref{Rn2} and \ref{Rn3}}
Let $u_{i}=l_1+\cdots+l_{i}$ for $i=1,2,\cdots n-1$. We have
\begin{eqnarray}\label{start}
R_{n}^{(m)}(p)&=&\frac{1}{mp}\sum\limits_{\begin{smallmatrix}l_1+l_2+\cdots+l_n=mp \\ l_{1},l_{2},\cdots, l_{n}\in \mathcal{P}_{p} \end{smallmatrix}}{\frac{l_1+l_2+\cdots+l_n}{l_1l_2\cdots l_n}}  \nonumber\\
&=&\frac{n}{mp}\sum\limits_{\begin{smallmatrix} u_{n-1}<mp \\ l_1, \cdots,l_{n-1},u_{n-1} \in \mathcal{P}_{p}\end{smallmatrix}}{\frac{1}{l_1l_2\cdots l_{n-1} }} \nonumber\\
&=&\frac{n}{mp}\sum\limits_{\begin{smallmatrix} u_{n-1}<mp \\ l_1, \cdots,l_{n-1},u_{n-1} \in \mathcal{P}_{p}\end{smallmatrix}}{\frac{l_1+l_2+\cdots+l_{n-1}} {l_1l_2\cdots l_{n-1}u_{n-1} }} \nonumber\\
&=&\frac{n(n-1)}{mp}\sum\limits_{\begin{smallmatrix} u_{n-2}<u_{n-1}<mp \\ l_1, \cdots,l_{n-2},u_{n-1}-u_{n-2},u_{n-1} \in \mathcal{P}_{p}\end{smallmatrix}}{\frac{1} {l_1l_2\cdots l_{n-2}u_{n-1} }} =\cdots \nonumber\\
&=&\frac{n!}{mp}\sum\limits_{\begin{smallmatrix} u_1<\cdots <u_{n-1}<mp \\ u_1, u_2-u_1,\cdots,u_{n-1}-u_{n-2},u_{n-1} \in \mathcal{P}_{p}\end{smallmatrix}}{\frac{1} {u_1u_2\cdots u_{n-1} }}.
\end{eqnarray}
\begin{proof}[Proof of Theorem \ref{Rn2}]
Let $m=2$ in (\ref{start}). We have
\begin{eqnarray}\label{sumformula}
R_{n}^{(2)}(p)&=&\frac{n!}{2p}\Bigg(\sum\limits_{\begin{smallmatrix} u_1<\cdots<u_{n-1}<2p\\ u_1,\cdots,u_{n-1}\in \mathcal{P}_{p} \end{smallmatrix}}{\frac{1}{u_1u_2\cdots u_{n-1}}}-\sum\limits_{a=1}^{n-2}\sum\limits_{\begin{smallmatrix} u_1<\cdots<u_{n-1}<2p\\ u_{a+1}-u_{a}=p, u_1,\cdots,u_{n-1}\in \mathcal{P}_{p} \end{smallmatrix}}{\frac{1}{u_1u_2\cdots u_{n-1} }} \nonumber \\
&&+\sum\limits_{a=2}^{n-2}\sum\limits_{\begin{smallmatrix} u_1<\cdots<u_{n-1}<2p\\ u_{a}=p,u_{j}\in \mathcal{P}_{p}, \forall j\ne a \end{smallmatrix}}{\frac{1}{u_1u_2\cdots u_{n-1}}} \Bigg)\\
&=&T_1-T_2+T_3, \nonumber
\end{eqnarray}
here we denote the three terms in the right hand side of (\ref{sumformula}) by $T_{1}$, $T_{2}$ and $T_{3}$ respectively. We will deal with these terms one by one.

For the first term, since $n-1$ is odd, by (\ref{HU}) and Lemma \ref{U} we deduce that
\begin{equation}\label{T1}
T_1=\frac{n!}{2p}H_{2p}(\{1\}^{n-1})=\frac{n}{2p}U_{2p}(\{1\}^{n-1})\equiv 0 \pmod{p}.
\end{equation}

For the second term,  according to (\ref{HU}) and Lemma \ref{U}, we have
\begin{eqnarray*}
&&\sum\limits_{a=1}^{n-2}\sum\limits_{\begin{smallmatrix}u_1<u_1<\cdots <u_{n-1}<2p\\ u_{a+1}-u_{a}=p \end{smallmatrix}}{\frac{1}{u_1u_2\cdots u_{n-1}} }  \quad (\textrm{replace $u_{k+1}$ by $u_{k}+p$ for $k\ge a$} )\\
&=& \sum\limits_{a=1}^{n-2}\sum\limits_{u_1<\cdots <u_{n-2}<p}{\frac{1}{u_1\cdots u_a(u_a+p)(u_{a+1}+p)\cdots (u_{n-2}+p) }}\\
&\equiv &\sum\limits_{a=1}^{n-2}\sum\limits_{u_1<\cdots <u_{n-2}<p}{\frac{1}{u_1\cdots u_{a}^2u_{a+1}\cdots u_{n-2}}\Big(1-\frac{p}{u_a}-\frac{p}{u_{a+1}}-\cdots -\frac{p}{u_{n-2}}\Big)} \\
&\equiv & \sum\limits_{a=1}^{n-2}H_{p}(\{1\}^{a-1},2,\{1\}^{n-2-a})-p\sum\limits_{a=1}^{n-2}H_{p}(\{1\}^{a-1},3,\{1\}^{n-2-a})\\
&& -p\sum\limits_{1\le a <b\le n-2}{H_{p}(\{1\}^{a-1},2,\{1\}^{b-a-1},2,\{1\}^{n-2-b})}\\
&\equiv& \frac{n-2}{(n-2)!}U_{p}(2,\{1\}^{n-3})-p\cdot \frac{n-2}{(n-2)!}U_{p}(3,\{1\}^{n-3})-p\frac{{n-2 \choose 2}}{(n-2)!}U_{p}(2,2,\{1\}^{n-4}) \\
&\equiv &0 \pmod{p^2}.
\end{eqnarray*}
Therefore, we have
\begin{equation}\label{T2}
T_{2}\equiv 0 \pmod{p}.
\end{equation}

Now we calculate the third term. For $2\le a \le n-2$, from (\ref{HU}) and Lemma \ref{U} we obtain that
\begin{displaymath}
\begin{split}
&\sum\limits_{\begin{smallmatrix} u_1<\cdots<u_{n-1}<2p\\ u_{a}=p,u_1,\cdots,u_{n-1}\in \mathcal{P}_{p} \end{smallmatrix}}{\frac{1}{u_1u_2\cdots u_{n-1}}} \quad (\textrm{replace $u_{k}$ by $x_{k-a}+p$ for $k>a$})\\
&=\frac{1}{p}\Bigg(\sum\limits_{u_1<\cdots<u_{a-1}<p}{\frac{1}{u_1u_2\cdots u_{a-1}} } \Bigg)\Bigg(\sum\limits_{x_{1}<\cdots <x_{n-1-a}<p}{\frac{1}{(x_{1}+p)\cdots (x_{n-1-a}+p)}} \Bigg)\\
&\equiv \frac{1}{p}H_{p}(\{1\}^{a-1}) \sum\limits_{x_1<\cdots <x_{n-1-a}<p}{\frac{1}{x_1\cdots x_{n-1-a}}\Big(1-\frac{p}{x_1}-\cdots -\frac{p}{x_{n-1-a}}\Big) } \\
&\equiv \frac{1}{p}H_{p}(\{1\}^{a-1})\Big(H_{p}(\{1\}^{n-1-a})-p\cdot \frac{n-1-a}{(n-1-a)!}U_{p}(2,\{1\}^{n-2-a}) \Big)\\
&\equiv \frac{1}{p}\cdot \frac{U_{p}(\{1\}^{a-1}) }{(a-1)! }\cdot \frac{U_{p}(\{1\}^{n-1-a})}{(n-1-a)!} \\
&\equiv \left\{\begin{array}{ll}
\frac{B_{p-a}B_{p-n+a}}{a(n-a)}p & \textrm{if $a$ is odd}\\
0 & \textrm{if $a$ is even}
\end{array} \right. \pmod{p^2}.
\end{split}
\end{displaymath}
Therefore, we have
\begin{equation}\label{T3}
T_3\equiv \frac{n!}{2p}\sum\limits_{\begin{smallmatrix}a=2\\ a \,\, \textrm{odd} \end{smallmatrix}}^{n-2}{\frac{B_{p-a}B_{p-n+a}}{a(n-a)}p }\\
\equiv \frac{n!}{2}\sum\limits_{\begin{smallmatrix}a=2\\ a \,\,\textrm{odd} \end{smallmatrix}}^{n-2}{\frac{B_{p-a}B_{p-n+a}}{a(n-a)} } \pmod{p}.
\end{equation}
Substituting (\ref{T1})--(\ref{T3}) into (\ref{sumformula}), we complete the proof of Theorem \ref{Rn2}.
\end{proof}
Now we are able to determine $S_{n}^{(2)}(p)$ modulo $p$.
\begin{cor}\label{Sn2}
Let $n$ be an even integer and $p>n+2$ be a prime. We have
\begin{displaymath}
\sum\limits_{\begin{smallmatrix}l_1+l_2+\cdots +l_n=2p\\ l_1,l_2,\cdots ,l_n<p \end{smallmatrix}}\frac{1}{l_1l_2\cdots l_n} \equiv \frac{n!}{2}\sum\limits_{\begin{smallmatrix} a=2 \\ a \,\mathrm{odd}\end{smallmatrix}}^{n-2}\frac{B_{p-a}B_{p-n+a}}{a(n-a)} \pmod{p}.
\end{displaymath}
\end{cor}
\begin{proof}
From (\ref{zhou}) we know $R_{n}^{(1)}(p)\equiv S_{n}^{(1)}(p)\equiv 0$ (mod $p$). By Lemma \ref{RSrelation}, we obtain
\[S_{n}^{(2)}(p)\equiv R_{n}^{(2)}(p)-nR_{n}^{(1)}(p) \equiv R_{n}^{(2)}(p) \pmod{p}.\]
The corollary then follows from Theorem \ref{Rn2}.
\end{proof}

\begin{proof}[Proof of Theorem \ref{Rn3}]
Let $m=3$ in (\ref{start}). We have
\begin{equation}\label{Rn3start}
\begin{split}
R_{n}^{(3)}(p)& =\frac{n!}{3p}\Bigg(\sum\limits_{\begin{smallmatrix} u_1<\cdots <u_{n-1}<3p \\ u_1, \cdots, u_{n-1} \in \mathcal{P}_{p} \\ u_2-u_1, \cdots, u_{n-1}-u_{n-2}\in \mathcal{P}_{p}\end{smallmatrix}}+\sum\limits_{i=2}^{n-4}\sum\limits_{j=i+2}^{n-2}\sum\limits_{\begin{smallmatrix}1\le u_1<\cdots <u_{n-1}<3p \\ u_{i}=p,u_{j}=2p, u_{k}\in \mathcal{P}_{p}, \forall k\ne i,j\\u_2-u_1, \cdots, u_{n-1}-u_{n-2} \in \mathcal{P}_{p} \end{smallmatrix}}\\
&\quad +\sum\limits_{j=2}^{n-2}\sum\limits_{\begin{smallmatrix}u_1<\cdots <u_{n-1}<3p, u_{j}=p \\ u_2-u_2, \cdots, u_{n-1}-u_{n-2}\in \mathcal{P}_{p} \\u_{k}\in \mathcal{P}_{p}\,(k \ne j)\end{smallmatrix}}+\sum\limits_{j=2}^{n-2}\sum\limits_{\begin{smallmatrix}u_1<\cdots <u_{n-2}<3p, u_{j}=2p \\ u_2-u_2, \cdots, u_{n-1}-u_{n-2}\in \mathcal{P}_{p} \\u_{k}\in \mathcal{P}_{p}\,(k \ne j)\end{smallmatrix}}{\frac{1}{u_1u_2\cdots u_{n-2}}} \Bigg)\\
& =\frac{n!}{3p}\Big(T_1+T_2+T_3+T_4\Big),
\end{split}
\end{equation}
here we denote the four sums in the bracket by $T_{1}$, $T_{2}$, $T_{3}$ and $T_{4}$, respectively.

For the first sum in (\ref{Rn3start}), we have
\begin{equation}\label{T1start}
\begin{split}
T_{1}=\sum\limits_{\begin{smallmatrix}u_{1}<\cdots <u_{n-1}<3p\\ u_1,\cdots, u_{n-1}\in \mathcal{P}_{p}  \end{smallmatrix}}-\sum\limits_{i=1}^{n-2}\sum\limits_{a=1}^{2}\sum\limits_{\begin{smallmatrix} u_1<\cdots <u_{n-1}<3p \\u_{i+1}-u_{i}=ap \\u_{1}, \cdots, u_{n-1} \in \mathcal{P}_{p} \end{smallmatrix}}+\sum\limits_{1\le i <j \le n-2}\sum\limits_{\begin{smallmatrix} u_1<\cdots <u_{n-1}<3p \\u_{i+1}-u_{i}=u_{j+1}-u_{j}=p \\ u_{1},u_{2},\cdots,u_{n-1}\in \mathcal{P}_{p}\end{smallmatrix}}\frac{1}{u_1u_2\cdots u_{n-1}}
\end{split}
\end{equation}
For $a=1$ or 2, by (\ref{HU}) and Lemma \ref{U} we have
\begin{eqnarray}\label{aterm}
&&\sum\limits_{i=1}^{n-2}\sum\limits_{\begin{smallmatrix}u_1<\cdots <u_{n-1}<3p \\u_{i+1}-u_{i}=ap \\ u_{1},u_{2},\cdots,u_{n-1}\in \mathcal{P}_{p}\end{smallmatrix}}\frac{1}{u_1u_2\cdots u_{n-1}} \quad (\textrm{replace $u_{k}$ by $u_{k-1}+ap$ for $k>i$}) \nonumber\\
&=&\sum\limits_{i=1}^{n-2}\sum\limits_{\begin{smallmatrix}u_1<\cdots  <u_{n-2}<(3-a)p \\u_1,\cdots ,u_{n-2} \in \mathcal{P}_{p}  \end{smallmatrix}}{\frac{1}{u_1\cdots u_{i}(u_{i}+ap)(u_{i+1}+ap)\cdots (u_{n-2}+ap)}}  \nonumber\\
&\equiv& \sum\limits_{i=1}^{n-2}\sum\limits_{\begin{smallmatrix}u_1<\cdots <u_{n-2}<(3-a)p\\ u_1,\cdots, u_{n-2}\in \mathcal{P}_{p}\end{smallmatrix}}\frac{1}{u_1\cdots u_{i-1}u_{i}^{2}u_{i+1}\cdots u_{n-2}}\Big(1-\frac{ap}{u_{i}}-\cdots -\frac{ap}{u_{n-2}}\Big)  \nonumber\\
&\equiv& \frac{n-2}{(n-2)!}U_{(3-a)p}(2,\{1\}^{n-3})-ap\frac{(n-2)U_{(3-a)p}(3,\{1\}^{n-3})}{(n-2)!}  \nonumber\\
&&-ap\frac{{n-2 \choose 2}}{(n-2)!}U_{(3-a)p}(2,2,\{1\}^{n-4})  \nonumber \\
&\equiv& 0 \pmod{p^2}.
\end{eqnarray}

According to $j=i+1$ or not, we split the third sum in (\ref{T1start}) into two parts. For the case $j\ne i+1$, replacing $u_{k}$ by $u_{k-1}+p$ for $i<k\le j$ or $u_{k-2}+2p$ for $k>j$, applying (\ref{HU}) and Lemma \ref{U}, we deduce that
\begin{eqnarray}\label{minisum2}
&&\sum\limits_{\begin{smallmatrix}1 \le i <j \le n-2 \\j \ne i+1  \end{smallmatrix}}\sum\limits_{\begin{smallmatrix}u_1<\cdots <u_{n-1}<3p\\ u_{i+1}-u_{i}=u_{j+1}-u_{j}=p \end{smallmatrix}}\frac{1}{u_1u_2\cdots u_{n-1}}\nonumber \\
&=&\sum\limits_{\begin{smallmatrix}1 \le i <j \le n-2 \\j\ne i+1 \end{smallmatrix} }\sum\limits_{u_1<\cdots <u_{n-3}<p }\frac{1}{u_1\cdots u_{i}(u_{i}+p)\cdots (u_{j-1}+p)(u_{j-1}+2p)\cdots (u_{n-3}+2p)} \nonumber\\
&& (\text{\rm{replace $j$ by $j+1$}}) \nonumber\\
&=&\sum\limits_{1 \le i <j \le n-3 }\sum\limits_{u_1<\cdots <u_{n-3}<p }\frac{1}{u_1\cdots u_{i}(u_{i}+p)\cdots (u_{j}+p)(u_{j}+2p)\cdots (u_{n-3}+2p)}  \nonumber\\
&\equiv& \sum\limits_{1 \le i <j \le n-3 }\sum\limits_{u_1<\cdots <u_{n-3}<p }{\frac{1}{u_1\cdots u_{i-1}u_{i}^{2}u_{i+1}\cdots u_{j-1}u_{j}^{2}u_{j+1}\cdots u_{n-3}}} \nonumber\\
&&{\cdot \Big(1-\frac{p}{u_i}-\cdots-\frac{p}{u_{j}}-\frac{2p}{u_j}-\cdots -\frac{2p}{u_{n-3}} \Big) }  \nonumber\\
&\equiv& \frac{{n-3 \choose 2}}{(n-3)!}U_{p}(2,2,\{1\}^{n-5})-p\sum\limits_{1 \le i <j \le n-3 }{H_{p}(\{1\}^{i-1},3,\{1\}^{j-i-1},2,\{1\}^{n-3-j})}\nonumber\\
&&-p\sum\limits_{1\le i <k<j\le n-3}H_{p}(\{1\}^{i-1},2,\{1\}^{k-i-1},2,\{1\}^{j-k-1},2,\{1\}^{n-j-3}) \nonumber \\
&& -3p\sum\limits_{1\le i <j\le n-3 }H_{p}(\{1\}^{i-1},2,\{1\}^{j-1-i},3,\{1\}^{n-3-j}) \nonumber \\
&& -2p\sum\limits_{1\le i <j<k \le n-3}H_{p}(\{1\}^{i-1},2,\{1\}^{j-1-i},2,\{1\}^{k-j-1},2,\{1\}^{n-3-k})  \nonumber\\
&\equiv& -p\sum\limits_{1\le i<j\le n-3}{ \Big(H_{p}(\{1\}^{i-1},3,\{1\}^{j-i-1},2,\{1\}^{n-3-j})}  \nonumber\\
&&{+3H_{p}(\{1\}^{i-1},2,\{1\}^{j-i-1},3,\{1\}^{n-3-j} )   \Big)  }\pmod{p^2},
\end{eqnarray}
where the last congruence equality follows from the facts that
\begin{displaymath}
\begin{split}
&\quad \sum\limits_{1\le i <j<k \le n-3}H_{p}(\{1\}^{i-1},2,\{1\}^{j-1-i},2,\{1\}^{k-j-1},2,\{1\}^{n-3-k}) \\
&=\frac{{n-3 \choose 3}}{(n-3)!}U_{p}(2,2,\{1\}^{n-5}) \equiv 0 \pmod{p}.
\end{split}
\end{displaymath}

Similarly, for the case $j=i+1$, replacing $u_{k}$ by $u_{k-1}+p$ for $k=i+1$ or $u_{k-2}+2p$ for $k>i+1$, we have
\begin{eqnarray}\label{minisum3}
&&\sum\limits_{\begin{smallmatrix}1 \le i <j \le n-2 \\j = i+1  \end{smallmatrix}}\sum\limits_{\begin{smallmatrix}u_1<\cdots <u_{n-1}<3p\\ u_{i+1}-u_{i}=u_{j+1}-u_{j}=p \end{smallmatrix}}\frac{1}{u_1u_2\cdots u_{n-1}}  \nonumber \\
&=&\sum\limits_{1\le i \le n-3}\sum\limits_{\begin{smallmatrix}u_1<\cdots <u_{n-3}<p \end{smallmatrix}}{\frac{1}{u_{1}\cdots u_{i}(u_{i}+p)(u_{i}+2p)(u_{i+1}+2p)\cdots (u_{n-3}+2p)}}  \nonumber\\
&\equiv& \sum\limits_{1\le i \le n-3}\sum\limits_{\begin{smallmatrix}u_1<\cdots <u_{n-3}<p \end{smallmatrix}}{\frac{1}{u_{1}\cdots u_{i}^{3}u_{i+1}\cdots u_{n-3}}\Big(1-\frac{3p}{u_{i}}-\frac{2p}{u_{i+1}}-\cdots -\frac{2p}{u_{n-3}}  \Big)   }  \nonumber\\
&\equiv &\frac{U_{p}(3,\{1\}^{n-4})}{(n-3)!}-3p\frac{U_{p}(4,\{1\}^{n-4})}{(n-3)!}-2p\sum\limits_{1\le i<j \le n-3}H_{p}(\{1\}^{i-1},3,\{1\}^{j-1-i},2,\{1\}^{n-3-j})  \nonumber\\
&\equiv &-2p\sum\limits_{1\le i<j \le n-3}H_{p}(\{1\}^{i-1},3,\{1\}^{j-1-i},2,\{1\}^{n-3-j}) \pmod{p^2}.
\end{eqnarray}

Substituting (\ref{aterm})--(\ref{minisum3}) into (\ref{T1start}), we obtain
\begin{equation}\label{T1final}
\begin{split}
T_{1}&\equiv \frac{U_{3p}(\{1\}^{n-1})}{(n-1)!}-3p \sum\limits_{1\le i <j \le n-3}\Big(H_{p}(\{1\}^{i-1},3,\{1\}^{j-i-1},2,\{1\}^{n-3-j})\\
&\qquad +H_{p}(\{1\}^{i-1},2,\{1\}^{j-i-1},3,\{1\}^{n-3-j}) \Big)\\
&\equiv -3p \frac{2{n-3 \choose 2}}{(n-3)!}U_{p}(2,3,\{1\}^{n-5}) \\
&\equiv 0 \pmod{p^2}.
\end{split}
\end{equation}

For the second sum in (\ref{Rn3start}), we have
\begin{eqnarray*}
T_{2}&=&\frac{1}{2p^2}\sum\limits_{i=2}^{n-4}\sum\limits_{j=i+2}^{n-2}\sum\limits_{\begin{smallmatrix}u_1<\cdots <u_{i-1}<p\\<u_{i+1}<\cdots <u_{j-1}<2p \\<u_{j+1}<\cdots <u_{n-1}<3p \end{smallmatrix}}{\frac{1}{u_1\cdots u_{i-1}u_{i+1}\cdots u_{j-1} u_{j+1}\cdots u_{n-1}}} \\
&=&\frac{1}{2p^2}\sum\limits_{i=2}^{n-4}\sum\limits_{j=i+2}^{n-2}H_{p}(\{1\}^{i-1})\Big(\sum\limits_{x_{1}<\cdots <x_{j-1-i}<p}{\frac{1}{(p+x_1)\cdots (p+x_{j-1-i})} }\Big) \\
&&\cdot \Big(\sum\limits_{y_1<\cdots <y_{n-1-j}<p}{\frac{1}{(2p+y_1)\cdots (2p+y_{n-1-j})}}\Big)\\
&\equiv &\frac{1}{2p^2}\sum\limits_{i=2}^{n-4}\sum\limits_{j=i+2}^{n-2}H_{p}(\{1\}^{i-1})\sum\limits_{x_1<\cdots <x_{j-1-i}<p}\frac{1}{x_1\cdots x_{j-1-i}}\Big(1-\frac{p}{x_1}-\cdots -\frac{p}{x_{j-1-i}}\Big) \\
&&\cdot \sum\limits_{y_1<\cdots<y_{n-1-j}<p}{\frac{1}{y_1\cdots y_{n-1-j}}\Big(1-\frac{2p}{y_1}-\cdots-\frac{2p}{y_{n-1-j}}\Big)} \\
&\equiv & \frac{1}{2p^2}\sum\limits_{i=2}^{n-4}\sum\limits_{j=i+2}^{n-2}H_{p}(\{1\}^{i-1})\Bigg(H_{p}(\{1\}^{j-1-i})-p\frac{U_{p}(2,\{1\}^{j-2-i})}{(j-2-i)!}\Bigg)\\
&&\cdot \Bigg(H_{p}(\{1\}^{n-1-j})-2p\frac{U_{p}(2,\{1\}^{n-2-j})}{(n-2-j)!}\Bigg)\\
&\equiv &\frac{1}{2p^2}\sum\limits_{i=2}^{n-4}\sum\limits_{j=i+2}^{n-2}H_{p}(\{1\}^{i-1})H_{p}(\{1\}^{j-1-i})H_{p}(\{1\}^{n-1-j}) \pmod{p^2}.
\end{eqnarray*}
Since $(i-1)+(j-1-i)+(n-1-j)=n-3$ is odd, so at least one of $i-1$, $j-1-i$ and $n-1-j$ must be odd. Hence we deduce that
\begin{equation}\label{SumT2}
T_{2}\equiv 0 \pmod{p^2}.
\end{equation}

For the third sum in (\ref{Rn3start}), we have
\begin{equation}\label{SumT31}
\begin{split}
T_{3}&=\frac{1}{p}\sum\limits_{j=2}^{n-2}\Bigg(\sum\limits_{u_1<\cdots <u_{j-1}<p}\frac{1}{u_1\cdots u_{j-1}}\Bigg)\\
&\qquad \cdot \Bigg(\sum\limits_{\begin{smallmatrix} x_1<\cdots <x_{n-1-j}<2p\\ x_2-x_1, \cdots, x_{n-1-j}-x_{n-2-j}\in \mathcal{P}_{p} \\x_{1}, x_{2}, \cdots, x_{n-1-j}\in \mathcal{P}_{p} \end{smallmatrix}}\frac{1}{(p+x_1)\cdots (p+x_{n-1-j})}\Bigg).
\end{split}
\end{equation}
Note that
\begin{eqnarray}\label{SumT32}
&&\sum\limits_{\begin{smallmatrix} x_1<\cdots <x_{n-1-j}<2p\\ x_2-x_1, \cdots, x_{n-1-j}-x_{n-2-j} \in \mathcal{P}_{p}\\x_{1}, x_{2},\cdots, x_{n-1-j}\in \mathcal{P}_{p} \end{smallmatrix}}\frac{1}{(p+x_1)\cdots (p+x_{n-1-j})} \nonumber\\
&=&\sum\limits_{\begin{smallmatrix}x_1<\cdots <x_{n-1-j}<2p\\x_1,\cdots,x_{n-1-j}\in\mathcal{P}_{p} \end{smallmatrix}}-\sum\limits_{a=1}^{n-2-j}\sum\limits_{\begin{smallmatrix}x_1<\cdots<x_{n-1-j}<2p \\x_{a+1}-x_{a}=p\\ x_1,\cdots, x_{n-1-j}\in \mathcal{P}_{p} \end{smallmatrix}}\frac{1}{(p+x_1)\cdots (p+x_{n-1-j})}  \nonumber\\
&\equiv& \sum\limits_{\begin{smallmatrix}x_1<\cdots <x_{n-1-j}<2p\\x_1,\cdots,x_{n-1-j}\in\mathcal{P}_{p} \end{smallmatrix}}{\frac{1}{x_1\cdots x_{n-1-j}}\Big(1-\frac{p}{x_1}-\cdots -\frac{p}{x_{n-1-j}}\Big)  }  \nonumber \\
&&-\sum\limits_{a=1}^{n-2-j}\sum\limits_{x_1<\cdots <x_{n-2-j}<p}{\frac{1}{x_1\cdots x_{a}(x_{a}+p)(x_{a+1}+p)\cdots (x_{n-2-j}+p)}} \nonumber \\
&\equiv& H_{2p}(\{1\}^{n-1-j})-\frac{p}{(n-2-j)!}U_{2p}(2,\{1\}^{n-2-j}) \nonumber \\
&&-\sum\limits_{a=1}^{n-2-j}\sum\limits_{x_1<\cdots <x_{n-2-j}<p}{\frac{1}{x_1\cdots x_{a-1}x_{a}^{2}x_{a+1}\cdots x_{n-2-j}}\Big(1-\frac{p}{x_a}-\cdots -\frac{p}{x_{n-2-j}} \Big) }  \nonumber\\
&\equiv& H_{2p}(\{1\}^{n-1-j})-\sum\limits_{a=1}^{n-2-j}\Big(H_{p}(\{1\}^{a-1},2,\{1\}^{n-2-j-a})-pH_{p}(\{1\}^{a-1},3,\{1\}^{n-2-j-a}) \nonumber\\
&& -p\sum\limits_{b=a+1}^{n-2-j}H_{p}(\{1\}^{a-1},2,\{1\}^{b-a-1},2,\{1\}^{n-2-j-b}) \Big)   \nonumber\\
&\equiv&\frac{U_{2p}(\{1\}^{n-1-j})}{(n-1-j)!}-\frac{U_{p}(2,\{1\}^{n-3-j})}{(n-3-j)!}+p\cdot \frac{U_{p}(3,\{1\}^{n-j-3})}{(n-3-j)!} \nonumber \\
&& +p\cdot \frac{{n-2-j \choose 2}}{(n-2-j)!}U_{p}(2,2,\{1\}^{n-4-j})  \nonumber\\
&\equiv& \frac{U_{2p}(\{1\}^{n-1-j})}{(n-1-j)!}-\frac{U_{p}(2,\{1\}^{n-3-j})}{(n-3-j)!} \pmod{p^2}.
\end{eqnarray}
Substituting (\ref{SumT32}) into (\ref{SumT31}), we obtain that
\begin{equation}\label{SumT3}
T_{3}\equiv \sum\limits_{\begin{smallmatrix}j=2 \\ j \,\, \mathrm{odd}\end{smallmatrix}}^{n-2}(n+1-j)\frac{B_{p-j}B_{p-n+j}}{j(n-j)}p \pmod{p^2}.
\end{equation}

In the same way, we can show that
\begin{equation}\label{SumT4}
T_{4}\equiv \frac{1}{2}\sum\limits_{\begin{smallmatrix}j=2 \\ j \,\, \mathrm{odd}\end{smallmatrix}}^{n-2}(j+1)\frac{B_{p-j}B_{p-n+j}}{j(n-j)}p \pmod{p^2}.
\end{equation}

Finally, substituting (\ref{T1final}), (\ref{SumT2}), (\ref{SumT3}) and (\ref{SumT4}) into (\ref{Rn3start}), we complete the proof of Theorem \ref{Rn3}.
\end{proof}
Similar to the proof of Corollary \ref{Sn2}, we can give the determination of $S_{n}^{(3)}(p)$ modulo $p$.
\begin{cor}\label{Sn3}
Let $n>4$ be an even integer and $p>n$ be a prime. We have
\begin{displaymath}
\sum\limits_{\begin{smallmatrix}l_1+l_2+\cdots +l_n=3p\\ l_1,l_2,\cdots ,l_n<p \end{smallmatrix}}\frac{1}{l_1l_2\cdots l_n} \equiv -\frac{n!}{6}\sum\limits_{\begin{smallmatrix} a=2 \\ a \,\mathrm{odd}\end{smallmatrix}}^{n-2}(n+a-3)\frac{B_{p-a}B_{p-n+a}}{a(n-a)} \pmod{p}.
\end{displaymath}
\end{cor}

\section{Proofs of Theorems \ref{thm1} and \ref{thm2} }
\begin{proof}[Proof of Theorem \ref{thm1}]
By (\ref{zhou}), Corollaries \ref{Sn2} and \ref{Sn3}, we obtain
\begin{equation}\label{S6pcong1}
S_{6}^{(1)}(p)\equiv 0 \pmod{p}, \quad S_{6}^{(2)}(p)\equiv 40B_{p-3}^{2} \pmod{p}, \quad S_{6}^{(3)}(p) \equiv -80B_{p-3}^{2} \pmod{p}.
\end{equation}
From \cite[Lemma 1(i)]{Wang2015}, for any $r\ge 1$ we have
\begin{equation}\label{internal}
S_{6}^{(1)}(p^r)\equiv S_{6}^{(5)}(p^r) \pmod{p^r}, \quad S_{6}^{(2)}(p^r)\equiv S_{6}^{(4)}(p^r) \pmod{p^r}.
\end{equation}
Therefore we have
\begin{equation}\label{S6pcong2}
S_{6}^{(4)}(p)\equiv 40B_{p-3}^{2} \pmod{p}, \quad S_{6}^{(5)}(p)\equiv 0 \pmod{p}.
\end{equation}

Suppose $1\le m \le 5$ and $r\ge 1$. For any 6-tuple $({{l}_{1}},\cdots ,{{l}_{6}})$ of integers satisfying
\[{{l}_{1}}+\cdots +{{l}_{6}}=m{{p}^{r+1}}, \quad 1 \le l_{i} <p^{r+1}, \quad l_{i} \in \mathcal{P}_{p}, \quad 1 \le i \le 6,\]
we rewrite them as
\[{{l}_{i}}={{x}_{i}}{{p}^{r}}+{{y}_{i}},\quad 0\le {{x}_{i}}<p,\quad 1\le {{y}_{i}}<{{p}^{r}}, \quad {{y}_{i}}\in {\mathcal{P}_{p}}, \quad 1 \le i \le 6.\]
Since
\[\big(\sum\limits_{i=1}^{6}{{{x}_{i}}}\big){{p}^{r}}+\sum\limits_{i=1}^{6}{{{y}_{i}}}=m{{p}^{r+1}},\]
we know there exists $1\le a \le 5$ such that
\[\left\{ \begin{array}{ll}
  {{x}_{1}}+\cdots +{{x}_{6}}=mp-a \\
 {{y}_{1}}+\cdots +{{y}_{6}}=ap^{r} \\
\end{array} \right..\]
By Lemma \ref{Cam}, we have $C_{a}^{(m)} \equiv 0$ (mod $p$) for $1\le a \le 5$. Hence for $1\le j \le 6$, we have
\[\sum_{\substack{x_1+\dots+x_6=mp-a \\ 0 \le x_i<p, 1 \le i \le 6}} x_j
= \frac{1}{6} \sum_{\begin{smallmatrix}x_1+\dots+x_6=mp-a \\ 0 \le x_{i} < p, 1 \le i \le 6\end{smallmatrix}} (x_1+x_2+\dots+x_6)
= \frac{mp-a}{6} C^{(m)}_{a} \equiv 0  \pmod{p}. \]

We have
\begin{eqnarray*}
   S_{6}^{(m)}({{p}^{r+1}}) & =&\sum\limits_{\begin{smallmatrix}
 {{l}_{1}}+\cdots +{{l}_{6}}=m{{p}^{r+1}} \\
 {{l}_{i}}\in {\mathcal{P}_{p}},{{l}_{i}}<{{p}^{r+1}}
\end{smallmatrix}}{\frac{1}{{{l}_{1}}{{l}_{2}}\cdots {{l}_{6}}}} \\
 & =&\sum\limits_{a=1}^{5}{\sum\limits_{\begin{smallmatrix}
 {{x}_{1}}+\cdots +{{x}_{6}}=mp-a \\
 0\le {{x}_{i}}<p
\end{smallmatrix}}{\sum\limits_{\begin{smallmatrix}
 {{y}_{1}}+\cdots +{{y}_{6}}=a{{p}^{r}} \\
 {{y}_{i}}\in {\mathcal{P}_{p}},{{y}_{i}}<{{p}^{r}}
\end{smallmatrix}}{\frac{1}{({{x}_{1}}{{p}^{r}}+{{y}_{1}})\cdots ({{x}_{6}}{{p}^{r}}+{{y}_{6}})}}}} \\
& \equiv &\sum_{a=1}^{5}\sum_{\begin{smallmatrix}
 x_1+\cdots +x_6=mp-a \\
 0\le x_i<p\end{smallmatrix}}\ \sum_{\begin{smallmatrix}
  y_1+\cdots +y_6=ap^r \\
 y_i\in \mathcal{P}_p,\, y_i<p^r\end{smallmatrix}} \left(1-\frac{x_1}{y_1}p^r-\cdots-\frac{x_6}{y_6}p^r\right)
 \frac1{ y_1\cdots y_6}   \\
 & \equiv &{{C}_{1}^{(m)}}S_{6}^{(1)}({{p}^{r}})+{{C}_{2}^{(m)}}S_{6}^{(2)}({{p}^{r}})+{{C}_{3}^{(m)}}S_{6}^{(3)}({{p}^{r}})\\
 && +{{C}_{4}^{(m)}}S_{6}^{(4)}({{p}^{r}}) +C_{5}^{(m)}S_{6}^{(5)}(p^r) \pmod{p^{r+1}}.
\end{eqnarray*}

Since $C_{a}^{(m)} \equiv 0$ (mod $p$), from (\ref{internal}) we deduce that
\begin{equation}\label{add1}
\begin{split}
   S_{6}^{(m)}({{p}^{r+1}})  &\equiv ({{C}_{1}^{(m)}}+{{C}_{5}^{(m)}})S_{6}^{(1)}({{p}^{r}})+({{C}_{2}^{(m)}}+{{C}_{4}^{(m)}})S_{6}^{(2)}({{p}^{r}}) \\    & \quad +C_{3}^{(m)}S_{6}^{(3)}(p^r) \pmod {p^{r+1}}.
\end{split}
\end{equation}

Let $r=1$ in (\ref{add1}). From (\ref{S6pcong1}), (\ref{S6pcong2}) and Lemma \ref{Cam} we deduce that
\begin{equation}\label{S6p2}
\begin{split}
S_{6}^{(1)}({{p}^{2}}) &\equiv \frac{2p}{5}S_{6}^{(1)}(p)-\frac{p}{10}S_{6}^{(2)}(p)+\frac{p}{30}S_{6}^{(3)}(p) \equiv -\frac{20}{3}pB_{p-3}^2 \pmod{p^2},\\
 S_{6}^{(2)}({{p}^{2}}) &\equiv -\frac{8p}{5}S_{6}^{(1)}(p)+\frac{2p}{5}S_{6}^{(2)}(p)-\frac{2p}{15}S_{6}^{(3)}(p) \equiv \frac{80}{3}pB_{p-3}^2 \pmod{p^2},\\
 S_{6}^{(3)}({{p}^{2}}) &\equiv \frac{12p}{5}S_{6}^{(1)}(p)-\frac{3p}{5}S_{6}^{(2)}(p)+\frac{p}{5}S_{6}^{(3)}(p) \equiv -40pB_{p-3}^2 \pmod{p^2}.
\end{split}
\end{equation}

Now suppose for some $r \ge 2$ we have
\begin{equation}\label{S6induction}
\begin{split}
S_{6}^{(1)}({{p}^{r}}) &\equiv  -\frac{20}{3}p^{r-1}B_{p-3}^2 \pmod{p^r}, \\
S_{6}^{(2)}({{p}^{r}})& \equiv \frac{80}{3}p^{r-1}B_{p-3}^2 \pmod{p^r}, \\
S_{6}^{(3)}({{p}^{r}})&\equiv  -40p^{r-1}B_{p-3}^2 \pmod{p^r}.
\end{split}
\end{equation}
From (\ref{add1}) and Lemma \ref{Cam} we deduce that
\begin{displaymath}
\begin{split}
S_{6}^{(1)}({{p}^{r+1}}) &\equiv \frac{2p}{5}S_{6}^{(1)}(p^r)-\frac{p}{10}S_{6}^{(2)}(p^r)+\frac{p}{30}S_{6}^{(3)}(p^r) \equiv -\frac{20}{3}p^rB_{p-3}^2 \pmod{p^{r+1}},\\
 S_{6}^{(2)}({{p}^{r+1}}) &\equiv -\frac{8p}{5}S_{6}^{(1)}(p^r)+\frac{2p}{5}S_{6}^{(2)}(p^r)-\frac{2p}{15}S_{6}^{(3)}(p^r) \equiv \frac{80}{3}p^rB_{p-3}^2 \pmod{p^{r+1}},\\
 S_{6}^{(3)}({{p}^{r+1}}) &\equiv \frac{12p}{5}S_{6}^{(1)}(p^r)-\frac{3p}{5}S_{6}^{(2)}(p^r)+\frac{p}{5}S_{6}^{(3)}(p^r) \equiv -40p^rB_{p-3}^2 \pmod{p^{r+1}}.
\end{split}
\end{displaymath}
By induction on $r$, we complete our proof.
\end{proof}

\begin{proof}[Proof of Theorem \ref{thm2}]
Let $n=mp^{r}$, where $p$ does not divide $m$.
For any 6-tuple $({{l}_{1}},\cdots ,{{l}_{6}})$ of integers satisfying ${{l}_{1}}+\cdots +{{l}_{6}}=n$, $l_{i} \in \mathcal{P}_{p}$, $1 \le i \le 6$, we rewrite them as
\[{{l}_{i}}={{x}_{i}}{{p}^{r}}+{{y}_{i}}, \quad x_{i} \ge 0, \quad 1\le {{y}_{i}}<{{p}^{r}}, \quad {{y}_{i}}\in {\mathcal{P}_{p}}, \quad 1 \le i \le 6.\]
Since
\[\big(\sum\limits_{i=1}^{6}{{{x}_{i}}}\big){{p}^{r}}+\sum\limits_{i=1}^{6}{{{y}_{i}}}=m{{p}^{r}},\]
we know there exists $1\le a \le 5$ such that
\[\left\{ \begin{array}{ll}
  {{x}_{1}}+\cdots +{{x}_{6}}=m-a \\
 {{y}_{1}}+\cdots +{{y}_{6}}=ap^{r} \\
\end{array} \right..\]
For $1 \le a \le 5$, the equation $x_{1}+x_{2}+x_{3}+x_{4}+x_{5}+x_{6}=m-a$ has $\binom {m-a+5}{5}$ solutions $(x_{1},x_{2},x_{3},x_{4},x_{5},x_{6})$ of nonnegative integers. Hence
\begin{equation}\label{thm2rec}
\begin{split}
 & \quad \sum\limits_{\begin{smallmatrix}
 {{l}_{1}}+\cdots +{{l}_{6}}=m{{p}^{r}} \\
  {{l}_{1}},\cdots ,{{l}_{6}}\in {\mathcal{P}_{p}}
\end{smallmatrix}}{\frac{1}{{{l}_{1}}{{l}_{2}}\cdots {{l}_{6}}}} \\
 & =\sum\limits_{a=1}^{5}{\sum\limits_{\begin{smallmatrix}
 {{x}_{1}}+\cdots +{{x}_{6}}=m-a \\
  {{x}_{i}} \ge 0
\end{smallmatrix}}{\sum\limits_{\begin{smallmatrix}
 {{y}_{1}}+\cdots +{{y}_{6}}=a{{p}^{r}} \\
 {{y}_{i}}\in {\mathcal{P}_{p}},{{y}_{i}}<{{p}^{r}}
\end{smallmatrix}}{\frac{1}{({{x}_{1}}{{p}^{r}}+{{y}_{1}})\cdots ({{x}_{6}}{{p}^{r}}+{{y}_{6}})}}}} \\
&\equiv \sum\limits_{a=1}^{5}{\binom{m-a+5}{5}S_{6}^{(a)}(p^r)} \pmod{p^{r}}.
\end{split}
\end{equation}

According to $r=1$ or $r \ge 2$, we split our proof into two cases.

(i) If $r=1$, then from (\ref{S6pcong1}), (\ref{S6pcong2}), and (\ref{thm2rec}) we obtain
\begin{displaymath}
\begin{split}
&\quad \sum\limits_{\begin{smallmatrix}
 {{l}_{1}}+{{l}_{2}}+\cdots +{{l}_{6}}=n \\
 {{l}_{1}},\cdots ,{{l}_{6}}\in {\mathcal{P}_{p}}
\end{smallmatrix}}{\frac{1}{{{l}_{1}}{{l}_{2}}{{l}_{3}}{{l}_{4}}{{l}_{5}}l_{6}}} \\
 & \equiv \bigg(\binom{m+4}{5}+\binom{m}{5}\bigg)S_{6}^{(1)}(p) +\bigg(\binom{m+3}{5}+\binom{m+1}{5}\bigg)S_{6}^{(2)}(p) \\
&\quad +\binom{m+2}{5}S_{6}^{(3)}(p) \\
 & \equiv \frac{20}{3}(m^{3}-m)B_{p-3}^{2}  \pmod{p}.
 \end{split}
\end{displaymath}
Since $m=\frac{n}{p}$, we complete the proof of (i).

(ii) If $r \ge 2$, since we have already proved that (\ref{S6induction}) is true for all $r\ge 2$, we deduce that  $S_{6}^{(2)}({{p}^{r}})\equiv -4S_{6}^{(1)}(p^{r})$ (mod ${{p}^{r}}$) and $S_{6}^{(3)}(p^r) \equiv 6 S_{6}^{(1)}(p^r)$ (mod $p^r$). Again by (\ref{internal}), we have
\[ S_{6}^{(5)}({{p}^{r}}) \equiv S_{6}^{(1)}({{p}^{r}})  \pmod{{{p}^{r}}}, \quad  S_{6}^{(4)}({{p}^{r}})  \equiv S_{6}^{(2)}({{p}^{r}}) \equiv  -4S_{6}^{(1)}(p^{r})\pmod{{p}^{r}}.\]
Hence from (\ref{thm2rec}) we obtain
\begin{displaymath}
\begin{split}
 &\quad \sum\limits_{\begin{smallmatrix}
 {{l}_{1}}+{{l}_{2}}+\cdots +{{l}_{6}}=n \\
 {{l}_{1}},\cdots ,{{l}_{6}}\in {\mathcal{P}_{p}}
\end{smallmatrix}}{\frac{1}{{{l}_{1}}{{l}_{2}}{{l}_{3}}{{l}_{4}}{{l}_{5}}l_{6}}} \\
 & \equiv \bigg(\binom{m+4}{5}+\binom{m}{5}-4\binom{m+3}{5}-4\binom{m+1}{5}+6\binom{m+2}{5}\bigg)S_{6}^{(1)}(p^{r}) \\
 & \equiv m S_{6}^{(1)}(p^{r})  \pmod{p^{r}}.
 \end{split}
\end{displaymath}
By Theorem \ref{thm1}, we have $S_{6}^{(1)}(p^{r}) \equiv -\frac{5!}{18}p^{r-1}B_{p-3}^{2}$ (mod $p^{r}$). This completes the proof of (ii).
\end{proof}

\section{Concluding Remarks}

As some examples for Theorems \ref{Rn2} and \ref{Rn3}, for any prime $p>10$ we have

\[R_{4}^{(2)}(p) \equiv 0 \pmod{p}, \quad R_{4}^{3}(p) \equiv 0 \pmod{p}, \]
\[R_{6}^{(2)}(p) \equiv \frac{6!}{18}B_{p-3}^{2} \pmod{p},   \quad R_{6}^{(3)}(p)\equiv \frac{2\cdot 6!}{9}B_{p-3}B_{p-5} \pmod{p},\]
\[R_{8}^{(2)}(p) \equiv \frac{8!}{15}B_{p-3}B_{p-5} \pmod{p},  \quad R_{8}^{(3)}(p) \equiv \frac{8!}{3}B_{p-3}B_{p-5} \pmod{p}, \]
\[R_{10}^{(2)}(p) \equiv \frac{10!}{1050}\big(50B_{p-3}B_{p-7}+21B_{p-5}^{2}\big) \pmod{p}, \]
\[R_{10}^{(3)}(p)\equiv \frac{10!}{175}\big(50B_{p-3}B_{p-7}+21B_{p-5}^{2}\big)  \pmod{p}.\]

If we wan to apply the method in \cite{Wang2015} and the current paper to find the modulo $p^r$ determination for $S_{n}(p^r)$, the key step is to find the modulo $p$ determination of $S_{n}^{(m)}(p)$ for $1\le m < n$. Since $S_{n}^{(1)}(p)$ has already been figured out by (\ref{zhou}), we only need to consider the case $m\ge 2$. Equivalently, Lemma \ref{RSrelation} tells us that we only need to determine $R_{n}^{(m)}(p)$ modulo $p$ for $m\ge 2$. However, as $m$ increases, the computation becomes much more complicated. This is the major obstacle in solving the cases for large number of variables. Nevertheless, from the process of determining $R_{n}^{(2)}(p)$ and $R_{n}^{(3)}(p)$ presented in this paper, we believe that for any integers $n\ge 3$ and $1\le m < n$, there exists some rational numbers $c_{a_1,a_2,\cdots,a_k}$ and $d_{a_1,a_2,\cdots,a_k}$ such that for any prime $p>n+2$ and $r\ge 2$, we have
\[R_{n}^{(m)}(p) \equiv \sum\limits{c_{a_1,a_2,\cdots,a_k}B_{p-a_{1}}B_{p-a_{2}}\cdots B_{p-a_{k}} \pmod{p},}\]
\[S_{n}(p^r) \equiv p^{r-1}\sum\limits{d_{a_1,a_2,\cdots,a_k}B_{p-a_{1}}B_{p-a_{2}}\cdots B_{p-a_{k}} \pmod{p^{r}},}\]
here the sum on the right hand side runs over all possible tuples $(a_1,a_2,\cdots, a_k)$ of odd integers such that
\[a_1+a_2+\cdots +a_k=n.\]

For $n\le 6$, this has already been verified from Theorems \ref{Rn2} and \ref{Rn3} and the work of \cite{WangCai,Wang2015}. For $n\ge 7$, we are able to guess what the congruence should look like. For example, the congruences for $n=8,9$ should be of the forms
\[S_{8}(p^r) \equiv ap^{r-1}B_{p-3}B_{p-5} \pmod{p}, \quad a\in \mathbb{Q},\]
\[S_{9}(p^r) \equiv p^{r-1}\big(b_{1}B_{p-3}^3+b_{2}B_{p-9}\big) \pmod{p^r}, \quad b_{1},b_{2} \in \mathbb{Q}.\]