\section{Introduction}
\label{sec:intro}

Ultraperipheral collisions (UPCs) of relativistic ions correspond to large
impact parameters between the nuclei exceeding the sum of their radii, so
that short-range strong interactions between the ions are suppressed and
reactions proceed rather via the emission of quasi-real photons by the
colliding ions. Thus, UPCs allow one to study photon-photon and photon-hadron
(proton, nucleus) interactions at high energies
\cite{Baltz:2007kq}.
During the last decade, UPCs have become an active field of research, driven
by experimental results obtained at the Relativistic Heavy Ion Collider (RHIC)
and the Large Hadron Collider (LHC), for a recent experimental review see,
e.g., \cite{Klein:2017vua}. Notable examples of various UPC processes and
their analyses include the two-photon production of dilepton pairs
\cite{Klein:2018cjh,Goncalves:2018gca}; light-by-light scattering $\gamma
\gamma \to \gamma \gamma$ and searches for potential physics beyond the
Standard Model \cite{dEnterria:2013zqi,Fichet:2014uka,Klusek-Gawenda:2016euz};
an electromagnetic double-scattering contribution to dimuon pair production
in photon-photon scattering \cite{vanHameren:2017krz}; exclusive
photoproduction of charmonia in proton-proton \cite{Aaij:2013jxj,Aaij:2014iea},
proton-nucleus \cite{TheALICE:2014dwa} and nucleus-nucleus
\cite{Abbas:2013oua,Abelev:2012ba,Adam:2015sia,Khachatryan:2016qhq} UPCs and
of bottomonia in proton-proton \cite{Aaij:2015kea} and proton-nucleus UPCs
\cite{Chudasama:2016eck}; new constraints on the small-$x$ gluon distribution
in the proton \cite{Jones:2013pga,Guzey:2013qza} and heavy nuclei
\cite{Adeluyi:2012ph,Guzey:2013xba} and the dynamics of strong interactions at
high energies in the color dipole framework \cite{Lappi:2013am,%
Goncalves:2014wna,Mantysaari:2017dwh}; and exclusive photoproduction of $\rho$
mesons on nuclei \cite{Adler:2002sc,Abelev:2007nb,Agakishiev:2011me,%
Adam:2015gsa} as well as tests of models of nuclear shadowing
\cite{Rebyakova:2011vf,Frankfurt:2015cwa}.

Focusing on UPC studies of nuclear structure in QCD at the LHC, coherent
$J/\psi$ photoproduction in Pb-Pb UPCs at $\sqrt{s_{NN}} = 2.76$ TeV
\cite{Abbas:2013oua,Abelev:2012ba,Adam:2015sia,Khachatryan:2016qhq} revealed
a significant nuclear suppression of the measured rapidity distributions. In
the framework of the leading logarithmic approximation of perturbative QCD
\cite{Ryskin:1992ui}, it can be interpreted as evidence of large nuclear
gluon shadowing, $R_g=f_{g/A}(x,\mu^2)/[Af_{g/N}(x,\mu^2)] \approx 0.6$ at $x=
10^{-3}$ and $\mu^2=3$ GeV$^2$ ($f_{g/A}$ and $f_{g/N}$ are gluon densities in
Pb and the proton, respectively). This value of $R_g$ agrees with predictions
of the leading twist nuclear shadowing model \cite{Frankfurt:2011cs}, which
are characterized by small theoretical uncertainties in this kinematic region.
It is also broadly consistent with the EPS09 \cite{Eskola:2009uj}, nCTEQ15
\cite{Kovarik:2015cma}, and EPPS16 \cite{Eskola:2016oht} nuclear parton
distribution functions (nPDFs), which however have significant uncertainties
in this kinematic regime. Note that in the collinear factorization framework,
next-to-leading order (NLO) perturbative QCD corrections to the cross section
of $J/\psi$ photoproduction are large \cite{Ivanov:2004vd,Jones:2015nna} and
the relation between the gluon parton distribution function (PDF) and the
gluon generalized parton distribution (GPD) is model-dependent, which makes it
challenging to interpret the UPC data on $J/\psi$ photoproduction on nuclei in
terms of the NLO gluon nPDF.

The program of UPC measurements continues with Run 2 at the LHC, where besides
photoproduction of vector mesons, inclusive dijet photoproduction in Pb-Pb
UPCs $AA \to A+{\rm 2 jets}+X$ has also recently been measured by the ATLAS
collaboration \cite{Atlas} (for leading-order QCD predictions for rates of
this process, see \cite{Strikman:2005yv}). The cross section of this process
is sensitive to quark and gluon nPDFs $f_{j/A}(x,\mu^2)$ in a wide range of
the momentum fraction $x$ and the resolution scale $\mu > {\cal O}(20)$ GeV,
where one still expects sizable nuclear modifications of the PDFs. In addition,
imposing the requirement that the target nucleus stays intact, one can study
diffractive dijet photoproduction in UPCs $AA \to A+{\rm 2 jets}+X+A$. Studies
of this process may shed some light on the mechanism of QCD factorization
breaking in diffractive photoproduction and, for the first time, give access
to nuclear diffractive PDFs \cite{Guzey:2016tek,Basso:2017mue}. While further progress in
constraining nPDFs will benefit from studies of high-energy hard processes
with nuclei in proton-nucleus ($pA$) scattering at the LHC
\cite{Salgado:2011wc} and lepton-nucleus ($eA$) scattering at a future
Electron-Ion Collider (EIC) \cite{Accardi:2012qut} and LHeC
\cite{AbelleiraFernandez:2012cc}, UPCs at the LHC present an important and
complimentary method of obtaining new constraints on nPDFs in a wide
kinematic range already now.

In this work, we make predictions for the cross section of inclusive dijet
photoproduction in Pb-Pb UPCs at the LHC using NLO perturbative QCD
\cite{Klasen:2002xb} and nCTEQ15 nPDFs. We show that our approach provides a
good description
of various cross section distributions measured
by the ATLAS collaboration \cite{Atlas}.
Our analysis also shows that the dijet photoproduction cross
section in the considered kinematics is sensitive to nuclear modifications of
the PDFs. As a function of the momentum fraction $x_A$, the ratio of the cross
sections calculated with nPDFs and in the impulse approximation behaves
similarly to $R_g$ for a given $\mu$ and deviates from
unity by $10-20$\% for the central nCTEQ15 fit. The calculations using EPPS16
nPDFs and predictions of the leading twist nuclear shadowing model give
similar results. This suggests that inclusive dijet photoproduction on nuclei
can be used to reduce uncertainties in the determination of nPDFs, which are
currently significant and comparable in size to the magnitude of the
calculated nuclear modifications of the dijet photoproduction cross section.

The remainder of this paper is structured as follows. In Sec.
\ref{sec:formalism}, we outline the formalism of dijet photoproduction
in UPCs using NLO perturbative QCD. We present and discuss our results for
the LHC in Sec.\ \ref{sec:results} and draw conclusions in Sec.\
\ref{sec:conclusions}.

\section{Photoproduction of dijets in UPCs in NLO perturbative QCD}
\label{sec:formalism}

Typical leading-order (LO) Feynman diagrams for dijet photoproduction in UPCs
of nuclei $A$ and $B$ are shown in Fig.\ \ref{fig:photo_jets_upc}, where the
graphs ($a$) and ($b$) correspond to the direct and resolved photon
contributions, respectively. Note that beyond LO, the separation of the direct
and resolved photon contributions depends on the factorization scheme and
scale
(see the discussion below).
\begin{figure}[t]
\begin{center}
\epsfig{file=photo_jets_upc.eps,scale=0.8}
\caption{Typical leading-order Feynman graphs for dijet photoproduction in UPCs
 of hadrons $A$ and $B$. Graphs (a) and (b) correspond to the direct and
 resolved photon contributions, respectively.}
\label{fig:photo_jets_upc}
\end{center}
\end{figure}

Using the Weizs\"acker-Williams method, which allows one to treat the electromagnetic field of an ultra-relativistic
ion as a flux of equivalent quasi-real photons~\cite{Baltz:2007kq,Bertulani:2005ru}, and
the collinear factorization framework for photon-nucleus scattering,
the cross section of the UPC
process $AB\to A+{\rm 2 jets}+X$ is given by \cite{Klasen:2002xb}
\begin{eqnarray}
&& d\sigma(AB \to A+{\rm 2 jets}+X) = \nonumber\\
&& \sum_{a,b} \int^{y_{\rm max}}_{y_{\rm min}} dy \int^1_0 dx_{\gamma} \int^{x_{A,{\rm max}}}_{x_{A,{\rm min}}} dx_A f_{\gamma/A}(y)f_{a/\gamma}(x_{\gamma},\mu^2)f_{b/B}(x_A,\mu^2) d\hat{\sigma}(ab \to {\rm jets})\,,
\label{eq:cs}
\end{eqnarray}
where $a,b$ are parton flavors; $f_{\gamma/A}(y)$ is the flux of equivalent photons
emitted by ion $A$, which depends on the photon light-cone momentum fraction $y$;
$f_{a/\gamma}(x_{\gamma},\mu^2)$ is the PDF of the photon, which depends on the
momentum fraction $x_{\gamma}$ and the factorization scale $\mu$; $f_{b/B}(x_A,\mu^2)$
is the nuclear PDF with $x_A$ being the corresponding parton momentum fraction;
and $d\hat{\sigma}(ab \to {\rm jets})$ is the elementary cross section for production
of two- and three-parton final states emerging as jets in hard scattering of partons
$a$ and $b$. The sum over $a$ involves quarks and gluons for the resolved photon contribution
and the photon for the direct photon contribution dominating at $x_{\gamma} \approx 1$.
At LO, the direct photon contribution has support exactly only at $x_{\gamma} = 1$,
i. e., $f_{a_/\gamma}=\delta(1-x_\gamma)$.
At NLO, the virtual and real corrections are calculated with massless
quarks in dimensional regularization, ultraviolet (UV) divergences are renormalized in
the $\overline{\rm MS}$ scheme, infrared (IR) divergences are canceled and factorized
into the proton and photon PDFs, respectively. For the latter, this implies
a transformation from the DIS$_\gamma$ into the $\overline{\rm MS}$
scheme.
The integration limits are determined by the rapidities and transverse momenta
 of the produced jets, see Sec.~\ref{sec:results}.
Note that Eq.~(\ref{eq:cs}) is based on
the clear separation of scales, which characterize the long-distance electromagnetic interaction
and the short-distance strong interaction.
It generalizes
the NLO perturbative QCD formalism
of collinear factorization for jet photoproduction in lepton-proton
scattering developed in Refs.~\cite{Klasen:1995ab,Klasen:1996it,Klasen:1997br,Klasen:2002xb}, which
successfully described HERA $ep$ data on dijet photoproduction~\cite{Klasen:2011ax}.
Hence, Eq.~(\ref{eq:cs}) involves universal nuclear PDFs $f_{b/B}(x_A,\mu^2)$, which can be accessed
in a variety of hard processes involving nuclear targets~\cite{Eskola:2009uj,Kovarik:2015cma,Eskola:2016oht},
and the universal photon PDFs $f_{a/\gamma}(x_{\gamma},\mu^2)$, which are determined by $e^{+}e^{-}$ data, for a review, see \cite{Klasen:2002xb}.
Hence, the interplay between the direct and resolved photon contributions in Eq.~(\ref{eq:cs}) is also universal
and controlled by the standard $\mu^2$ evolution equations of photon PDFs and the choice of the factorization scheme.

