  In this paper we show, that the $J/\psi$ production in DIS, is the main
source of information about the events with two parton shower production. We
attempt to develop our theoretical acumen of this process, to a level
compatible with the theoretical description of inclusive DIS. We revisit the
problem of the linear evolution equation for the double gluon densities, and
include Bose-Einstein enhancement to these equations. We find that the
Bose-Einstein correlations lead to an increase of the anomalous dimension,
which turns out to be suppressed as $1/(N^2_c -1)^2$, in agreement with the
estimates for the twist four anomalous dimension. We believe that understanding
what happens to these contributions at ultra high energies, is a key question
for an effective theory, based on high energy QCD. We derive the evolution
equation for the scattering amplitude of two dipoles with a nucleus, taking
into account the shadowing corrections, and investigate the analytical
solutions in two distinct kinematic regions: deep in the saturation region, and
in the vicinity of the saturation scale. The suggested non-linear evolution
equation is a direct generalization of the Balitsky-Kovchegov equation, which
has to be solved with the initial condition that depends on the saturation
scale $Q_s(Y=Y_0,b)$. With the goal of finding a new small parameter, it is
instructive to compare the solution of the non-linear equation with the
qusi-classical approximation, in which in the initial condition we replace
$Q_s(Y=Y_0,b)$ by $Q_s(Y,b)$. Our final result is that the shadowing
corrections in the elastic amplitude generate the survival probability, which
suppresses the growth of the amplitude with energy, caused by the Bose-Einstein
enhancement.