In our analysis, we used the following input for Eq.~(\ref{eq:cs}).
For photon PDFs $f_{a/\gamma}(x_{\gamma},\mu^2)$, we used
the GRV HO parametrization~\cite{Gluck:1991jc}, which we transformed from the DIS$_{\gamma}$ to the $\overline{{\rm MS}}$ factorization scheme. These photon PDFs have been tested profoundly at HERA and the Large Electron-Positron (LEP) collider at CERN and are very robust,
in particular at high $x_{\gamma}$ (dominated by the pQCD photon-quark splitting), which is correlated with the low-$x_A$ gluons
and sea quarks in Pb that present one of the points of interest of the present study.
For nuclear PDFs $f_{b/B}(x_A,\mu^2)$, we employed the nCTEQ15 parametrization~\cite{Kovarik:2015cma}.
The photon flux $f_{\gamma/A}(y)$ produced by a relativistic point-like charge $Z$ is given by
 the standard expression
 \begin{equation}
 f_{\gamma/A}(y)=\frac{2 \alpha_{\rm e.m.}Z^2}{\pi}\frac{1}{y} \left[\zeta K_0(\zeta)K_1(\zeta)-\frac{\zeta^2}{2}(K_1^2(\zeta)-K_0^2(\zeta)) \right] \,,
 \label{eq:flux}
 \end{equation}
where $\alpha_{\rm e.m.}$ is the fine-structure constant; $K_{0,1}$ are modified Bessel functions of the second kind;
$\zeta=y m_p b_{\rm min}$ with $m_p$ being the proton mass and $b_{\rm min}$ the minimal distance between two nuclei.
For Pb-Pb UPCs, Eq.~(\ref{eq:flux}) with $b_{\rm min}=14.2$ fm reproduces very well the photon flux calculated taking into account the nuclear form factor and the suppression of strong interactions at impact parameters $b < b_{\rm min}$, see the discussion in~\cite{Nystrand:2004vn}.

The NLO calculation of the dijet photoproduction cross section using Eq.~(\ref{eq:cs}) is numerically implemented in an NLO
parton-level Monte Carlo~\cite{Klasen:1995ab,Klasen:1996it,Klasen:1997br,Klasen:2002xb}, which has been successfully tested in many
different environments (HERA, LEP, Tevatron). It implements the anti-k$_T$ algorithm (but we have at most two partons in the jet)
and all the kinematic conditions and cuts used in the ATLAS analysis~\cite{Atlas} that are explicitly explained in the
following section. Hadronization corrections and underlying event (UE) subtractions are not part of our analysis, but they are expected
to be performed with PYTHIA simulations by the experiment once the data are final (as has been done at HERA).

\section{Predictions for dijet photoproduction in Pb-Pb collisions at the LHC}
\label{sec:results}

The main goal of the present paper is the first NLO QCD calculation of the cross section of inclusive dijet photoproduction
in Pb-Pb UPCs and the conclusion whether it can describe the results of the ATLAS measurement~\cite{Atlas}.
The ATLAS analysis was performed using
the following conditions and selection criteria:
\begin{itemize}
\item
the anti-k$_T$ algorithm with the jet radius $R=0.4$;
\item
  the leading jet has $p_{T,1} > 20$ GeV, while the other jets have a different cut on
  $p_{T,i\neq1}> 15$ GeV as
required \cite{Klasen:1995xe}, which corresponds to $35 < H_T < 400$ GeV, where $H_T=\sum_ip_{T,i}$;
\item
all jets have rapidities $|\eta_i| < 4.4$;
\item
the combined mass of all reconstructed jets is $ 35 < m_{\rm jets} < 400$ GeV;
\item
the parton momentum fraction on the photon side $z_{\gamma}=y x_{\gamma}$, $10^{-4} < z_{\gamma} < 0.05$;
\item
the parton momentum fraction on the nucleus side $x_A$, $5 \times 10^{-4} < x_A < 1$.
\end{itemize}

The ATLAS results are presented as distributions in terms of the total jet transverse momentum $H_T=\sum_ip_{T,i}$
and the photon $z_{\gamma}$ and nucleus $x_A$ light-cone momentum fractions
\begin{equation}
z_{\gamma}=\frac{m_{\rm jets}}{\sqrt{s_{NN}}}e^{y_{\rm jets}} \quad , \quad
x_{A}=\frac{m_{\rm jets}}{\sqrt{s_{NN}}}e^{-y_{\rm jets}} \,,
\label{eq:estimators}
\end{equation}
where
\begin{equation}
m_{\rm jets}=\left[\left(\sum_i E_i\right)^2-\left|\sum_i \vec{p}_i \right|^2 \right]^{1/2} \quad , \quad
y_{\rm jets}=\frac{1}{2} \ln \left(\frac{\sum_i E_i+p_{z,i}}{\sum_i E_i-p_{z,i}} \right) \,.
\label{eq:estimators2}
 \end{equation}
 In Eqs.~(\ref{eq:estimators2}), the index $i$ runs over all measured jets; $E_i$ and $\vec{p}_i$ denote the jet energy and momentum, respectively.
 Note that at LO, the kinematics of $2 \to 2$ parton scattering and the momentum fractions $z_{\gamma}$ and $x_A$
  can be exactly reconstructed from the dijet measurement.
 At NLO, Eqs.~(\ref{eq:estimators}) serve as hadron-level estimators of the momentum fractions entering Eq.~(\ref{eq:cs});
 for brevity, we use the same notations in Eqs.\ (\ref{eq:cs}) and (\ref{eq:estimators2}).

\begin{figure}[t]
\begin{center}
\epsfig{file=HT_unc_final_full_new.eps,scale=1.15}
\caption{NLO QCD predictions for the cross section of dijet photoproduction in Pb-Pb UPCs at $\sqrt{s_{NN}}=5.02$ TeV in the ATLAS kinematics as a function of $H_T$ for different bins of $x_A$. The central values and the corresponding shaded
uncertainty bands are obtained using nCTEQ15 nPDFs.
The crosses are the ATLAS data points that we extracted from \cite{Atlas}.}
\label{fig:HT}
\end{center}
\end{figure}

Figures~\ref{fig:HT}, \ref{fig:xA}, \ref{fig:zgamma}, and \ref{fig:xA2} show our results for the cross section of dijet photoproduction in Pb-Pb UPCs
in the ATLAS kinematics (see above) as a function of $H_T$, $x_A$, and $z_{\gamma}$
for different bins of these variables. They correspond to Figs.\ 12-15 of Ref.\ \cite{Atlas}.
In each bin, our predictions are obtained using the central fit of nCTEQ15 nPDFs~\cite{Kovarik:2015cma}. The shaded bands
quantify the uncertainty of our results $\Delta \sigma$ due to the uncertainty of nCTEQ15 nPDFs. It is calculated by adding in quadrature the individual uncertainties corresponding to each of 32 error sets
\begin{equation}
\Delta \sigma =\frac{1}{2} \sqrt{ \sum_{k=1,\ {\rm odd}}^{31} \left(\sigma(f_k)-\sigma(f_{k+1})\right)^2} \,,
\label{eq:error}
\end{equation}
where $\sigma(f_k)$ is the cross section calculated using the $f_k$ nCTEQ15 error nPDFs.

\begin{figure}[t]
\begin{center}
\epsfig{file=xA_unc_final_full_new.eps,scale=1.15}
\caption{NLO QCD predictions for the cross section of dijet photoproduction in Pb-Pb UPCs at $\sqrt{s_{NN}}=5.02$ TeV in the ATLAS kinematics as a function of $x_A$ for different bins of $H_T$.
The crosses are the ATLAS data points that we extracted from \cite{Atlas}. }
\label{fig:xA}
\end{center}
\end{figure}

\begin{figure}[t]
\begin{center}
\epsfig{file=zgamma_unc_final_full.eps,scale=1.15}
\caption{NLO QCD predictions for the cross section of dijet photoproduction in Pb-Pb UPCs at $\sqrt{s_{NN}}=5.02$ TeV in the ATLAS kinematics as a function of $z_{\gamma}$ for different bins of $H_T$.}
\label{fig:zgamma}
\end{center}
\end{figure}

\begin{figure}[t]
\begin{center}
\epsfig{file=xA2_unc_final_full.eps,scale=1.15}
\caption{NLO QCD predictions for the cross section of dijet photoproduction in Pb-Pb UPCs at $\sqrt{s_{NN}}=5.02$ TeV in the ATLAS kinematics as a function of $x_A$ for different bins of $z_{\gamma}$.}
\label{fig:xA2}
\end{center}
\end{figure}

A comparison of our results shown in Figs.~\ref{fig:HT}, \ref{fig:xA}, \ref{fig:zgamma}, and \ref{fig:xA2}
with Figs.~12-15 of Ref.\ \cite{Atlas} demonstrates that our calculations describe well
the corresponding distributions.
To illustrate this point, in Figs.~\ref{fig:HT} and \ref{fig:xA}, which present the phenomenologically important
distributions in $H_T$ and $x_A$,
respectively, we also explicitly show by crosses the ATLAS data points, which we extracted from \cite{Atlas} using the WebPlotDigitizer
tool~\cite{Digitizer}. One can readily see from these figures that the results of our calculations describe both the shape and
normalization of the data rather well. The description of the remaining two distributions is also adequate.
Note that the ATLAS data is preliminary and has not been corrected (unfolded) for the detector response.

\begin{figure}[t]
\begin{center}
\epsfig{file=xA_dopol.eps,scale=0.7}
\caption{NLO QCD predictions for the ratio of the cross section of dijet photoproduction
to that calculated in IA in Pb-Pb UPCs at $\sqrt{s_{NN}}=5.02$ TeV in the ATLAS kinematics
as a function of $x_A$ for different $H_T$ bins.
 The shaded bands give the uncertainty of nCTEQ15 nPDFs. }
\label{fig:xA_dopol}
\end{center}
\end{figure}

To assess the impact of measurements of the inclusive dijet photoproduction cross section on nPDFs,
we focus on the $x_A$ distribution.
Our results are shown in Figs.~\ref{fig:xA_dopol} and \ref{fig:xA3}. In Fig.~\ref{fig:xA_dopol}, we show
the ratio of the cross section calculated using nCTEQ15 nPDFs in lead to the
one calculated in the impulse approximation (IA), where nuclear PDFs are assumed not to include any nuclear modifications and are
given by the weighted sum of free proton and neutron PDFs,
$f_{b/A}^{\rm IA}=Zf_{b/p}+(A-Z)f_{b/n}$.
The panels in this figure corresponds to nine bins in $H_T$ presented in Fig.~\ref{fig:xA} (the numerator of the presented
ratio is given by the curves in Fig.~\ref{fig:xA}).
One can see that the cross section ratio as a function of $x_A$
behaves similarly to the ratio $R_g = f_{g/A}(x,\mu^2)/[Af_{g/N}(x,\mu^2)]$ of the nuclear and nucleon gluon distributions.
It dips below unity for $x_A < 0.01$ due to nuclear shadowing and then becomes enhanced around $x_A=0.1$ due to the assumed gluon
antishadowing. For $x_A > 0.3$, the cross section ratio shows again a suppression due to
the EMC effect encoded in the nPDFs.
While the scaling violations (the $H_T$ dependence) of the shown ratios are difficult to see, they seem to be positive
at small $x \sim 0.01$ and negative at large $x$, as they should
be, cf.~\cite{Kovarik:2015cma,Eskola:2016oht}.

Note that in spite of large values of the resolution scale probed in the considered kinematics, $\mu > {\cal O}(20)$ GeV,
one can see that one is still sensitive to nuclear modifications of the PDFs at the $10-20$\% level for the central value of our predictions.
One should also note that the uncertainty due to nPDFs, which is given by the shaded bands,
is significant and comparable to the size of the discussed nuclear modifications.
This can be viewed as an opportunity to reduce uncertainties of nPDFs using data on cross section of inclusive dijet photoproduction in nuclei in global QCD fits of nPDFs.

In Fig.~\ref{fig:xA3}, we show our predictions for the dijet cross section as a function of $x_A$
integrated over $H_T$ and $z_{\gamma}$.
The top panel presents separately the resolved (green, dot-dashed) and the direct
(blue, dashed) photon contributions to the cross section as well as their sum (red, solid).
As can be expected, because of the
correlation
between $x_A$ and $z_{\gamma}$, see Eq.~(\ref{eq:estimators}), the resolved photon contribution dominates for $x_A > 0.01$. We find that for small $x_A < 0.01$, the two contributions are comparable with
the direct contribution being somewhat larger.
While this behavior is qualitatively similar to the results of the LO analysis in the framework of PYTHIA 8 with EPPS16 nPDFs~\cite{Helenius:2018bai}, the relative contribution of the resolved photon term is larger at NLO, but this statement depends of course on the choice of the photon factorization scheme and scale.

\begin{figure}[t]
\begin{center}
\epsfig{file=xA3_unc_final_full.eps,scale=1.0}
\caption{NLO QCD predictions for the cross section of dijet photoproduction in Pb-Pb UPCs at $\sqrt{s_{NN}}=5.02$ TeV in the ATLAS kinematics as a function of $x_A$. Top: The resolved (green, dot-dashed) and direct
(blue, dashed) photon contributions and their sum (red, solid). Middle: The ratio to the impulse approximation. Bottom: The ratio of cross sections calculated using the nCTEQ15 and EPPS16 nPDFs.
The shaded bands show the uncertainty of nCTEQ15 nPDFs.}
\label{fig:xA3}
\end{center}
\end{figure}

The middle panel of Fig.~\ref{fig:xA3} presents the ratio to the impulse approximation.
Similarly to the trend already observed in Fig.~\ref{fig:xA_dopol} and discussed above,
this ratio as a function of $x_A$ behaves similarly to $R_g$.
This behavior is similar to the one observed in the case of dijet photoproduction in the kinematics of an Electron-Ion Collider
(EIC)~\cite{Klasen:2018gtb}.

Note that in the future and if/when experimental data become available, to study nuclear modifications of nPDFs one can directly
form the ratio of dijet cross sections measured in Pb-Pb and proton-proton ($pp$)
UPCs as a function of $x_A$. While the systematics are highly correlated between bins in $x_A$, the information on nPDFs is
in the shape of the cross section ratio, see  the middle panel of Fig.~\ref{fig:xA3} and its discussion above.
Also, while the central Pb-Pb and $pp$ collisions are very different, Pb-Pb and $pp$ UPCs have comparable multiplicities.
Thus, one can expect that the systematic uncertainties largely cancel in the nucleus-to-proton cross section ratio.

Finally, the bottom panel of Fig.~\ref{fig:xA3} presents the ratio of the dijet cross section
calculated using nCTEQ15 nPDFs to the one calculated with the central value of EPPS16 nPDFs. The shaded band
quantifies the uncertainty of the nCTEQ15 fit. One can see from the panel that the two parameterizations
of nPDFs give similar predictions, which differ by at most 5\% for all but one values of $x_A$.
We have also explicitly checked that the use of nPDFs calculated in the model of leading twist nuclear
shadowing~\cite{Frankfurt:2011cs} gives similarly close predictions for the dijet photoproduction cross section.

In our calculations, following the standard prescription for setting the hard scale in QCD calculations, we used
$\mu = 2 E_{T,1}$ in Eq.~(\ref{eq:cs}). In detail, we performed calculations
using $\mu = (E_{T,1}/4, E_{T,1}/2, E_{T,1}, 2\,E_{T,1},4\,E_{T,1})$ both at NLO and LO and found that (i) the integrated cross section of inclusive dijet photoproduction at NLO as a function of $\mu$ is approximately
constant is the vicinity of $\mu = 2 E_{T,1}$, (ii) while the NLO cross section slightly increases with an increase of $\mu$ up to $2E_{T,1}$ and then starts to decrease again,
the LO cross section steeply decreases monotonically, and (iii) the values of the two cross sections are close around $\mu = 2 E_{T,1}$.
Therefore, $\mu = 2 E_{T,1}$ in Eq.~(\ref{eq:cs}) corresponds to the choice, which is most numerically stable
against higher-order corrections.

In this work, we used the framework of collinear factorization and NLO perturbative QCD
to examine the sensitivity of the dijet photoproduction cross section to nuclear modifications of PDFs.
Alternatively, one can use this process to look for signs of the
BFKL and gluon saturation dynamics in the high-energy ($k_{T}$) factorization approach \cite{Kotko:2017oxg}.

\section{Conclusions}
\label{sec:conclusions}

In this work, we calculated the cross section of inclusive dijet photoproduction in Pb-Pb UPCs at the LHC
using NLO perturbative QCD and nCTEQ15 nPDFs.
We showed that our approach provides a good description of various cross section distributions measured by the ATLAS collaboration.
We found that the calculated dijet photoproduction cross section is sensitive to nuclear modifications of the PDFs. In particular,
as a function of the nucleus momentum fraction $x_A$,
the ratio of the cross sections calculated with nPDFs and in the impulse approximation
 behaves similarly to $R_g$ for given $\mu$ and deviates from unity by $10-20$\% for the central nCTEQ15 fit.
The calculations using EPPS16 nPDFs and predictions of the leading
twist nuclear shadowing model give similar results.
Therefore, inclusive dijet photoproduction on nuclei has the potential to reduce
uncertainties in determination of nPDFs, which are comparable to the magnitude
of the calculated nuclear modifications of the dijet photoproduction cross section.
Our present analysis is a step in this direction.

\acknowledgments

VG would like to thank A.~Angerami, I.~Helenius, K.~Kovarik, H.~Paukkunen, and M.~Strikman  for useful discussions
and the Institut f\"ur Theoretische Physik, Westf\"alische Wilhelms-Universit\"at M\"unster for hospitality.
VG's research is supported in part by
RFBR, research project 17-52-12070. The authors gratefully acknowledge financial
support of DFG through the grant KL 1266/9-1 within the framework of the joint German-Russian project ``New constraints on nuclear parton distribution
functions at small $x$ from dijet production in $\gamma A$ collisions at the LHC".


\section{Introduction}
Systems of an antikaon and few nucleons are of considerable interest in few-body physics because the low-energy $\bar{K}N$ interaction is considered to be sufficiently attractive to generate a quasibound state, the $\Lambda(1405)$, below $\bar{K}N$ threshold~\cite{Dalitz:1959dn,Dalitz:1960du,Dalitz:1967fp}. One of the recent theoretical achievements is the determination of the precise meson-baryon scattering amplitude, including the pole position of the $\Lambda(1405)$~\cite{Tanabashi:2018oca}. Currently, a whole set of experimental data near the $\bar{K}N$ threshold is successfully described in the framework of chiral SU(3) coupled-channels dynamics~\cite{Ikeda:2011pi,Ikeda:2012au}.

In order to study the $\bar{K}$ few-nucleon systems with established techniques for rigorous few-body calculations, it is desirable to construct a local meson-baryon potential. A detailed strategy for constructing such a local potential equivalent to chiral SU(3) dynamics has been developed in Ref.~\cite{Hyodo:2007jq}. It is also shown that one can renormalize the effects of $\pi Y$ ($\equiv\pi\Sigma, \pi\Lambda$) channels to obtain a single-channel $\bar{K}N$ potential. By applying this strategy to the scattering amplitude in Refs.~\cite{Ikeda:2011pi,Ikeda:2012au}, an elaborate single-channel $\bar{K}N$ potential has been derived in Ref.~\cite{Miyahara:2015bya}. This potential is applied to the investigations of few-body kaonic nuclei~\cite{Ohnishi:2017uni} and kaonic deuterium~\cite{Hoshino:2017mty}.

While the single-channel $\bar{K}N$ potential exactly reproduces the two-body scattering amplitude, in the application to the few-body systems, dynamical effects of the coupled-channels (such as $\pi Y$ plus nucleons) are not properly included. For states near the $\bar{K}N$ threshold, the neglect of such components does not affect observables very much. However, the results in Ref.~\cite{Ohnishi:2017uni} suggest that the quasibound state of $\bar{K}$ and six nucleons appears near the threshold of $\pi \Sigma$ plus five nucleons. In such cases, the explicit treatment of coupled-channels dynamics is expected to give sizable contributions.

Here we construct the coupled-channels $\bar{K}N$-$\pi\Sigma$-$\pi\Lambda$ potential~\cite{Miyahara:2018onh} equivalent to the scattering amplitude in Refs.~\cite{Ikeda:2011pi,Ikeda:2012au}. As an application within the two-body sector, we study the meson-baryon fractions in the $\Lambda(1405)$ by analyzing the wave function of the quasibound state.

\section{Coupled-channels meson-baryon potential}
We start from the coupled-channels Schr\"odinger equation for a nonrelativistic energy $E$,
\begin{align}
\left[ -\frac{\nabla^2}{2\mu_i}\delta_{ij} + \Delta M_i\, \delta_{ij} + V^\text{equiv}_{ij}(\bm{r},E) \right]\psi_{j}(\bm{r}) = E \psi_{i}(\bm{r})~,
\label{eq:Sch_psi}
\end{align}
where $\mu_{i}$, $\Delta M_{i}$, and $\psi_{i}$ are the reduced mass, the energy difference from the reference threshold, and the wave function in channel $i$, respectively. The coupled-channels potential $V_{ij}^{\text{equiv}}$ has a matrix form in channel basis with $i=\bar{K}N$, $\pi\Sigma$, and $\pi\Lambda$. In contrast to the single-channel $\bar{K}N$ potential which has a complex strength
due to the absorption into the $\pi Y$ channels, the strengths of $V_{ij}^{\text{equiv}}$ are real valued. In addition, $V_{ij}^{\text{equiv}}$ has an energy dependence, which originates from the energy dependence of the chiral interaction and from the renormalization of higher energy channels such as $\eta\Lambda$ and $K\Xi$.

Our task is to determine $V^\text{equiv}_{ij}$ such that the scattering amplitude from Eq.~\eqref{eq:Sch_psi} reproduces the one in Refs.~\cite{Ikeda:2011pi,Ikeda:2012au}. The potential $V^\text{equiv}_{ij}$ should be related to the interaction kernel in chiral SU(3) dynamics, but there is no direct way of converting the interaction kernel, given that the framework is different. We therefore introduce several matching conditions of $V^\text{equiv}_{ij}$ and the interaction kernel, not only on the real energy axis but also in the complex energy plane, in order to systematically determine $V^\text{equiv}_{ij}$. For the details of the potential construction procedure, see Ref.~\cite{Miyahara:2018onh}.

For practical applications, it is useful to represent $V_{ij}^{\text{equiv}}$ in a conveniently parametrized form. We use the following parametrization:
\begin{align}
V_{ij}^\text{equiv}(\bm{r},E) &= e^{-r^2(1/2b_{i}^2+1/2b_{j}^2)}
\sum_{\alpha=0}^{\alpha_{\rm max}} K_{\alpha,{ij}}\,
\left(\frac{E}{100 \text{ MeV}}\right)^\alpha\, .
\label{eq:Vfit}
\end{align}
The spatial distribution is assumed to be a Gaussian form with the range parameter $b_{i}$. The energy dependence is parametrized by a polynomial of order $\alpha_{\rm max}$ with real coefficients $K_{\alpha,{ij}}$. These parameters are determined by the matching conditions mentioned above, and the results are tabulated in Ref.~\cite{Miyahara:2018onh}. One remarkable fact is that the energy dependence of $V_{ij}^\text{equiv}$ can be accurately parametrized by a second-order polynomial ($\alpha_{\rm max}=2$). This is in sharp contrast to the single-channel version in Ref.~\cite{Miyahara:2018onh} where one needs $\alpha_{\rm max}$ as large as 10 for sufficient accuracy. In other words, the explicit inclusion of the $\pi Y$ channels permits a parametrization of the potential strength with a more natural, much weaker energy dependence.

\section{Compositeness of the $\Lambda(1405)$}
In general, a resonance eigenstate is expressed by a pole of the scattering amplitude in the complex energy plane. In the case of the $\Lambda(1405)$, the resonance is expressed by two poles~\cite{Jido:2003cb,Tanabashi:2018oca}. From the coupled-channels potential, we find poles at
$\sqrt{s} = 1424-27i$ MeV and $\sqrt{s} = 1380-81i$
with $\sqrt{s}=E+M_{N}+m_{K}$, consistently with the original scattering amplitude~\cite{Ikeda:2011pi,Ikeda:2012au}. The wave function of the $\Lambda(1405)$ can be obtained by evaluating Eq.\,\eqref{eq:Sch_psi} at these pole energies. Starting from the wave function it is possible to extract the properties of the resonance. It should be kept in mind that due to the unstable nature of resonances, the expectation value of operators (such as norm of the state) becomes complex.

Here we discuss the compositeness $X_{i}$ and the ``elementarity'' $Z$~\cite{Hyodo:2013nka} which are related to the norm of the wave function in channel $i$:
\begin{align}
  X_{i}&=\bra{\psi^{\dag}_{i}}\kket{\psi_{i}} , \quad
  Z=1-\sum_{i}X_{i},
\end{align}
where $\bra{\psi^{\dag}_{i}}$ is a left eigenstate $\bra{\psi^{\dag}_{i}}H=\bra{\psi^{\dag}_{i}}E$. When we deal with resonances, $\bra{\psi_{i}^{\dag}}$ is different from $\bra{\psi_{i}}$. The compositeness $X_{i}$ represents the fraction of the channel $i$ component in the total wave function. When the potential is energy dependent, the  sum of $X_{i}$s does not become unity, and the elementarity $Z$ is introduced. Using the Feshbach projection method, $Z$ can be expressed by the expectation value of the energy derivative of the potential~\cite{Miyahara:2018onh}. Originally, $Z$ has been introduced as a field renormalization constant of the bare state representing the elementary component of the composite system~\cite{Weinberg:1965zz}. In the present context, it is understood as the contributions not explicitly included in the model space.

\begin{table}[b]
\caption{Compositeness $X_{i}$ and elementarity $Z$ for high-mass pole of the $\Lambda(1405)$ coupled-channels system.}
\label{tbl:compositeness}
\begin{center}
\begin{tabular}{lccc}
\hline
Method & $X_{\pi\Sigma}$ & $X_{\bar{K} N}$ & $Z$\\[2pt]
\hline
Coupled-channels potential & $-0.02-0.25i$ & \ $1.01-0.13i$ \ & $0.01+0.37i$ \\
Single-channels potential~\cite{Miyahara:2015bya}  &               & \ $1.01-0.07i$ \ &  \\
Residue of the pole~\cite{Sekihara:2014kya} & $-0.19-0.22i$ & \ $1.14+0.01i$ \ & $0.05+0.21i$ \\[2pt]
\hline
\end{tabular}
\end{center}
\end{table}

We show the results of $X_{\pi\Sigma}$, $X_{\bar{K}N}$, and $Z$ for high-mass pole at $1424-27i$ MeV in Table~\ref{tbl:compositeness}. The results are consistent with those obtained by the single-channel potential\,\cite{Miyahara:2015bya}, as well as with the evaluation by the residue of the pole~\cite{Sekihara:2014kya}. In all cases, $X_{\bar{K}N}$ is close to unity and the others are almost zero. Hence, we conclude the $\bar{K}N$ dominance of the high-mass pole of the $\Lambda(1405)$.

\section{Summary}
We have constructed a coupled-channels meson-baryon ($\bar{K}N\leftrightarrow \pi Y$) potential based on chiral SU(3) dynamics. The evaluation of the compositeness tells us that the high-mass pole of the $\Lambda(1405)$ is dominated by the $\bar{K}N$ component. This potential will be useful in order to shed new light on meson-baryon few-body systems with strangeness.


\section{INTRODUCTION}

Binocular rivalry is the striking phenomenon that ensues when the two eyes view markedly different stimuli. The observer perceives only one stimulus at a time, and perception alternates between the two stimuli at irregular intervals\cite{levelt}. The perceived durations of the images are stochastic and uncorrelated with previous perceived durations\cite{kalar}. One of the initial studies conducted on this phenomenon was by Levelt in 1965. His four propositions went onto become the cornerstone of many of the mathematical models developed later on the same. They can be broadly summarised as:
\begin{enumerate}
    \item Increasing stimulus strength for one eye will increase the
perceptual predominance of that eye’s stimulus.
    \item Increasing stimulus strength for one eye will not affect the
average perceptual dominance duration of that eye’s stimulus.
Instead, it will reduce the average perceptual dominance duration of the other eye’s stimulus.
    \item Increasing stimulus strength for one eye will increase the
perceptual alternation rate.
    \item Increasing stimulus strength in both eyes while keeping it
equal between eyes will increase the perceptual alternation rate.
\end{enumerate}
This switching of dominance between the two perceptions is primarily due to the inhibition the dominant side exerts over the other. This reciprocal inhibition, is constantly accompanied by a slow negative feedback exerted by the suppressed side. This negative feedback acts through spike frequency adaptation and synaptic depression. \\
Spike frequency adaptation is a result of a decrease in the activity of the dominant side, eventually leading to the suppressed neuronal population expressing itself. Whereas synaptic depression results when the strength of the connectivity between neurons has reduced, thereby leading to decreased inhibitory effects of the dominant population on the suppressed one. These principles have been incorporated into developing various mathematical models, such as  \cite{kalar}, \cite{wilson}, \cite{laing}, \cite{main}, \cite{blake}, \cite{gberg}, \cite{lago}, \cite{lehky}, \cite{mats} and \cite{stol}. This article proceeds to study the dynamics predicted by Wilson's model, Laing and Chow's model and it's adaptation only variant as well as Kalarickal and Marshall's model.

\section{WILSON'S MODEL}
One of the distinguishing features of Wilson's model is the existence of separate neuronal populations to exert the inhibitory effect. That is, the neuronal population responsible for perception is physically different from the neuronal population exerting the inhibition on the suppressed side. The rivalry dynamics are modelled as below:

\begin{equation*}
    \tau\dot{E_1} = -E_1 + \frac{100(V_1-gI_2)_{+}^{2}}{(10+H_1)^2+(V_1-gI_2)_{+}^{2}}
\end{equation*}
\begin{equation*}
    \tau_{H}\dot{H_1} = -H_1 + hE_1
\end{equation*}
\begin{equation*}
    \tau_{I}\dot{I_1} = -I_1 + E_1
\end{equation*}
\begin{equation*}
    \tau\dot{E_2} = -E_2 + \frac{100(V_2-gI_1)_{+}^{2}}{(10+H_2)^2+(V_2-gI_1)_{+}^{2}}
\end{equation*}
\begin{equation*}
    \tau_{H}\dot{H_2} = -H_2 + hE_2
\end{equation*}
\begin{equation*}
    \tau_{I}\dot{I_2} = -I_2 + E_2
\end{equation*}
$V_i$ are the input stimuli, $g$ represents the cross inhibition parameter - it scales the strength of the inhibition exerted by the inhibitory neuronal population, $E_i$ represents the population firing rate, an indication of the activity of the neuronal population (higher firing rate implies greater dominance in perception), $I_i$ is the inhibitory firing rate of the $i$th inhibition population on the other, $H_i$ is the adaptation variable, $\tau$, $\tau_H$ and $\tau_I$ correspond to the time constant of the firing rate of the excitatory populations, adaption process and firing rate of the inhibitory population respectively \cite{wilson}. The asymptotic firing rates is determined by Naka-Rushton function, where  $(V_1-gI_2)_{+} = (V_1-gI_2)$ if $(V_1-gI_2) > 0$ else $V_1-gI_2)_{+} = 0$.
\\The above equations are simulated on MATLAB for various stimulus strengths and these parametric values: $g$ = 0.45, $\tau$ = 20$ms$, $\tau_H$ = 900$ms$ (adaptation is a relatively slow process), $\tau_I$ = 11$ms$ and $h$ = 0.47. The simulation results are as follows:

\begin{figure}[H]
\centering
\includegraphics[width=8cm, height=2cm]{WilBR1.jpg}\vspace{2mm}
\includegraphics[width=8cm, height=2cm]{WilBR2.jpg}\vspace{2mm}
\includegraphics[width=8cm, height=2cm]{WilBR3.jpg}\vspace{2mm}
\includegraphics[width=8cm, height=2cm]{WilBR4.jpg}\vspace{2mm}
\includegraphics[width=8cm, height=2cm]{WilBR5.jpg}
\caption{Wilson Model: Variation of perceptual alternation with stimulus strength.}
\label{fig:pic7}
\end{figure}

At high stimulus strengths, followed by an initial transient it is seen that both the stimuli equally dominate. There is an absence of alternation. However as the stimuli strength decreases, there is an onset of oscillations. The mean dominance time of these oscillations vary with stimulus strength. Increasing the stimulus strength decreases the mean dominance time of these oscillations, which is in compliance with the fourth proposition of Levelt, which states that increasing the stimulus strength while keeping them equal, increases the perceptual alternation rate. Upon further reducing the stimulus strength, there is a clear dominance of stimuli over the other, solely based on the initial conditions. There is no oscillation. This is commonly referred to as the "winner take all" regime. At an even lower stimulus strength, there is again an onset of perceptual oscillations. However these oscillations are different as when compared to the second regime. As the stimulus strength increases there is an increase in the mean dominance time of these oscillations, completely in contrast to what was stated in the fourth proposition. Upon further reduction in stimuli strength, it is seen that there is a fusion of these stimuli perceptions again. Hence it is evident that there is a strong absence of monotonous relationship between stimulus strength and the perceptual alternations. This is in strong contradiction of the initial fourth proposition of Levelt.
Although these simulations are for a certain set of parametric values, the nonmonotonicity of this behaviour is robust across many such feasible parametric values \cite{main}. \\
\vspace{-2mm}
\section{LAING AND CHOW'S MODEL}
In contrast to the above Wilson's model, Laing and Chow's dynamics describes a single neuronal population to exert the inhibitory effect and result in perception. Unlike other models, Laing and Chow's model also describes the evolution of adaptation and synaptic depression as first order equations.
\begin{equation*}
    \dot{u_1} = -u_1 + f(\alpha u_1g_1 - \beta u_2g_2 - a_1 +I_1)
\end{equation*}
\begin{equation*}
    \tau_{a}\dot{a_1} = -a_1 + \phi_a f(\alpha u_1g_1 - \beta u_2g_2 - a_1 +I_1)
\end{equation*}
\begin{equation*}
    \tau_{d}\dot{g_1} = 1 - g_1 - g_1\phi_d f(\alpha u_1g_1 - \beta u_2g_2 - a_1 +I_1)
\end{equation*}
\begin{equation*}
    \dot{u_2} = -u_2 + f(\alpha u_2g_2 - \beta u_1g_1 - a_2 +I_2)
\end{equation*}
\begin{equation*}
    \tau_{a}\dot{a_2} = -a_2 + \phi_a f(\alpha u_2g_2 - \beta u_1g_1 - a_2 +I_2)
\end{equation*}
\begin{equation*}
    \tau_{d}\dot{g_2} = 1 - g_2 - g_2\phi_d f(\alpha u_2g_2 - \beta u_1g_1 - a_2 +I_2)
\end{equation*}
$I_i$ are the input stimuli, $\beta$ represents the cross inhibition parameter, $u_i$ represents the average population firing rate, an indication of the activity of the neuronal population, $g_i$ is the synaptic depression, $a_i$ is the adaptation variable, $\tau_a$ and $\tau_d$ correspond to the time constant of adaption and synaptic depression processes respectively \cite{laing}. The gain function $f$ is taken to be the Heaviside step function, i.e., $f(x) = 1$ for $x ≥ 0$ and $f(x) = 0$ for $x < 0$. $\alpha$, $\phi_a$ and $\phi_d$ are scaling parameters. All the dynamic variables are normalised to hold values between $0$ and $1$. \\
The following simulation results are obtained for these parametric values: $\alpha = 0.35$, $\phi_a = 0.6$, $\phi_d = 0.6$, $\beta = 0.7$, $\tau_a$ = 20$ms$ and $\tau_d$ = 40$ms$.

\begin{figure}[H]
\centering
\includegraphics[width=8cm, height=2cm]{LCBR1.jpg}\vspace{1.5mm}
\includegraphics[width=8cm, height=2cm]{LCBR2.jpg}\vspace{1.5mm}
\includegraphics[width=8cm, height=2cm]{LCBR3.jpg}\vspace{1.5mm}
\includegraphics[width=8cm, height=2cm]{LCBR4.jpg}\vspace{1.5mm}
\includegraphics[width=8cm, height=2cm]{LCBR5.jpg}
\label{fig:pic7}
\caption{Laing and Chow's model: Variation of perceptual alternation with stimulus strength. }
\end{figure}

These dynamics do not predict fusion at low stimuli strength or equal dominance at high stimuli strengths. However, there is a winner take all regime (winner determined from the initial conditions) present at intermediate input values. Once again, at high stimuli strengths, as the strength is increased there is a decrease in the mean dominance time while at low stimuli strengths, there is an increase in the mean dominance time. This behaviour is in accordance with Wilson model predictions as well, and once again in contradiction of Levelt's fourth proposition.
\\
To specifically study the effect of adaptation and depression independently on excitatory firing rates of the neuronal populations, two models were derived from Laing and Chow's initial model, namely the adaptation and the depression model \cite{main}.

\subsection{Laing and Chow's Adaptation Model}
\vspace{-2mm}
\begin{equation*}
    \dot{u_1} = -u_1 + f(-\beta u_2 - ga_1 +I_1)
\end{equation*}
\begin{equation*}
    \tau_{a}\dot{a_1} = -a_1 + u_1
\end{equation*}
\begin{equation*}
    \dot{u_2} = -u_2 + f(-\beta u_1 - ga_2 +I_2)
\end{equation*}
\begin{equation*}
    \tau_{a}\dot{a_2} = -a_2 + u_2
\end{equation*}

\subsection{Laing and Chow's Depression model}
\vspace{-2mm}
\begin{equation*}
    \dot{u_1} = -u_1 + f(- \beta u_2g_2 +I_1)
\end{equation*}
\begin{equation*}
    \tau_{d}\dot{g_1} = 1 - g_1 - \gamma u_1g_1
\end{equation*}
\begin{equation*}
    \dot{u_2} = -u_2 + f(- \beta u_1g_1 +I_2)
\end{equation*}
\begin{equation*}
    \tau_{d}\dot{g_2} = 1 - g_2 - \gamma u_2g_2
\end{equation*}

However, in the above two models, $f$ is no longer a Heaviside step function, but is a sigmoid instead. That is $f(x) = 1/{1 + exp[−(x − \theta)/k]}$, where $1/k$ is its slope and $\theta$ is its threshold.\\
The adaptation model is simulated for these parameters: $k = 0.1$, $\theta = 0.2$, $g = 0.5$, $\tau_a$ = 100$ms$, and $\beta = 0.9$. Here $g$ is a scaling parameter of the adaptation on neuronal firing rate.

\begin{figure}[H]
\centering
\includegraphics[width=8cm, height=2cm]{LCABR1.jpg}\vspace{2mm}
\includegraphics[width=8cm, height=2cm]{LCABR2.jpg}\vspace{2mm}
\includegraphics[width=8cm, height=2cm]{LCABR3.jpg}
\end{figure}
\begin{figure}[H]
\centering
\includegraphics[width=8cm, height=2cm]{LCABR4.jpg}
\caption{Laing and Chow's Adaptaion model: Variation of perceptual alternation with stimulus strength.}
\label{fig:pic7}
\end{figure}

The simulation results of the adaptation model adhere to the previously predicted nonmonotonicity. However, there is an absence of the winner take all regime. This model predicts oscillations for all input stimuli strength. \\
Wilson's model and Laing's models were relatively new in explaining the process of binocular rivalry. A slightly older model is the one proposed by Kalarickal and Marshall. This model shares a few fundamental principles with the above models, however its predictions are quite different.
\vspace{2mm}
\section{Kalarickal and Marshall's Model}
This model was proposed prior to Wilson's and Laing and Chow's model. This model accounts for adaptation but not synaptic depression and has properties common with Laing's model. $x_i$ is the activation level, $I_i$ is the input stimulus, $W_i^{+}$ is the synaptic weight of the excitatory input pathway to neuron $i$, $W_{ij}^{-}$ is the synaptic weight of the lateral inhibitory pathway from neuron $i$ to neuron $j$, and $y_{ij}$ is the adaptation level of that pathway \cite{kalar}. It is important to notice that this model incorporates noise as well, that is $b(t)$.
\begin{equation*}
    \dot{x_1} = -x_1 + (1 - x_1)W_1^+I_1 - (c_1 + x_1)W_{21}^-y_{21}max(x_2,0)
\end{equation*}
\begin{equation*}
    \dot{y_{21}} = c_2((1-y_{21}) - c_3max(x_2,0)W_{21}^-y_{21}) + b(t)
\end{equation*}
\begin{equation*}
    \dot{x_2} = -x_2 + (1 - x_2)W_2^+I_2 - (c_1 + x_2)W_{12}^-y_{12}max(x_1,0)
\end{equation*}
\begin{equation*}
    \dot{y_{12}} = c_2((1-y_{12}) - c_3max(x_1,0)W_{12}^-y_{12}) + b(t)
\end{equation*}
The noise variable $b(t)=m$ if $r_{ij}(t) < p$ and $b(t) = -m$ otherwise, where $r_{ij}(t)$ is a uniformly distributed random variable in [0,1), $p$ is a positive constant in [0,1], and $m$ is a small positive constant. For the purpose of simulating the above equations the values of the parameters were taken as $W_1^+ = W_2^+ = 0.25$, $c_1 = 0.01$, $c_2 = 0.008$, $c_3 = 0.083$, $W_{12}^- = W_{21}^- = 250$, $p = 0.5$ and $m = 0.0025$.

\begin{figure}[H]
\centering
\includegraphics[width=8cm, height=2cm]{kalarBR1.jpg}\vspace{2mm}
\includegraphics[width=8cm, height=2cm]{KalarBR2.jpg}
\end{figure}
\begin{figure}[H]
\centering
\includegraphics[width=8cm, height=2cm]{KalarBR3.jpg}\vspace{2mm}
\includegraphics[width=8cm, height=2cm]{KalarBR4.jpg}\vspace{2mm}
\includegraphics[width=8cm, height=2cm]{kalarBR5.jpg}
\caption{Kalarickal and Marshall's model: Variation of perceptual alternation with stimulus strength.}
\label{fig:pic7}
\end{figure}

Kalarickal's model possesses a network architecture similar to that of Laing and Chow's model, that is, a single neuronal population being responsible for visual perception as well as exerting inhibitory effect on the other population. However, its predictions are quite contrary to the above two models. This model predicts no fusion at any stimuli strength. However, at low values there is a dominance of one perception over the other, once again subject to initial conditions. Also, these oscillations are not as uniform as in the previous cases, due to noise being incorporated into the modelling. It also leads to significantly long transition period prior oscillations. However, the trend of the mean dominance time with stimulus strength is robust in the presence of noise as well. It is evident from the simulations, that as stimuli strength increases there is a decrease in the mean dominance time. And this trend is persistent till extremely low stimuli strengths, in agreement with Levelt's fourth proposition.
\\
Although the models have some differences in their trends, there seems to be an agreement over the existence of oscillations, whose behaviour complies with Levelt's fourth proposition in some range of stimulus strength. The reasons for these oscillations is considered in some detail in the next section, along with discussion on the remaining three of Levelt's propositions.
\vspace{2mm}

\section{DISCUSSION}

When the dominant side is inhibiting the suppressed population from expressing, it is constantly experiencing a slow negative feedback reducing its dominance. And simultaneously, the suppressed population experiences a form of feedback recovery. One possible cause of this flipping in dominance can result from when the negative feedback on the dominant side increases significant enough for its activity to drop below that of the suppressed population, there by letting the latter express itself. This phenomenon is known as "release". In such a phenomenon it is likely to observe a more significant dip in the activity of the active population than a surge in the activity of the suppressed one. Another possible mechanism is known as the "escape" mechanism. This occurs when the feedback recovery of the suppressed population slowly grows enough to result in the overall input becoming positive, thereby resulting in switching of perception. This results in a significant rise in the activity of the suppressed population as when compared to the drop in the activity of the dominant side \cite{main}. The behaviour of these populations just before switching of the dominance is heavily dictated by the gain function incorporated into the modelling, which in the above models has been either Heaviside step function, the sigmoid function or the Naka-Rushton function. The reason for the not allowing negative values for these gain functions, is its resulting repercussion on the population firing rate. Negative values of the gain function results in negative population firing rates, which physically has no significance.
\\
A salient feature to be observed in the above models is the presence or absence of recurrent excitation. Wilson's, Laing's Adaptation and Kalarickal's models do not incorporate any recurrent excitation where as Laing and Chow's model does. It's evident, that even in its absence, the predictions of Wilson's model are coherent with Laing's to a decent degree. Upon incorporating recurrent excitation into Wilson's model, there was not much change observed in its qualitative behaviour, there by leading to the conclusion that recurrent excitation isn't necessary for the existence of oscillations. In fact, it has been previously studied that it results in system instability and predicts unreported behaviour.\\
Another important consideration is oscillation variation with different parameters. It has been previously found that the qualitative behaviour of these simulations doesn't get altered with differences in the time constants of the slow processes such as adaptation (in case of Wilson's model \cite{main}).  \\
The most important parameter is the one that dictates the strength of cross inhibition. The following simulations were conducted with varying values of the cross inhibition parameter for the same stimulus strength for two models, namely Wilson's and Laing and Chow's.\\

\textit{Laing and Chow's model}
\vspace{-3mm}
\begin{figure}[H]
\centering
\includegraphics[width=8cm, height=2cm]{LCCP1.jpg}\vspace{2mm}
\includegraphics[width=8cm, height=2cm]{LCCP2.jpg}\vspace{2mm}
\includegraphics[width=8cm, height=2cm]{LCCP3.jpg}
\end{figure}
\begin{figure}[H]
\centering
\includegraphics[width=8cm, height=2cm]{LCCP4.jpg}
\caption{Laing and Chow's model: Variation of perceptual alternation with cross inhibition parameter}
\label{fig:pic7}
\end{figure}

\textit{Wilson's Model}
\begin{figure}[h!]
\centering
\includegraphics[width=8cm, height=2cm]{WCP1.jpg}\vspace{2mm}
\includegraphics[width=8cm, height=2cm]{WCP2.jpg}\vspace{2mm}
\includegraphics[width=8cm, height=2cm]{WCP3.jpg}\vspace{2mm}
\includegraphics[width=8cm, height=2cm]{WCP4.jpg}
\caption{Wilson's model: Variation of perceptual alternation with cross inhibition parameter}
\label{fig:pic7}
\end{figure}\\
Both the models present coherent results with respect to the variation of mean dominance time with cross inhibition parameter. As its value is reduced, there is a reduction in the mean dominance time of the perceptual alternations. This phenomenon can be explained from the release and escape perspective presented at the beginning of this section. Lower the value of the cross inhibition parameter, lower the inhibition to be overcome by the suppressed population and lower the negative feedback required by the dominant side to get suppressed. This behaviour is monotonic across multiple stimuli strengths. \\
\\
So far the discussion has been with respect to binocular rivalry in response to stimuli of equal strength. In case of the Wilson model, when simulated for asymmetric stimuli the following is obtained.
\begin{figure}[H]
\centering
\includegraphics[width=8cm, height=2cm]{LVL1.jpg}
\end{figure}
\begin{figure}[H]
\centering
\includegraphics[width=8cm, height=2cm]{LVL2.jpg}\vspace{2mm}
\includegraphics[width=8cm, height=2cm]{LVL3.jpg}\vspace{2mm}
\includegraphics[width=8cm, height=2cm]{LVL4.jpg}
\caption{Wilson's model: Variation of perceptual alternation with asymmetrical stimuli strength.}
\label{fig:pic7}
\end{figure}
It can be seen that as the strength of one stimulus is increased with respect to the other, the dominance duration of the stronger stimulus increases, while that of the constant stimuli remains the same. This leads to an overall decrease in the rate of the perceptual alternations. This reduction in the alternation rate is highly sensitive to the difference in the stimuli strength. These findings, along with the non monontonic behaviour of mean duration of dominance with varying but equal stimulus strength are in contradiction with the initial propositions of Levelt. They can be altered to accommodate the above as follows.
\\\\
\textit{Modified Levelt's Propositions}
\begin{itemize}
    \item Increasing stimulus strength for one eye will increase the perceptual predominance of that eye's stimulus.

    \item Increasing the difference in stimulus strength between the two eyes will primarily act to increase the average perceptual dominance duration of the stronger stimulus.

    \item Increasing the difference in stimulus strength between the two eyes will reduce the perceptual alternation rate.

    \item  Increasing stimulus strength in both eyes while keeping it equal between eyes will generally increase the perceptual alternation rate, but this effect may reverse at lower stimulus strengths \cite{levelt}.
\end{itemize}
\vspace{2mm}

\section{Future Work}
The study of rivalry may shed light on neural mechanisms underlying perceptual selection and on the resolution of ambiguous sensory information. Another important question is whether the characteristics reviewed above apply to forms of perceptual bistability outside of vision. An affirmative answer would demonstrate even greater generality of the neural properties that are uncovered by bistable perception paradigms. Similar models can then be used to explain the dynamics of CPG \cite{main} (Central Pattern Generators - biological neural circuits that produce rhythmic outputs in the absence of rhythmic input.) This method of study of rivalry can also be extended to test the presence or absence of auditory or olfactory rivalry. \\
An alteration in the binocular rivalry of an individual can be used in diagnosis of mental disorders such as schizophrenia, bipolar disorder and depression. It has been found that patients suffering with bipolar disorder exhibit lower binocular rivalry rate \cite{disord}. Physiological information on binocular rivalry might aid in diagnosing these diseases more efficiently as well as benefit in the design of relevant treatment and drugs.\\
The theoretically important possibility that Proposition IV may not hold true at low stimulus strengths has received only limited experimental verification \cite{levelt}. An advancement in experimental methods and techniques might give better verifiable results. \\
Recent developments show that a single neural structure is not responsible for the perceptual alternations in each of these phenomena. Further probing can be done into models which account for the spatial distribution of neural structures too. Most of the commonly employed testing equipment involve functional MRIs \cite{fmri} and EEGs. Electrocorticography or intracranial EEGs can aid in better isolation of areas within the brain for binocular rivalry.

\vspace{2mm}


\section{Introduction}

Deep neural networks (DNNs) has achieved many state-of-the-art results
in various fields \cite{lecun2015deep}, such as object recognition,
language translation and game-play. A fully understanding of why DNNs
can achieve such good results remains elusive. The often used DNNs
equip much more parameters than the number of the training data. As
Von Neumann said “With four parameters I can fit an elephant, and
with five I can make him wiggle his trunk”. It is no surprise that
such DNNs can well fit the training data. However, counter-intuitive
to the traditional learning theory, such DNNs often do not overfit
(DNNs often generalize well to the test data which are not seen during
the training), which is often referred to as ``apparent paradox''
\cite{zhang2016understanding}.

A series of recent works \cite{xu_training_2018,xu2018understanding,rahaman2018spectral,arpit2017closer,poggio2018theory,perez2018deep,jakubovitz2018generalization,wu2017towards},
both experiments and theories, have gained us more understanding to
this paradox. In this work, we focus on the Fourier analysis of DNNs
\cite{xu_training_2018,xu2018understanding,rahaman2018spectral}.
Although the often used dataset, such as MNIST and CIFAR, are relative
simple compared with practical dataset, the input dimension (the pixel
number of each input image) is still very high for a quantitative
analysis. A good starting point to understand this apparent paradox
is to find an example that is simple enough for analysis but also
preserves this interesting paradox. The simple example turns out to
be the fitting of a function with one-dimension (1-d) input and 1-d
output \cite{xu_training_2018}. A prompt example to understand the
apparent paradox is that a very high-order polynomial fitting for
randomly sampled data points often overfit the training data, that
is, high oscillation occurs around the sample boundary (Runge's phenomenon);
however, a DNN with small weight initialization, no matter how large
size the DNN is, often learns the training data with a relative flat
function \cite{wu2017towards}. Starting from such 1-d functions,
experimentally \cite{xu_training_2018,rahaman2018spectral} and theoretically
\cite{xu2018understanding}, there exists a \emph{Frequency Principle}
(F-Principle) that DNNs often first quickly capture low-frequency
components while keeping high-frequency ones small, and then relatively
slowly captures high-frequency components. By F-Principle, the high-frequency
components of the DNN output is controlled by the training data. High
oscillation which exists in the Runge's phenomenon is then absent
in the DNN fitting. F-Principle also holds well in the often used
dataset \cite{xu_training_2018}, that is, MNIST and CIFAR-10. Theoretical
work indicates that the key ingredient underlying the F-Principle
in general DNN fitting problems is that the power spectrum of the
activation function decays in the Fourier space, where the power-decay
property is easy to be satisfied, such as sigmoid function and rectified
linear unit.

Previous studies focused on the DNN with mean square error \cite{xu_training_2018,xu2018understanding,rahaman2018spectral}.
It is yet to study whether F-Principle applies in the DNN with other
types of loss functions. This is important since loss function varies
in different problems, such as image classification and solving differential
equations \cite{weinan2018deep}. In this work, we perform a theoretical
analysis to show that for a general loss function, e.g., cross entropy,
the F-Principle qualitative holds in the DNN training, which is also
verified by experiments. The first experiment is a classification
problem with the loss function of cross entropy. The second experiment
is to apply DNN to solve Poisson equation by using Dirichlet's principle.

The DNN is a powerful tool to solve differential equations \cite{weinan2018deep,khoo2017solving,weinan2017deep},
especially for high-dimensional problems. It is well-known that different
frequencies converge with different speeds in solving differential
equations by numerical schemes. For example, for the Jacobi method,
low frequency converges much slower than high frequency. Multigrid
method is designed to speed up the convergence, which explicitly first
captures low-frequency parts \cite{briggs2000multigrid}. In addition,
manual frequency marching from low frequency to high frequency has
achieved great success in designing numerical schemes in various problems,
such as inverse scattering problems \cite{bao2015inverse} and Cryo-EM
reconstruction problems \cite{barnett2017rapid}. By showing F-Principle
in solving Poisson's equations, we emphasize that the DNN structure,
which implicitly endows low frequency with high priority, could be
a powerful tool to the problems that benefit from a fast converging
of low frequency. For example, we propose an ideal that combines DNN
and conventional methods (e.g., Jacobi method or Gauss-Seidel method),
in which DNN is in charge of capturing low-frequency parts and conventional
methods are in charge of capturing high-frequency parts. This idea
is exemplified by solving a 1-d Poisson's equation.

\section{F-Principle with general loss function}

Consider a general DNN, and denote its output as $\Upsilon_{\theta}(x)$,
where $\theta$ stands for the DNN parameters and $x$ stands for
the input. Represent $\Upsilon_{\theta}(x)$ with orthonormal basis
$\{p_{k}(x)\}$:

\begin{equation}
\Upsilon_{\theta}(x)=\sum_{k}c_{\theta,k}p_{k}(x),
\end{equation}
where $c_{\theta,k}$ is the coefficient of mode $k$ depending on
$\theta$. Denote the loss at sample $x$ as
\begin{equation}
l_{x}=l(\Upsilon_{\theta}(x)).
\end{equation}
The total loss $L$ is
\begin{equation}
L=\sum_{x}l_{x}=\sum_{x}l(\Upsilon_{\theta}(x)).
\end{equation}
Consider the gradient of the total loss with respect to parameter
$\theta$:
\begin{align}
\frac{\partial L}{\partial\theta} & =\sum_{x}\frac{\partial l(\Upsilon_{\theta}(x))}{\partial\theta}\\
 & =\sum_{x}\frac{\partial l(\Upsilon_{\theta}(x))}{\partial\Upsilon}\frac{\partial\Upsilon_{\theta}(x)}{\partial\theta}\\
 & =\sum_{x}\frac{\partial l(\Upsilon_{\theta}(x))}{\partial\Upsilon}\sum_{k}p_{k}(x)\frac{\partial c_{\theta,k}}{\partial\theta}\\
 & =\sum_{k}\frac{\partial c_{\theta,k}}{\partial\theta}\sum_{x}\frac{\partial l(\Upsilon_{\theta}(x))}{\partial\Upsilon}p_{k}(x).
\end{align}
Let
\begin{equation}
d_{k}\triangleq\sum_{x}\frac{\partial l(\Upsilon_{\theta}(x))}{\partial\Upsilon}p_{k}(x).
\end{equation}
$d_{k}$ is the coefficient of $\partial l(\Upsilon_{\theta}(x))/\partial\Upsilon$
at the component of $p_{k}$. Consider that $\{p_{k}(x)\}$ is Fourier
basis. According to Riemann-Lebesgue lemma, if a function is an integrable
function on an interval, then the Fourier coefficients of this function
tend to 0 as the order $k$ tends to infinity. Therefore, when the
activation function and the target function both are integrable functions
on the considered interval, $c_{\theta,k}$ and $d_{k}$ tend to 0
as the order $k$ tends to infinity. Denote
\begin{equation}
\mathcal{L}_{k}\triangleq\sum_{k}\frac{\partial c_{\theta,k}}{\partial\theta}d_{k}.
\end{equation}
We have
\begin{equation}
\frac{\partial L}{\partial\theta}=\sum_{k}\mathcal{L}_{k}.
\end{equation}
Therefore, we can decompose $\partial L/\partial\theta$ into a summation
of $\mathcal{L}_{k}$, which tends to 0 as the order $k$ tends to
infinity. This analysis implies that for any loss function, the change
of any parameter $\theta$ at each training step is affected more
by lower frequencies, which would rationalize the F-Principle in general
loss functions, as examined in the following experiments.

\section{Experiment: cross entropy loss}

The loss function of cross entropy is widely used in classification
problems. We use experiments to show that F-Principle holds in the
DNN training with this loss function.

\subsection{Toy data}

Consider a target function $y(x)=(y_{1}(x),y_{2}(x))$, where
\begin{equation}
y_{1}(x)=\begin{cases}
1 & x\geq0\\
0 & x<0
\end{cases},
\end{equation}
\begin{equation}
y_{2}(x)=\begin{cases}
0 & x>0\\
1 & x\leq0
\end{cases}.
\end{equation}
This fitting problem is a toy classification problem. In the DNN in
this problem, the output layer has two neurons with softmax as activation
function. The output is denoted as $\Upsilon(x)=(\Upsilon_{1}(x),\Upsilon_{2}(x))$.
The loss function is
\begin{equation}
L=\sum_{j=1}^{2}L_{j}=-\sum_{j=1}^{2}\sum_{x}y_{j}(x)\log\Upsilon_{j}(x)+(1-y_{j}(x))\log(1-\Upsilon_{j}(x)).
\end{equation}

For illustration, we focus on $y_{1}(x)$, which is shown in Fig.\ref{fig:Poisson-2}a.
Next, we examine the convergence of different frequencies. In a finite
interval, the frequency components of a target function are quantified
by Fourier coefficients computed from Discrete Fourier Transform (DFT).
Note that because the frequency in DFT is discrete, we can refer to
a frequency component by its index instead of its physical frequency.
The Fourier coefficient of $y_{1}(x)$ for the $\gamma$-th frequency
component is denoted by $F[y_{1}](\gamma)$ (a complex number in general).
$|F[y_{1}](\gamma)|$ is the corresponding amplitude, where $|\cdot|$
denotes the absolute value. Note that we call $\gamma$ the \emph{frequency
index}. $|F[y_{1}](\gamma)|$ is shown in Fig. \ref{fig:Poisson-2}b.
To examine the convergence behavior of different frequency components
during the training of a DNN, we compute the relative difference of
the DNN output $\Upsilon_{1}(x)$ and $y_{1}(x)$ in frequency domain
at each recording step, i.e.,
\begin{equation}
\Delta_{F}(\gamma)=\frac{|F[\Upsilon_{1}](\gamma)-F[y_{1}](\gamma)|}{|F[\Upsilon_{1}](\gamma)|}.\label{eq:dfreq-1}
\end{equation}

During the training, $\Upsilon_{1}(x)$ captures $y_{1}(x)$ from
low to high frequency in a clear order, as shown in Fig.\ref{fig:Poisson-2}c.
\begin{center}
\begin{figure*}
\begin{centering}
\subfloat[]{\begin{centering}
\includegraphics[scale=0.27]{pic/CELoss/24405/y_train_true}
\par\end{centering}
}\subfloat[]{\begin{centering}
\includegraphics[scale=0.27]{pic/CELoss/24405/ffttrain-impPeak}
\par\end{centering}
}\subfloat[DNN ]{\begin{centering}
\includegraphics[scale=0.27]{pic/CELoss/24405/PeakDifferencePcolor-impPeak}
\par\end{centering}
}
\par\end{centering}
\caption{F-Principle with cross entropy loss. The first output dimension of
the target function is shown in (a) and its Fourier coefficient amplitude
as a function of frequency index is shown in (b). Frequency peaks
are marked by black dots. (c) Relative difference at different recording
steps for different selected frequency indexes. The training data
are evenly sampled in $[-1,1]$ with sample size $201$. We use a
DNN with width 400-400-200-100 with full batch training, the output
layer has two neurons with softmax as activation function, and learning
rate is $2\times10^{-4}$. The parameters of the DNN are initialized
following a Gaussian distribution with mean $0$ and standard deviation
$0.1$.\label{fig:Poisson-2} }
\end{figure*}
\par\end{center}

\subsection{MNIST data}

To verify that the F-Principle holds in the image classification problems
(MNIST) with the loss function of cross entropy, we perform Fourier
analysis in the first principle component of the input space. The
procedure is as follows.

The training set is a list of images with labels: $\{\vec{x}_{k},\vec{y}_{k}\}_{k=0}^{n-1}$.
Each image is represented by a vector $\vec{x}_{k}\in\mathbb{R}^{N_{in}}$,
where $N_{in}=784$ is the pixel number of an image. $\vec{y}_{k}$
is an one-hot vector indicating the label. The dimensions for the
input layer and the output layer\emph{ }are $N_{in}$ and $10$, respectively.
First, we compute the first principle direction. Transform each image
by
\begin{equation}
\vec{x}_{j}^{\prime}=\vec{x}_{j}-\frac{1}{n}\sum_{k=0}^{n-1}\vec{x}_{k},\quad j=0,1,\cdots,n-1.
\end{equation}
Denote all images by $X=[\vec{x}_{0}^{\prime},\vec{x}_{1}^{\prime},\cdots,\vec{x}_{n-1}^{\prime}]\in\mathbb{R}^{N_{in}\times n}$.
The covariance matrix is $C_{x}=XX^{T}$. The eigenvector of the maximal
eigenvalue of $C_{x}$ can be obtained, denoted by $\vec{p}_{1}\in\mathbb{R}^{N_{in}}$,
i.e., the first principle direction. The projection of each image
in the $\vec{p}_{1}$ direction is $x_{k}^{\prime}\triangleq\vec{p}_{1}^{T}\vec{x}_{k}$.
We rescale $x_{k}^{\prime}$ to $x_{k}$ such that $x_{k}\in[0,1]$
by
\begin{equation}
x_{k}=\frac{(x_{k}^{\prime}-\min_{j}x_{j}^{\prime})}{\max_{l}(x_{l}^{\prime}-\min_{j}x_{j}^{\prime})}.
\end{equation}
Then, the sample set is $S=\{(x_{0},\vec{y}_{0}),(x_{1},\vec{y}_{1}),\cdots,(x_{n-1},\vec{y}_{n-1})$.
For illustration, we only consider the first component of $\vec{y}$,
i.e., $\vec{y}^{(1)}$. Note that now $\{x_{k}\}_{k=0}^{n-1}$ is
a non-uniform sampling. Using non-uniform FFT (NUFFT), we can obtain
\begin{equation}
F[\vec{y}^{(1)}][\gamma_{k}]=\sum_{j=0}^{n-1}\vec{y}_{k}^{(1)}\exp\left(-2\pi ix_{j}k\right).
\end{equation}
Then, the sampling on the Fourier domain is
\begin{equation}
S_{\gamma}=\{(\gamma_{0},F[\vec{y}^{(1)}](\gamma_{0})),(\omega_{1},F[\vec{y}^{(1)}](\gamma_{1})),\cdots,(\gamma_{n-1},F[\vec{y}^{(1)}](\gamma_{n-1}))\}.
\end{equation}
After each training step, we feed each $\vec{x}_{k}$ into the DNN
and obtain the DNN output $\Upsilon(\vec{x}_{k})$:
\begin{equation}
S_{T}=\{(x_{0},\Upsilon(\vec{x}_{0})),(x_{1},\Upsilon(\vec{x}_{1})),\cdots,(x_{n-1},\Upsilon(\vec{x}_{n-1}))\}.
\end{equation}
Using non-uniform FFT (NUFFT), similarly, we can obtain the sampling
of the DNN's first dimension output on the Fourier domain
\begin{equation}
S_{\Upsilon,\gamma}=\{(\gamma_{0},F[\Upsilon^{(1)}](\gamma_{0})),(\omega_{1},F[\Upsilon^{(1)}](\gamma_{1})),\cdots,(\gamma_{n-1},F[\Upsilon^{(1)}](\gamma_{n-1}))\},
\end{equation}
which is shown in Fig.\ref{fig:Poisson-2-1}a. We then examine the
relative error of certain selected important frequency components
(marked by black squares). As shown in the first column in Fig.\ref{fig:Poisson-2-1}b,
we can observe that the DNN tends to capture low-frequency components
first.
\begin{center}
\begin{figure*}
\begin{centering}
\subfloat[]{\begin{centering}
\includegraphics[scale=0.27]{pic/CELoss/MNISTDNN/DNNmpca_1233732/fftTrainCleanPeak}
\par\end{centering}
}\subfloat[]{\begin{centering}
\includegraphics[scale=0.27]{pic/CELoss/MNISTDNN/DNNmpca_1233732/PeakDifferencePcolorPCA_pca_1-impPeak}
\par\end{centering}
}
\par\end{centering}
\caption{F-Principle with cross entropy loss on MNIST dataset. The Fourier
coefficient amplitude of the first output dimension of the target
function is shown in (a). Frequency peaks are marked by black dots.
(b) Relative difference at different recording steps for different
selected frequency indexes. The training data are $10000$ test samples
of MNIST dataset. We use a DNN with width 400-200 with batch size
as 128, the output layer has 10 neurons with softmax as activation
function, and learning rate is $10^{-5}$. The parameters of the DNN
are initialized following a Gaussian distribution with mean $0$ and
standard deviation $0.2$. \label{fig:Poisson-2-1} }
\end{figure*}
\par\end{center}

\section{Experiment: Poisson's equations}

Consider one-dimension (1-d) Poisson's equation \cite{weinan2018deep,evans2010partial}:
\begin{equation}
-\triangle u(x)=g(x),\quad x\in\Omega=(-1,1)\label{eq:Poisson1}
\end{equation}
\[
u(x)=0,\quad x=-1,1.
\]
The Poisson's equation can be solved by numerical schemes (e.g., Jacobi
method) or DNN. As well known, high frequency converges faster in
the Jacobi method. In the following, we would show that high frequency
converges slower when the DNN is applied to solve the above Poisson's
equation.

\subsection{Central differencing scheme and Jacobi method}

$[-1,1]$ is uniformly discretized into $n+1$ points with step $\Delta x=2/n$,
i.e., $x_{0},x_{1},\cdots,x_{n}$. The Poisson's equation in Eq. (\ref{eq:Poisson1})
can be solved by central differencing scheme:
\begin{equation}
-\triangle u_{i}=-\frac{u_{i+1}-2u_{i}+u_{i-1}}{(\triangle x)^{2}}=g(x_{i}),\quad i=1,2,\cdots,n.
\end{equation}
Write the above in the matrix form:
\begin{equation}
Au=g,\label{eq:auf}
\end{equation}
where
\begin{equation}
A=\left(\begin{array}{cccccc}
2 & -1 & 0 & 0 & \cdots & 0\\
-1 & 2 & -1 & 0 & \cdots & 0\\
0 & -1 & 2 & -1 & \cdots & 0\\
\vdots & \vdots & \cdots &  &  & \vdots\\
0 & 0 & \cdots & 0 & -1 & 2
\end{array}\right)_{(n-1)\times(n-1)},
\end{equation}

\begin{equation}
u=\left(\begin{array}{c}
u_{1}\\
u_{2}\\
\vdots\\
u_{n-2}\\
u_{n-1}
\end{array}\right),\quad g=(\triangle x)^{2}\left(\begin{array}{c}
g_{1}\\
g_{2}\\
\vdots\\
g_{n-2}\\
g_{n-1}
\end{array}\right),\quad x_{i}=2\frac{i}{n}.
\end{equation}
If $n$ is not a large number, Eq. (\ref{eq:auf}) can be solved by
performing the inverse of $A$. When $n$ is a very large number,
this problem can be solved by iterative schemes. For example, we illustrate
the Jacobi method. Let $A=D-L-U$, where $D$ is diagonal, and $L$
and $U$ are the strictly lower and upper parts of $-A$, respectively.
Then, we can obtain
\begin{equation}
u=D^{-1}(L+U)u+D^{-1}g.
\end{equation}
The Jacobi iteration is
\begin{equation}
u^{l+1}=D^{-1}(L+U)u^{l}+D^{-1}g.
\end{equation}
We perform error analysis of the above iteration process. Denote $u^{*}$
as the true value obtained by directly performing inverse of $A$
in Eq. (\ref{eq:auf}). The error at step $l+1$ is $e^{l+1}=u^{l+1}-u^{*}$.
Then, $e^{l+1}=R_{J}e^{l}$, where $R_{J}=D^{-1}(L+U)$. The converging
speed of $e^{l}$ is determined by the eigenvalues of $R_{J}$, that
is,
\begin{equation}
\lambda_{k}=\lambda_{k}(R_{J})=\cos\frac{k\pi}{n},\quad k=1,2,\cdots,n-1,
\end{equation}
and the corresponding eigenvector $v_{k}$ is
\begin{equation}
v_{k,j}=\sin\frac{jk\pi}{n},j=1,2,\cdots,n-1.
\end{equation}
Write
\begin{equation}
e^{l}=\sum_{k=1}^{n-1}\alpha_{k}^{l}v_{k},
\end{equation}
where $\alpha_{k}^{l}$ can be understood as the magnitude of $e^{l}$
in the direction of $v_{k}$. Then,
\begin{equation}
e^{l+1}=\sum_{k=1}^{n-1}\alpha_{k}^{l}R_{J}v_{k}=\sum_{k=1}^{n-1}\alpha_{k}^{l}\lambda_{k}v_{k}.
\end{equation}
\[
\alpha_{k}^{l+1}=\lambda_{k}\alpha_{k}^{l}.
\]
Therefore, the converging speed of $e^{l}$ in the direction of $v_{k}$
is controlled by $\lambda_{k}$. Since
\begin{equation}
\cos\frac{k\pi}{n}=-\cos\frac{(n-k)\pi}{n},
\end{equation}
the frequencies $k$ and $(n-k)$ are closely related and converge
with the same speed. Consider the frequency $k<n/2$, $\lambda_{k}$
is larger for lower frequency. Therefor, lower frequency converges
slower in the Jacobi method.

\subsection{DNN approach}

Similar as the loss function in Ref \cite{weinan2018deep}, we consider
the following loss function (energy method)
\begin{equation}
I(u)=\int_{\Omega}\left(\frac{1}{2}|\nabla_{x}u(x)|^{2}-g(x)u(x)\right)dx+\beta\int_{\partial\Omega}u(x)^{2}ds.\label{eq:Energy}
\end{equation}
It is equivalent to solve Poisson's equation by finding the function
that minimizes $I(u)$ (Dirichlet’s principle) \cite{evans2010partial}.
The last term in $I(u)$ is a penalty in order to satisfy the boundary
condition. The DNN structure is 1-d input (i.e., $x$) and 1-d output
(denoted as $\Upsilon(x)$) for solving Eq. (\ref{eq:Poisson1}).
$\beta$ is a constant.

The procedure is similar. We discretized $[-1,1]$ into $n+1$ even-space
points. In each training step, we compute $\Upsilon(x_{i})$ for $i=0,1,2,\cdots n$.
The gradient of $I(\Upsilon)$ with respect to parameter $\Theta$
is
\begin{equation}
\frac{dI(\Upsilon)}{d\Theta}=\sum_{i=0}^{n}\frac{dI(\Upsilon(x_{i}))}{d\Theta}.
\end{equation}
At each training step, we compare $\Upsilon(x)$ and $u^{*}(x)$ in
the Fourier domain. Note that $u^{*}(x)$ is the one obtained by directly
performing inverse of $A$ in Eq. (\ref{eq:auf}).

\subsection{Experiment}

Consider
\begin{equation}
g(x)=\sin(x)+4\sin(4x)-8\sin(8x)+16\sin(24x).
\end{equation}
As shown in Fig. \ref{fig:Poisson}a, after training, the DNN output
can well fit the solution $u^{*}$ obtained by directly performing
inverse of $A$ in Eq. (\ref{eq:auf}). As shown in Fig. \ref{fig:Poisson}b,
there are three peaks of $u^{*}$ in the Fourier domain.

To examine the convergence behavior of different frequency components
during the DNN training, we compute the relative difference of the
DNN output $\Upsilon(x)$ and $u^{*}(x)$ in frequency domain at each
recording step, i.e.,
\begin{equation}
\Delta_{F}(\gamma)=\frac{|F[u^{*}](\gamma)-F[\Upsilon](\gamma)|}{|F[u^{*}](\gamma)|}.\label{eq:dfreq}
\end{equation}
As shown in Fig. \ref{fig:Poisson}c, F-Principle holds well in solving
Poisson's equation \cite{xu_training_2018,xu2018understanding}. For
comparison, we also show that low frequency converges much slower
than high frequency in Jacobi method, as shown in Fig. \ref{fig:Poisson}d.
\begin{center}
\begin{figure}
\begin{centering}
\subfloat[]{\begin{centering}
\includegraphics[scale=0.27]{pic/Poisson_1_14128/y}
\par\end{centering}
}\subfloat[]{\begin{centering}
\includegraphics[scale=0.27]{pic/Poisson_1_14128/ffttrain-impPeak}
\par\end{centering}
}
\par\end{centering}
\begin{centering}
\subfloat[DNN ]{\begin{centering}
\includegraphics[scale=0.27]{pic/Poisson_1_14128/PeakDifferencePcolor-impPeak}
\par\end{centering}
}\subfloat[ Jacobi]{\begin{centering}
\includegraphics[scale=0.27]{pic/JacoPoisson_1_55911/PeakDifferencePcolor-impPeak}
\par\end{centering}
}
\par\end{centering}
\caption{Frequency domain analysis of the Poisson's equation in Eq. (\ref{eq:Poisson1})
with $g(x)=\sin(x)+4\sin(4x)-8\sin(8x)+16\sin(24x)$. (a) $u(x)$.
The true value is computed by central differencing scheme with directly
compute the inverse of coefficient matrix. (b) $|F[u^{*}]|$ (red
solid line) as a function of frequency index. Frequency peaks are
marked by black dots. (c, d) Relative difference at different recording
steps for different selected frequency indexes. The training data
and the test data are evenly sampled in $[-1,1]$ with sample size
$51$ and $401$, respectively. We use a DNN with width 4000-800 with
full batch training. The learning rate is $5\times10^{-6}$ at beginning
and halved every $10000$ training epochs. $\beta$ is $10$. Each
step consists of four epochs. (d) is the result of Jacobi iteration.
The parameters of the DNN are initialized following a Gaussian distribution
with mean $0$ and standard deviation $0.05$. \label{fig:Poisson} }
\end{figure}
\par\end{center}

\section{Combination of DNN and convention methods}

In light of the above numerical simulations, it is natural to consider
if we can combine DNN and Jacobi method to solve the Poisson's equation.
For simplicity, we call the combination method D-Jacobi method.

In the first part of D-Jacobi method, we solve the Poisson's equation
by DNN with $M$ steps. In the second part, we use the DNN output
at step $M$ as the initial value for the Jacobi method.

We solve the problem in Fig. \ref{fig:Poisson} by a laptop (Dell,
Precision 5510). As shown in Fig.\ref{fig:Poisson-1}a, the DNN loss
fluctuates after some running time. We use Jacobi method to solve
the problem after some time points, which are indicated by vertical
dashed lines. As shown in Fig.\ref{fig:Poisson-1}b (Fig.\ref{fig:Poisson-1}c),
green stars indicates the $|\Upsilon-u^{*}|_{\infty}$ of the DNN
output at different steps. Dashed lines indicates the evolution of
the Jacobi (Gauss-Seidel) method. As we can see, if the selected timing
is too early, it would still take long time to converge to a small
error, because the low frequencies are not converged, yet. If the
selected timing is too late, much time would be waste because the
DNN is hard to capture high frequencies and fluctuates a lot. The
selected timing of the green or the red one is a better choice. In
practice, a better way to select the timing is when the loss gets
flat and fluctuated for a short while.
\begin{center}
\begin{figure*}
\begin{centering}
\subfloat[Loss]{\begin{centering}
\includegraphics[scale=0.3]{pic/Combine/GS11828/Loss_ver}
\par\end{centering}
}\subfloat[Jacobi: $|u-u^{*}|_{\infty}$]{\begin{centering}
\includegraphics[scale=0.3]{pic/Combine/1_11828/ydiff}
\par\end{centering}
}\subfloat[GS: $|u-u^{*}|_{\infty}$]{\begin{centering}
\includegraphics[scale=0.3]{pic/Combine/GS11828/ydiff}
\par\end{centering}
}
\par\end{centering}
\caption{Combination methods for solving the Poisson's equation in Eq. (\ref{eq:Poisson1})
with $g(x)=\sin(x)+4\sin(4x)-8\sin(8x)+16\sin(24x)$. The abscissa
is the computer running time. (a) Loss is the form in Eq. (\ref{eq:Energy}).
We use Jacobi method to solve the problem after several time points,
which are indicated by vertical dashed lines. (b) Green stars indicates
the $|\Upsilon-u^{*}|_{\infty}$ at different steps. Dashed lines
indicates the evolution of the Jacobi method. (c) Gauss-Seidel method.
The training data are evenly sampled in $[-1,1]$ with sample size
$1001$. $u^{*}$ is computed by central differencing scheme with
directly computing the inverse of coefficient matrix. We use a DNN
with width 4000-500-400 with full batch training, and learning rate
is $5\times10^{-4}$. $\beta$ is $10$. The parameters of the DNN
are initialized following a Gaussian distribution with mean $0$ and
standard deviation $0.02$.\label{fig:Poisson-1} }
\end{figure*}
\par\end{center}

\section{Discussion}

In this work, we have shown that F-Principle holds well in the DNN
training with a general loss function, extending the study of F-Principle
in the loss function of mean square error in previous works \cite{xu_training_2018,xu2018understanding,rahaman2018spectral}.
Along with the previous study that F-Principle holds in both DNN and
convolutional neural networks with the activation function of either
tanh or Relu \cite{xu_training_2018,xu2018understanding}, these works
implicate that the F-Principle may provide understandings to the generalization
ability of general DNNs.

We also show that the generality of F-Principle in the DNN training
could potentially be useful in designing algorithms for solving practical
problems. To be specific, we apply DNN to solve 1-d Poisson's equation.
Compared with conventional numerical schemes, DNN could potentially
work better in rather high dimensions \cite{weinan2018deep}. In addition,
it does not requires discretization for the DNN method, which would
be much easier to be implemented. In future, it would be interested
to use DNN's F-Principle to develop numerical schemes for solving
various problems which would benefit from a fast converging of low
frequency.

\section*{Acknowledgments}

The author wants to thank Weinan E, Wei Cai for helpful discussion.
The author also wants to thank Tao Luo, Zheng Ma, Yanyang Xiao and
Yaoyu Zhang for the discussion of the F-Principle. This work was funded
by the NYU Abu Dhabi Institute G1301.

\subsection*{ \bibliographystyle{unsrt}