\section{Introduction}

The idea that physical processes can be considered as computations done on
some universal machines traces back to Turing and von
Neumann~\cite{Coveney:1996tx},
and the growth of the computational power allowed for further development of these concepts.
This resulted in a completely new approach to science in which the complexity of observed phenomena is closely related to the complexity of computational resources needed to simulate them~\cite{Wolfram:1985je}.
In addition, there are physical phenomena that simply cannot be traced with analytical tools,
which further motivated a computational approach to physics~\cite{Moore:1990bn}.
Moreover, the idea of quantum computation~\cite{Feynman:82} lead to a
discovery of a few problems that seem not
efficiently traceable on classical computers but efficiently
on a quantum version~\cite{Shor:1994tb,Grover:1996}.

Classical physics can be simulated on universal Turing machines, or other
computationally equivalent models~\cite{Wolfram:2001ue}.  On the other hand,
efficient simulation of quantum systems requires a replacement of
deterministic universal Turing machines with
quantum computers whose states are non-classically correlated. Such machines
can even simulate any local quantum system efficiently~\cite{Deutsch:1985ba,Lloyd:96}.
Can we experimentally distinguish between these two descriptions of the universe
using a logically self-contained computational approach?

In this paper, we show that
there are processes which cannot be simulated on local classical machines at all,
independently of the available classical resources.
We first introduce the notion of Kolmogorov complexity,
a measure of the classical complexity of a phenomena,
and later apply it to derive a bound on classical descriptions~\cite{Li:2004kz}. Next, we use the fact that Kolmogorov complexity can be approximated by compression algorithms~\cite{Cilibrasi:2005jna}.

We compress
experimental data obtained from polarisation
measurements on entangled photon pairs
and show the violation of a classical bound.

Let's consider the description of a machine, whether classical or quantum, that outputs a string $x$ made of 0's and 1's.
In the case of a Turing machine $U$, we can always write a program $\Lambda$ that generates $x$.
The simplest such program is obviously `PRINT $x$'.
However, this is not optimal: in many cases the program can be much shorter than the string itself.

This brings us to the concept of Kolmogorov complexity $K(x)$, the minimal
length $l(\Lambda)$ of all programs $\Lambda$ that reproduce a specific output
$x$. If $K(x)$ is of the order of
the length of the output
$l(x)$ then our algorithmic description of $x$ is inefficient, and $x$ is called algorithmically random~\cite{Li:2009wu}.
In most
cases $K(x)$ cannot be computed~\cite{Cover:2006ub}.
To circumvent this issue, we can estimate $K(x)$ with some efficient lossless compression algorithm $C(x)$~\cite{Cilibrasi:2005jna}.

We now extend this picture by considering two Turing machines
$U_A$ (Alice) and $U_B$ (Bob), which are spatially separated.
If these machines cannot communicate,
they generate two output strings that are independent,
although the programs fed into the machines can be correlated.
Moreover, the input programs are classical bit strings so the correlations
between them must be classical.

We determine the complexity of the generated strings using the
\emph{Normalized Information Distance} (NID)~\cite{Li:2004kz}.
This distance allows for a comparison of two data sets without detailed knowledge about their origin.
In practice, we evaluate an approximation to the NID, the \emph{Normalized Compression Distance} (NCD)~\cite{Cilibrasi:2005jna}, using a lossless compression software, in our case the LZMA Utilities, based on the Lempel-Ziv-Markov chain algorithm~\cite{ref:lzma}.

\section{Simulation by deterministic universal Turing machines}
\begin{figure}
    \begin{center}
        \includegraphics[width=.8\columnwidth]{figure1.eps}
        \caption{Measurement: $N$ particles enter a measuring device,
          characterized by two possible polarizer settings $a_0$ and $a_1$,
          which generates a bit string of $N$ measurement outcomes.
          Simulation: a universal Turing machine is fed with a
          program elements $\Lambda_i$ and information about the chosen setting
          $a_0$ or $a_1$. It delivers an output string of length $N$.}
        \label{fig0}
    \end{center}
\end{figure}
We
consider a model experiment, similar to the one used for testing the Bell
inequalities~\cite{Bell}: a source emits pairs of photons that travel to two
separate polarization analyzers M$_A$ (Alice) and M$_B$ (Bob).
Each analyzer has two outputs associated with bit values 0 and 1,
and can be
set along directions $a_0$ or $a_1$ for M$_A$, and $b_0$ or $b_1$ for M$_B$.
The record of the outputs from each analyzer forms a bit string (see Fig.~\ref{fig0}).

The output $x$ of
each individual analyzer can be described as the output of a Turing machine $U$, fed with the settings $a_j$ or $b_k$, and a program $\Lambda$.
The program will contain the information for
generating the correct output for
every detection event and for every setting.

If we consider a string of finite length $l(x)=N$, $\Lambda$ will have to describe the $4^N$ possible events.
The length of the shortest $\Lambda$ is equal to the Kolmogorov complexity of the generated string.

Next, we consider the simulation of the experiment with two local non-communicating machines $U_A$ and $U_B$ (see Fig.~\ref{fig:turing_pair}).
We feed a program $\Lambda$ to
both of them and obtain two output strings, $x$ and $y$, both of length $N$.
In this case, the program has to describe the behavior of all $2N$
events
for all possible settings $a_j$ and $b_k$, hence $16^N$ possible events.
\begin{figure}
    \begin{center}
        \includegraphics[width=.8\columnwidth]{figure2.eps}
        \caption{Local classical machines simulating the generation of strings $x$ and $y$ by correlation measurements on an entangled state.}
        \label{fig:turing_pair}
    \end{center}
\end{figure}
\subsection{Normalized Information Distance}
The Kolmogorov complexity of two bit strings $K(x, y)$ is the length of the shortest program generating them simultaneously. $K(x, y)$ can be shorter than $K(x)\,+\,K(y)$ if $x$ and $y$ are correlated - the more correlated they are, the simpler it is to compute one string knowing the other. This idea was further carried out by Cilibrasi and Vitanyi~\cite{Cilibrasi:2005jna} who constructed a distance measure between $x$ and $y$ called \emph{Normalized Information Distance} (NID),
\begin{equation}
    \text{NID}(x, y)\,=\,\frac{K(x,y)\,-\,\min\{K(x), K(y)\}}{\max\{K(x), K(y)\}}\,.
\end{equation}

The NID obeys all required properties of a metric, in particular, the triangle inequality
\begin{equation}\label{eq:triangle}
    \text{NID}(x, y)\,+\,\text{NID}(y, z)\,\geq\,\text{NID}(x, z)\,.
\end{equation}
The above inequality holds up to a correction of order $\log(l(x))$, which can
be neglected for sufficiently long strings~\cite{Cilibrasi:2005jna}.

\subsection{Information Inequality}
We consider the bit strings $x_{a_j}$ and $y_{b_k}$ generated by Alice and Bob
with fixed setting $a_j$ and $b_k$.
Equation~(\ref{eq:triangle}) then transforms into
\begin{equation}\label{eq:NID_a0b1}
    \text{NID}(x_{a_0}, y_{b_0})+\text{NID}(y_{b_0}, y_{b_1})\geq\text{NID}(x_{a_0}, y_{b_1})\,.
\end{equation}
However, $\text{NID}(y_{b_0}, y_{b_1})$ cannot be determined experimentally because the strings $y_{b_0}$ and $ y_{b_1}$ come from measurements of incompatible observables.
We therefore use the triangle inequality
\begin{equation}\label{eq:NID_b0b1}
    \text{NID}(x_{a_1}, y_{b_0})+\text{NID}(x_{a_1}, y_{b_1})\geq\text{NID}(y_{b_0}, y_{b_1})\,,
\end{equation}
and combine it with inequality~(\ref{eq:NID_a0b1})
to obtain a quadrangle inequality:
\begin{eqnarray}\label{eq:ineq}
    \text{NID}(x_{a_0}, y_{b_0})+\text{NID}(x_{a_1}, y_{b_0}) + &\text{NID}(x_{a_1}, y_{b_1})\geq \nonumber\\
    &\text{NID}(x_{a_0}, y_{b_1}).
\end{eqnarray}

Similar to various tests of Bell inequalities, we introduce a scalar quantity
$S'$ that quantifies the degree of violation of Eq.~(\ref{eq:ineq}):
\begin{eqnarray}\label{eq:ineq2}
   S' &=& \text{NID}(x_{a_0}, y_{b_1}) -  \text{NID}(x_{a_0},y_{b_0}) \nonumber \\
   &-& \text{NID}(x_{a_1},y_{b_0})\,-\,\text{NID}(x_{a_1},y_{b_1}) \leq 0
\end{eqnarray}

In order to experimentally test this inequality,
we have to address the following problem.
We can set up a source to generate entangled photon pairs in a state of our
choosing, but
we cannot control the nature of the measurement.
For every experimental run $i$ with the same preparation the resulting string $x_{i,a_j}$ can be different.
Consequently, the corresponding program $\Lambda_i$ is different for every experimental run.

It is reasonable to assume that for every two experimental runs $i$ and $i'$ the complexity of the generated strings  remains the same:
$K(x_{i,a_j})\,=\,K(x_{i',a_j})$ and $K(x_{i,a_j}, y_{i,b_k})\,=\,K(x_{i',a_j}, y_{i',b_k})$.
Without these assumptions
the same physical preparation of the experiment has different consequences and thus the notion of preparation loses its meaning. More generally, the predictive power of science can be expressed by saying that the same preparation results in the same complexity of observed phenomena.

\section{Estimation of Kolmogorov complexity}

In general the Kolmogorov complexity cannot be evaluated, but it can be estimated. One can adapt two conceptually different approaches.

\subsection{Statistical Approach}

This
approach takes into account the ensemble of all possible $N$-bit strings and asks
about their average Kolmogorov complexity. It can be shown that this average
equals  the Shannon entropy $H(X)$ of the ensemble~\cite{Cover:2006ub}, and thus
\begin{equation}\label{eq:nid_entropy}
    \langle \text{NID}(x, y)\rangle=\frac{H(x, y)-\min\{H(x), H(y)\}}{\max\{H(x), H(y)\}}\,.
\end{equation}

Inequality (\ref{eq:ineq})
becomes a type of entropic Bell inequality introduced by Braunstein and
Caves~\cite{Braunstein:1988en} if local entropies are maximal, i.e.,
$H(x)=H(y)=N$.
They showed that for a maximally entangled polarization state of two photons,
and polarizer angles obeying the constraints
\begin{eqnarray}
 \vec{a_0}\cdot\vec{b_1}&=&\cos 3\theta\,,\nonumber \\
 \vec{a_0}\cdot\vec{b_0}&=&\vec{a_1}\cdot\vec{b_0}=\vec{a_1}\cdot\vec{b_1}=\cos
 \theta\,,\label{eq:braunstein}
\end{eqnarray}
inequality~(\ref{eq:ineq}) is violated for an appropriate range of $\theta$.
Calculating the entropy $H(x, y)$ using
the probability distributions predicted by quantum mechanics,
it is possible to obtain the expected value of $S'$ as a function of $\theta$
(Fig.~\ref{fig:violation_versus_separation}a).
The maximal violation of this inequality is $S'=0.24$, with a separation of
$\theta\,=\,8.6^\circ$.

\subsection{Algorithmic approach}
It is possible
to avoid a statistical description of our experiment
following the ideas pioneered in~\cite{Cilibrasi:2005jna}. There, it was shown that the Kolmogorov complexity can be well approximated by
the application of
compression algorithms.
This approximation introduces the new distance called \emph{Normalized Compression Distance}~(NCD)
\begin{equation}\label{eq:NCD}
    \text{NCD}(x,y)=\frac{C(x,y)-\min\{C(x),C(y)\}}{\max\{C(x),C(y)\}}\,,
\end{equation}
where $C(x)$ is the length of the compressed string $x$, and $C(x,y)$ is the length of the compressed concatenated strings $x,y$.
Replacing NID with NCD in Eq.~(\ref{eq:ineq2}) leads to a new inequality:
\begin{eqnarray}\label{eq:ineq3}
   S'\rightarrow S&=&\text{NCD}(x_{a_0},y_{b_1}) - \text{NCD}(x_{a_0},y_{b_0}) \nonumber
   \\ &-&\text{NCD}(x_{a_1},y_{b_0}) - \text{NCD}(x_{a_1},y_{b_0}) \leq 0\,.
\end{eqnarray}
This expression can be tested experimentally
because the NCD distance measure is operationally defined.

\section{Choice of compressor}

Before moving to the experiment, we need to ensure the suitability of the
compression software we use to evaluate the NCD. For this, we
numerically simulate the outcome of an experiment, based on a distribution of
results predicted by quantum physics.
Among the packages we tested, we found that the LZMA Utility~\cite{ref:lzma}
approaches the Shannon limit~\cite{Shannon:1948iy} most closely.

The simulation also allows us to verify the angle that maximizes the
violation of Eq.~(\ref{eq:ineq3}) predicted from Eq.~(\ref{eq:braunstein}).
The results of the simulation are presented in Fig.~\ref{fig:violation_versus_separation}.
More details on the generation of the simulated data and the choice of the compressor are provided in the Appendix.

\section{Experiment}

In our experiment (see Fig.~\ref{fig:setup}), the output of a
grating-stabilized laser diode (LD, central wavelength 405\,nm) passes
through a single mode optical fiber (SMF) for spatial mode filtering, and is
focused to a beam waist of 80\,${\mu}$m into a 2\,mm thick BBO crystal. In
this crystal (cut for type-II phase-matching), photon pairs are generated via
spontaneous parametric down-conversion (SPDC) in a slightly non-collinear
configuration.  A half-wave plate ($\lambda$/2) and a pair of compensation
crystals (CC) take care of the temporal and transversal
walk-off~\cite{kwiat:95}.
\begin{figure}
    \begin{center}
        \includegraphics[width = 0.9\columnwidth]{figure3.eps}
        \caption{Schematic of the experimental set-up.
        Polarization correlations of entangled-photon pairs are measured by the polarization analyzers
        M$_A$ and M$_B$,
        each consisting of a half wave plate ($\lambda/2$) followed by a polarization beam splitter (PBS).
        All photons are detected by Avalanche photodetectors D$_H$ and D$_V$, and registered in a coincidence unit (CU).}
        \label{fig:setup}
    \end{center}
\end{figure}
Two spatial modes (labeled A and B) of down-converted light, defined
by the SMFs for 810\,nm, are matched to the pump mode to optimize
the
collection~\cite{kurtsiefer:01}. In type-II SPDC, each down-converted pair
consists of an ordinary and extraordinarily polarized photon, corresponding to
horizontal (H) and vertical (V) in our setup.  A pair of
polarization controllers (PC) ensures that the SMFs do not affect the
polarization of the collected photons. To arrive at an approximate singlet
Bell state, the phase $45^\circ$ between the two decay possibilities in the polarization state
\begin{equation}
    \ket{\psi}={1\over\sqrt{2}}\left(\ket{H}_A\ket{V}_B+e^{i45^\circ}\ket{V}\ket{H}_B\right),
    \label{eq:gen_ouput_state}
\end{equation}
is adjusted to $45^\circ$\,=\,$\pi$ by tilting the CC.

In the polarization analyzers (Fig.~\ref{fig:setup}), the photons from SPDC
are projected onto arbitrary linear polarization by $\lambda$/2 plates, set
to half of the analyzing angles $\theta_{A(B)}$, and
polarization beam splitter (PBS) in each analyzer.
Photons are detected by avalanche photo diodes (APDs), and corresponding detection events from the same pair identified by a coincidence unit (CU) if they arrive within $\approx\pm$3\,ns of each other.

The quality of polarization entanglement is tested by probing the polarization
correlations in a basis complementary to the intrinsic HV basis of the
crystal; for Bell states $\ket{\psi^\pm}$, strong polarization correlations
are e.g. expected in a $\pm$\,45$^\circ$ linear polarization basis.

With interference filters (IF) of 5\,nm bandwidth (FWHM) centered at 810\,nm,
we observe a visibility $V_{45}$ = 99.9$\pm$0.1\%. The visibility in
the natural H/V basis of the type-II down-conversion process also reaches
$V_{\rm HV}$ = 99.9$\pm$0.1\%. A separate test of a CHSH-type Bell inequality~\cite{clauser:69}
leads to a value of $S=2.826\pm0.0015$. This indicates a relatively high quality
of polarization entanglement; the uncertainties in the visibilities are obtained from propagated Poissonian counting statistics.

\subsection{Measurement and Data Post-processing}

We record two-fold coincidences of detection events between
detectors at A and B.
For each PBS, the transmitted output is associated with 0 and the reflected one with 1.
The
resulting binary strings $x$ from A, and $y$ from B are written into two
individual binary files.
From these, we calculate the NCD using Eq.~(\ref{eq:NCD}).
This procedure is repeated for each of the four settings ($a_0,b_0$), ($a_1,b_0$), ($a_1,b_1$), and ($a_0,b_1$) in order to obtain the value for $S$.

To remove the bias due to differences in the detection efficiency of the APDs
in the experiment, we also measure for each setting the associated orthogonal
ones (see Appendix for details).

\begin{figure}
    \begin{center}
        \includegraphics[width = 0.95\columnwidth]{figure4.eps}
        \caption{\label{fig:violation_versus_separation}
        Plots of $S$ versus angle of separation $\theta$.
        (a) Result obtained from Eq.~(\ref{eq:nid_entropy})
        (b) result obtained from using the LZMA compressor on a simulated data ensemble,
        (c) measurement of $S$ in the experiment shown in Fig.~\ref{fig:setup},
        and (d) longer measurement at the optimal angle $\theta=8.6^\circ$.}
    \end{center}
\end{figure}

\section{Results}

The inequality is experimentally tested by evaluating $S$ in
Eq.~(\ref{eq:ineq3}) for a range of $\theta$; the
obtained values [points (c), (d) in Fig.~\ref{fig:violation_versus_separation}]
are
consistently lower than the trace (a) calculated via entropy
using Eq.~(\ref{eq:nid_entropy}),
and than a simulation with the same compressor (b).
This is because the LZMA Utility is not working exactly at the Shannon limit,
and also due to imperfect state generation and detection.

As a consequence of Eq.~(\ref{eq:braunstein}), we expect the maximal violation for $\theta=8.6^\circ$.
For this particular angle we collected results from a large number of photon
pairs. Although we set out in this work to avoid a statistical argument in the
interpretation of measurement results, we do resort to statistical techniques to
assess the confidence in an experimental finding of a violation of inequality
Eq.~(\ref{eq:ineq3}).
To estimate an uncertainty of the experimentally obtained values for $S$,
this large data set
was subdivided
into files with length greater than $10^5$ bits.
The results from all the subdivided files are then averaged to obtain the
final result of $S(\theta=8.6^\circ)$=\,0.0494$\,\pm\,$0.0076, with the latter
indicating a relatively small standard deviation over these different subsets.

\section{Conclusion}
There is a trend to look at physical systems and processes as
programs run on a computer made of the constituents of our universe. We could show
that this is not possible if one uses a computation paradigm of a local
deterministic Turing machine. Although this has been already extensively
researched in quantum information theory, we present a complementary
algorithmic approach for an explicit, experimentally testable example.
This algorithmic approach is complementary to the orthodox Bell inequality approach to quantum nonlocality~\cite{Bell} that is statistical in its nature.

Any process that can be simulated on a local universal Turing machine can be encoded as a program that is fed into it. For every such a program there exists its shortest description called Kolmogorov complexity, which in most of the cases can only be  approximated using compression software. Moreover, such a description must obey distance properties as shown in~\cite{Li:2004kz,Cilibrasi:2005jna}.
By testing Eq.~(\ref{eq:ineq3}), we showed that this is unattainable in the
specific case of polarization-entangled photon pairs. Therefore, there
exist physical processes that cannot be simulated on local universal Turing
machines.

There are two fundamentally different notions of complexity in computer
science. On one hand, computational complexity, mainly researched on in quantum information science, studies how much resources are needed to solve a computational problem. These studies focus on complexity classes such as P, NP~\cite{Papadimitriou:1994vk}, and its main concern is,  given an input program, how efficiently it can be computed. On the other hand, algorithmic complexity deals with a problem of what the most efficient encoding of an input program is. This complementary problem to computation complexity has not yet received enough attention
in quantum information science,
and it would require a further work on quantum version of Kolmogorov complexity~\cite{Mueller:2007wc}.

We would like to stress that our analysis of the experimental data is purely
and consistently algorithmic. We do not resort to statistical methods that are alien to the concept of computation. If this approach can be extended to all quantum experiments, it would allow us to bypass the commonly used statistical interpretation of quantum theory.

\section{Acknowledgments}
We acknowledge the support of this work by the National Research Foundation \&
Ministry of Education in Singapore, partly through the Academic Research Fund MOE2012-T3-1-009. P.K. and D.K. are also supported by the Foundational Questions Institute (FQXi).
A.C. also thanks Andrea Baronchelli for the hints on the use of compression software.

\section{Appendix: Symmetrization of detector efficiencies}\label{app:symmetrization}

The experimental setup (Fig.~\ref{fig:setup}) uses four APDs:
$D_{HA}$, $D_{VA}$ (Alice), and $D_{HB}$, $D_{VB}$ (Bob) to register photon pair events
in the two spatial modes. By denoting events at $D_{H}$ and $D_{V}$ as 1 and
0, the four possible output patterns are 00, 01, 10, and 11, where the least
and most significant bit corresponds to the Alice and Bob mode,
respectively. Due to differences in the the losses in the transmitted and
reflected port of the PBS, efficiencies in coupling light into the APDs, and
the quantum efficiencies of APDs, the detection efficiencies for the four
output combinations are different.
The resulting effective pair efficiencies are then given by the product of the
contributing detection efficiencies $\eta_{VB}$, $\eta_{HB}$, $\eta_{VA}$, and $\eta_{HA}$.

This asymmetry will skew the statistics of the measurement results.
We symmetrize the effective pair efficiencies for each  $(\theta_{A}, \theta_{B})$
measuring also the following settings:
$(\theta_{A}+45^\circ,\theta_{B})$,
$(\theta_{A}, \theta_{B}+45^\circ)$,
and $(\theta_{A}+45^\circ,\theta_{B}+45^\circ)$.
This procedure
swaps the output ports of
the PBS at which each
polarization state is detected.
The resulting outcomes are then interleaved,
providing an uniform detection
probability
for the four possible outcomes.
The effective pair detection efficiency for all four combinations is then
$(\eta_{VB}\,\eta_{VA}\,+\,\eta_{VB}\,\eta_{HA}\,+\,\eta_{HB}\,\eta_{VA}\,+\,\eta_{HB}\,\eta_{HA})/4$.

\section{Appendix: Choice of compressor}\label{app:choice_compression}
In order to evaluate the NCDs of the binary strings, we need to choose a
compression algorithm that performs close to the Shannon limit. This is
necessary to ensure that it does not introduce any unintended
artifacts that lead to an overestimation of the violation. Preferably we
want to work in the regime where the obtained NCDs
always underestimate the violation. For this purpose, we characterized
four compression algorithms implemented by freely available compression
programs:
\textit{lzma}~\cite{ref:lzma}, \textit{bzip2}~\cite{ref:bzip2}, \textit{gzip}~\cite{ref:gzip},
and \textit{lzw}~\cite{Welch:84}. To eliminate
the overhead associated with the compression of ASCII text files, we
save data in a binary format.

For this characterization and a simulation of the experiment, we need to
generate a ``random'' string of bits (1, 0) or pairs of bits (00, 01, 10, and
11) of various length with various probability distributions.
We generate these strings using the \emph{MATLAB}~\cite{ref:matlab} function
\mbox{\emph{randsample()}} that uses the pseudo random number generator
\emph{mt19937ar} with a long period of $2^{19937} - 1$. It is based on the
Mersenne Twister~\cite{matsumoto:98}, with ziggurat~\cite{marsaglia:00} as the
algorithm that generates the required probability distribution. The complexity
of this (deterministic) source of pseudorandom numbers should be high enough
to {\emph not} be captured as algorithmic.

The first part of this characterization involves establishing the minimum string length required for the compression algorithms to perform consistently.
We start by generating binary strings, $x$, with equal probability of 1's and 0's, i.e. random strings, of varying length.
For each $x$, we evaluate the compression overhead $Q$ as
\begin{equation}\label{eq:compression_overhead}
    Q=\frac{C(x)-H(x)} {l(x)}\,.
\end{equation}
For a good compressor, we expect $Q$ to be close to 0.
From Fig.~\ref{fig:compression_ratio_singles}, it can be seen
that for all the compressors, $Q$ starts to converge after about
$10^5$ bits, setting the minimum string length required for the compressors to
work consistently.
The \textit{lzw} compressor fails this test, converging to
a $Q$ of 0.37 for long string,
while \textit{bzip2}, \textit{gzip}, and \textit{lzma} give a $Q$
below $10^{-1}$.

\begin{figure}
  \begin{center}
    \includegraphics[width = 0.95\columnwidth]{figure5.eps}
        \caption{\label{fig:compression_ratio_singles}
        Comparison of the compression overhead $Q$ obtained using four different compression algorithms on pseudo-random strings of varying lengths. The expected value for an ideal compressor is 0.
         From this characterization we can exclude \textit{lzw} as a useful compressor for our application.
        }
    \end{center}
\end{figure}

In the second part of this characterization,
test the compressors with strings with a known amount of correlation.
We generate a random string $x$ of length
$10^7$ using the same technique already described.
We then generate a second string $y$ of equal length
and with probability $p$ of being correlated
to $x$.
For $p=0$ the two strings are equal, i.e. perfectly correlated.
For $p=0.5$ they are uncorrelated.

The two strings $x$ and $y$ are then combined to form the string $xy$:
to avoid artifacts due to the limited data block size of the compression
algorithms, the elements of $x$ and $y$ are interleaved.
We then compress $xy$ and evaluate the compression overhead $Q$ as a function of $p$.
The results for different compressors are shown in
Fig.~\ref{fig:compressibility_versus_correlation}. Although
there are ranges of $p$
where \textit{bzip} and \textit{gzip} perform better than \textit{lzma}, the latter shows a more uniform performance over the entire interval of $p$.
It is reasonable to assume that
the use of \textit{lzma}
should reduce the possibility of artifacts in the estimation of the NCD also for the data obtained from the experiment.
\begin{figure}
  \begin{center}
 \includegraphics[width = 0.95\columnwidth]{figure6.eps}
        \caption{
          Compression overhead $Q$ for the string $xy$ as a function of the probability of pairwise correlation $p$ between the bits of the generating strings $x$ and $y$
          for three different compressors: \textit{bzip}, \textit{gzip}, and \textit{lzma}.
        \label{fig:compressibility_versus_correlation}}
    \end{center}
\end{figure}

\bibliographystyle{apsrev4-1}


\section{Introduction and main results}
Motivated by the study  of the  bilinear Hilbert transform
\begin{equation}\label{BHT}
H(f,g)(x)=p.v.\frac{1}{\pi}\int_{\mathbb{R}}f(x-y)g(x+y)\frac{dy}{y},
\end{equation}
by Lacey and Thiele \cite{LT1}, \cite{LT2}, a substantial amount of work has been produced in the area of multilinear singular integral and multiplier operators. In this paper we study certain kinds of multi(sub)linear square functions.

We introduce a     bilinear operator which is closely related to the bilinear Hilbert transform.
Let $w$ be a cube in $\mathbb{R}^d$, $d\geq 1$. Then for functions $f,g$ in the Schwartz class $\mathcal{S}(\mathbb{R}^d)$, we define the bilinear operator $S_w$ associated with the symbol $\chi_w(\xi-\eta)$ as follows:
$$
S_w(f,g)(x)=\int_{\mathbb{R}^d}\int_{\mathbb{R}^d}\hat{f}(\xi)\hat{g}(\eta)\chi_w(\xi-\eta)e^{2\pi i x\cdot(\xi+\eta)}d\xi d\eta.
$$
We notice that when $d=1$ and $w$ is the characteristic function of a half plane in $\mathbb R^2$, then
$S_w$ is related to the bilinear Hilbert transform in \eqref{BHT}.

Let $\{w_l\}_{l\in\mathbb{Z}^d}$ be a sequence of disjoint cubes in $\mathbb{R}^d$. Let $S_{w_l}$ be the bilinear operator associated with the symbol $\chi_{w_l}(\xi-\eta)$ as defined above. Then for $f,g\in\mathcal{S}(\mathbb{R}^d)$, the {\it non-smooth bilinear Littlewood-Paley square function} associated with the sequence $\{w_l\}_{l\in\mathbb{Z}^d}$ is defined as
\begin{equation}\label{BSSF}
S(f,g)(x)=\Big(\sum\limits_{l\in\mathbb{Z}^d}|S_{w_l}(f,g)(x)|^2\Big)^{1/2}.
\end{equation}
A   {\it smooth } version of the
bilinear Littlewood-Paley square function   in \eqref{BSSF} can be obtained if the
characteristic function of the cube $w$ is replaced by a smooth bump adapted to $w$. Precisely, we let
\begin{equation}\label{BSSF2}
T(f,g)(x)=\Big(\sum\limits_{l\in\mathbb{Z}^d}|T_{\phi_l}(f,g)(x)|^2\Big)^{1/2},
\end{equation}
where $T_{\phi_l}$ is a bilinear operator associated with the smooth function $\phi_l$ whose Fourier transform is supported in $ w_l$, i.e.,
$$
T_{\phi_l}(f,g)(x)=\int_{\mathbb{R}^d}\int_{\mathbb{R}^d}\widehat{f}(\xi)\widehat{g}(\eta)\widehat{\phi_l}(\xi-\eta)e^{2\pi ix\cdot(\xi+\eta)}d\xi d\eta.
$$

A particular case of the bilinear Littlewood-Paley square function $T$ in \eqref{BSSF2} was   studied by Lacey \cite{L} who proved the following: \medskip\\
\noindent\textbf{Theorem A} (\cite{L}).
Let $\phi$ be a smooth function on $\mathbb{R}^d$ such that $\hat{\phi}$ is supported in the unit cube of $\mathbb{R}^d$. For $l\in\mathbb{Z}^d$, let $\widehat{\phi_l}$ be the function defined by $\widehat{\phi_l}(\xi)=\widehat{\phi}(\xi-l)$. Then for $2\leq p,q\leq\infty$ with $1/p+1/q=1/2$, there is a constant $C$ such that for $f,g\in\mathcal{S}(\mathbb{R}^d)$,
$$
\|T(f,g)\|_{L^2(\mathbb{R}^d)}\leq C\|f\|_{L^p(\mathbb{R}^d)}\|g\|_{L^q(\mathbb{R}^d)}.
$$

\medskip
In addition, Lacey raised two questions related to this:
(a) Does  {Theorem A} hold  for $r\neq 2$?
(b) Are there analogues for  non-smooth square functions?
The boundedness of non-smooth square functions on $\mathbb{R}\times \mathbb{R}$   strengthens the boundedness of bilinear Hilbert transform; in this direction and in the spirit of question (b), Bernicot \cite{B}, Ratnakumar and Shrivastava \cite{RS2} have provided some interesting results. As far as  question (a),  answers were provided by   Bernicot and  Shrivastava \cite{BS}, Mohanty and Shrivastava \cite{MS}.
The proofs of these results were based on  rather complicated  time-frequency analysis. However, Ratnakumar and Shrivastava \cite{RS1}   provided  a proof for the boundedness of a smooth bilinear square functions, which is based on more elementary techniques. Motivated by the work of \cite{RS1} and the increasing interest in multilinear operators, the aim of this paper is to study the $L^p$ boundedness properties of smooth $m$-(sub)linear Littlewood-Paley square functions.

\begin{definition}
We define an {\it $m$-(sub)linear Littlewood-Paley square function} as:
\begin{equation}\label{kkdfjd}
T(f_1,\dots,f_m)(x)=\Big(\sum\limits_{\ell\in\mathbb{Z}^d}|T_{\ell}(f_1,\dots,f_m)(x)|^2\Big)^{1/2},
\end{equation}
where
\begin{equation}\label{eq1}
T_{\ell}(f_1,\dots,f_m)(x)=\int_{\mathbb{R}^d}f_1(x-\theta_1 y)\cdots f_m(x-\theta_m y)K_{\ell}(y)dy{ .}
\end{equation}
Here $\theta_j$, $j=1,\dots,m$ are fixed, distinct and nonzero real numbers, and $K_{\ell}$ are integrable  functions that satisfy  certain additional size conditions.
\end{definition}

Expressing each $T_\ell$ in $m$-linear Fourier multiplier form, we write
\begin{align*}
&T_{\ell}(f_1,\dots,f_m)(x)=\int_{(\mathbb{R}^d)^m}\Big[\prod_{i=1}^m( \widehat{f_i}(\xi_i)
e^{2\pi i x\cdot\xi_i})\Big]\widehat{K_{\ell}}(\theta_1\xi_1+\cdots+\theta_m\xi_m)d\xi_1\cdots d\xi_m.
\end{align*}
Then $T_\ell$ is associated with the multiplier (or symbol)
$\widehat{K_{\ell}}(\theta_1\xi_1+\cdots+\theta_m\xi_m)$.

We use the notation $r'=r/(r-1)$ for the dual exponent of $r\in [1,\infty]$ with $1'=\infty$ and $\infty'=1$. Our main result is the following:
\begin{theorem}\label{thm1}
Let $1\le p<\infty$ and let $K$ be a measurable function on $\mathbb{R}^d$ which satisfies
\begin{equation}\label{kkjrff}
B_p=\sum\limits_{u\in \mathbb{Z}^d}\big\|\chi_{Q_u}K\big\|_{L^{p' }(\mathbb{R}^d)}<\infty  \qquad \textup{when $1\le p<2$}
\end{equation}
 and
$$
B_2=\sum\limits_{u\in \mathbb{Z}^d}\big\|\chi_{Q_u}K\big\|_{L^2(\mathbb{R}^d)} <\infty  \qquad \textup{when $2\le p<\infty$,}
$$
where $Q_u=u+[0, 1)^d$ for $u \in \mathbb{Z}^d$.
For $\ell\in \mathbb{Z}^d$   define $ {K_{\ell}}(x)= {K}(x )e^{2\pi i  x \cdot  \ell }$ and let $T_{\ell}$ be the $m$-linear Fourier multiplier     in \eqref{eq1}, where the $\theta_j$ are nonzero and distinct.  Then for $2\le p_j<\infty$, $j=1,\dots,m$,
there exists a positive constant $C=C(d, \theta_j  ,p_j)$ such that  the square function $T$ in \eqref{kkdfjd} satisfies
\begin{equation}\label{1.2}
\|T(f_1,\dots,f_m)\|_{L^p(\mathbb{R}^d)}\leq BC\prod\limits_{j=1}^m\|f_j\|_{L^{p_j}(\mathbb{R}^d)}
\end{equation}
for all functions $f_j$ in $L^{p_j}(\mathbb{R}^d)$, where $p$ and $p_j$ are related via
\begin{equation}\label{RES}
 \frac{1}{p}=\sum\limits_{j=1}^m\frac{1}{p_j}
\end{equation}
and $B$ in \eqref{1.2} is either $B_2$ if $2\le p<\infty  $ or $B_p$ when $1\le p<2$.
\end{theorem}

The case $m=2$ of this theorem was obtained  by Ratnakumar and Shrivastava \cite{RS1}.
The new ingredient of this note is the extension of this result to the case where $m\ge 3$, where a multilinear version of Young's inequality is needed. Such an inequality was obtained by
Guliyev and Nazirova \cite{GN} for a certain range of exponents which is not sufficient for our work and for this reason, we
provide an extension of  this paper.   To state the  multilinear Young  inequality of Guliyev and Nazirova, we fix $\theta_i$ nonzero distinct real numbers and we define
$$
\vec{f}\otimes g(x)=\int_{\mathbb{R}^d}f_1(x-\theta_1y)\cdots f_m(x-\theta_my)g(y)dy{ ,}
$$
where $\vec{f}=(f_1,\dots , f_m)$ and $f_j$, $g$ are measurable functions so that the preceding integral converges.
We denote by $L^{p,\infty}(\mathbb{R}^d)$ the weak $L^p$ space of all measurable functions $f$ such that
$
\|f\|_{L^{p,\infty}}=\sup\limits_{\alpha>0}\big\{\alpha|\{x\in\mathbb{R}^d:|f(x)|>\alpha\}|^{1/p}\big\}<\infty.
$

\begin{theorem}\label{lem1} (\cite{GN})
Assume that $1<r <\infty$,  $1< p_j \le \infty$, $i=1,\dots,m$,   $1/p=1/{p_1}+\cdots+1/{p_m}$,  $1/q+1=1/p+1/r$, and
$\frac{r'}{1+r'}\le  p<r'$ (equivalently $1\le q<\infty$). Then for $g\in L^{r,\infty}(\mathbb{R}^d)$ and
$f_j\in L^{p_j}(\mathbb{R}^d) $  we have $\vec{f}\otimes g\in L^q(\mathbb{R}^d)$ with    norm inequality
$$
\|\vec{f}\otimes g\|_{L^q(\mathbb{R}^d)}\leq C(d,\theta_j,p_j,r )\prod\limits_{i=1}^m\|f_i\|_{L^{p_i}(\mathbb{R}^d)}\|g\|_{L^{r,\infty}(\mathbb{R}^d)}.
$$
 \end{theorem}

Unfortunately,  the proof of Theorem  \ref{lem1} in \cite{GN} in the  range    $\frac{r'}{r'+1} \le p\le 1$
 (equivalently $1\le q\le r$) contains an inaccurate deduction (and the assertion is incorrect in this case).
In the proof of Theorem \ref{thm1} we need the case where $\frac{r'}{1+r'}=  p $ (equivalently $q=1$) and for this reason we
provide an analogue in which the weak norm of $g$ is replaced by a strong norm.
So we fix the inaccurate deduction in \cite{GN} via the following result:

 \begin{theorem}\label{lem11}
Assume that $1<r \le  \infty$,  $1/p=1/{p_1}+\cdots+1/{p_m}$, $1\le p_j \le \infty$, $j=1,\dots,m$,   $1/q+1=1/p+1/r$, and
$\frac{r'}{1+r'}\le  p\le r'$ (equivalently $1\le q \le \infty$). Then for some
constant $C(d,\theta_j,p_j,r )$  the following inequalities are valid:\\
(a) if $1<r<\infty$, then when $1<   p<r'$ (equivalently $r< q<\infty$) we have
$$
\|\vec{f}\otimes g\|_{L^q(\mathbb{R}^d)}\leq C(d,\theta_j,p_j,r )\prod\limits_{i=1}^m\|f_i\|_{L^{p_i}(\mathbb{R}^d)}\|g\|_{L^{r,\infty}(\mathbb{R}^d)}.
$$
(b) if  $1<r\le \infty$, then when $\frac{r'}{1+r'}\le  p\le 1$ (equivalently $1\le q\le r$) we have
$$
\|\vec{f}\otimes g\|_{L^q(\mathbb{R}^d)}\leq C(d, \theta_j,p_j ,r)\prod\limits_{i=1}^m\|f_i\|_{L^{p_i}(\mathbb{R}^d)}\|g\|_{L^{r }(\mathbb{R}^d)}.
$$
\end{theorem}

 \begin{remark}
 We notice that  in Theorem \ref{lem11} we are assuming that $p_j\ge 1$ for all $j$
whereas in \cite{GN} it is assumed that
$p_j>1$ for all $j$.   This small detail, i.e., the case $p_j=1$ for some $j$, does not appear in case (a), since  we
always have $p\le p_j$ for all $j$ and in case (a) we have $p >1$. Thus the case where some $p_j$ equals $1$
is relevant only  in case (b), which is treated in the next section.
  \end{remark}

 \section{The proof of Theorem \ref{lem11}}

The case where $\frac{r'}{r'+1} \le p\le 1$, i.e.,
case (a) of Theorem \ref{lem11} is contained in \cite{GN}.
So we focus attention in  case (b)   of Theorem \ref{lem11}, i.e., when
$\frac{r'}{1+r'}\le  p\le r'$ (equivalently $1\le q\le r$). Below we assume that $1<r\le \infty$.

First we consider the case $p=1$ when $q=r$. Since $1/p_1+\cdots+1/p_m=1$,
  making use of H\"older's inequality and Fubini's theorem, we write
\begin{align*}
&\|\vec{f}\otimes g\|_{L^r(\mathbb{R}^d)}\\
&\leq \Big(\int_{\mathbb{R}^d}\Big(\int_{\mathbb{R}^d}|f_1(x-\theta_1 y)\cdots f_m(x-\theta_m y)g(y)|dy\Big)^rdx\Big)^{\frac 1 r}\\
&\leq \Big(\int_{\mathbb{R}^d}\Big\{\Big(\int_{\mathbb{R}^d}|g(y)|^r |f_1(x-\theta_1 y)\cdots f_m(x-\theta_m y)|dy\Big)^{\frac 1 r}\\
&\ \ \ \qquad \qquad\qquad \qquad \qquad \Big(\int_{\mathbb{R}^d}|f_1(x-\theta_1 y)\cdots f_m(x-\theta_m y)|dy\Big)^{\frac 1 {r'}}\Big\}^rdx\Big)^{\frac 1 r}\\
&\leq \bigg(\prod_{i=1}^m |\theta_i|^{-\frac{d}{p_i r'}}\bigg)\Big(\int_{\mathbb{R}^d} \int_{\mathbb{R}^d}|g(y)|^r |f_1(x-\theta_1 y)\cdots f_m(x-\theta_m y)|dy\,  dx\Big)^{\frac 1 r}\prod\limits_{i=1}^m\|f_i\|_{L^{p_i}}^{\frac 1{r'}}\\
&= \bigg(\prod_{i=1}^m |\theta_i|^{-\frac{d}{p_i r'}}\bigg) \Big(\int_{\mathbb{R}^d}|g(y)|^r\int_{\mathbb{R}^d}|f_1(x-\theta_1 y)\cdots f_m(x-\theta_m y)|dx\, dy\Big)^{\frac 1 r}\prod\limits_{i=1}^m\|f_i\|_{L^{p_i}}^{\frac 1{r'}}\\
&\leq \bigg(\prod_{i=1}^m |\theta_i|^{-\frac{d}{p_i r'}}\bigg)\|g\|_{L^r}\prod\limits_{i=1}^m\|f_i\|_{L^{p_i}}^{\frac 1r}\prod\limits_{i=1}^m\|f_i\|_{L^{p_i}}^{\frac 1{r'}}\\
&=\bigg(\prod_{i=1}^m |\theta_i|^{-\frac{d}{p_i r'}}\bigg)\|g\|_{L^r}\prod\limits_{i=1}^m\|f_i\|_{L^{p_i}},
\end{align*}
where in the second inequality we apply H\"older's inequality with respect to the measure $\prod\limits_{i=1}^m|f_i(x-\theta_i y)|dy$ to the function $y\mapsto|g(y)|$ and 1 with exponents $r$ and $r'$,  respectively. Notice that this proof also
works in the case where $r=\infty$ with a small notational modification.

Next we consider the case $p=\frac{r'}{r'+1}$, which is equivalent to $q=1$. In this case $r$ could be equal to $\infty$ and thus  $r'$ could take the value   $1$ when $p=1/2$.
Let us give the proof first in the case where $m=3$, where the notation is simpler. In this case when $r'>1$
notice that $1/p_1+1/p_2+1/p_3=1+1/r'<2$, thus in the  space  with coordinates $(x_1,x_2,x_3)$, the intersection of the plane $x_1+x_2+x_3=1+1/r'$ and of the boundary of the unit cube $[0,1]^3$
is a closed hexagon with vertices
at six points $V_1^{(3)}=(1/r',0,1)$, $V_2^{(3)}=(1/r',1,0)$, $V_3^{(3)}=(1,0,1/r')$, $V_4^{(3)}=(0,1,1/r')$, $V_5^{(3)}=(1,1/r',0)$, $V_6^{(3)}=(0,1/r',1)$.
In the case where $r=\infty$, i.e., $r'=1$,   the intersection of boundary of the
unit cube $[0,1]^3$ with the plane $x_1+x_2+x_3= 2$ is a closed triangle with vertices
the points $W_1^{(3)}=(0,1,1)$, $W_2^{(3)}=(1,0,1)$ and $W_3^{(3)}=(1,1,0)$.

We assert that $\vec{f}\otimes g$ maps $L^{p_1} \times L^{p_2}\times L^{p_3}\times L^r$ into $L^1$ when $(1/p_1,1/p_2,1/p_3)$ lies in the closed convex hull of the  six points $V_j^{(3)}$ when $r<\infty$ (or the three points $W_1^{(3)}$, $W_2^{(3)}$, $W_3^{(3)}$ when $r=\infty$).   This assertion will be a consequence of the fact that boundedness holds at the   vertices via multilinear complex interpolation (in particular by applying the theorem in  Zygmund  \cite[Chapter XII, (3.3)]{Zyg} or in Berg and L\"ofstrom   \cite[Theorem 4.4.2]{BL} or   the main result in \cite{GrMa}).

 We only prove boundedness for point $V_1^{(3)}$, since the method is similar for the remaining vertices. We have
\begin{align*}
&\|\vec{f}\otimes g\|_{L^1(\mathbb{R}^d)}\\
&\leq \int_{\mathbb{R}^d}\int_{\mathbb{R}^d}|f_1(x-\theta_1 y)||f_2(x-\theta_2 y)||f_3(x-\theta_3 y)||g(y)|dydx\\
& =\int_{\mathbb{R}^d}|g(y)|\int_{\mathbb{R}^d}|f_1(x-\theta_1 y)||f_2(x-\theta_2 y)||f_3(x-\theta_3 y)|dxdy\\
&\leq \Big(\int_{\mathbb{R}^d}|g(y)|^rdy\Big)^{ \frac 1r}\Big(\int_{\mathbb{R}^d}\Big(\int_{\mathbb{R}^d}|f_1(x-\theta_1 y)||f_2(x-\theta_2 y)||f_3(x-\theta_3 y)|dx\Big)^{r'}dy\Big)^{\frac 1{r'}}\\
&\leq \|g\|_{L^r}\|f_2\|_{L^\infty}\Big(\int_{\mathbb{R}^d}\Big(\int_{\mathbb{R}^d}|f_1(t)|
|f_3(t-(\theta_3-\theta_2)y)|dt\Big)^{r'}dy\Big)^{\frac 1{r'}}\\
&\leq |\theta_3-\theta_2|^{-\frac d{r'}} \|g\|_{L^r}\|f_1\|_{L^{r'}}\|f_2\|_{L^\infty}\|f_3\|_{L^1},
\end{align*}
where in the last inequality we use Young's inequality,  in view of the fact that    $p_1=r'$ and $p_3=1$ for the  point $V_1^{(3)}$,   which implies $1/p_1+1/p_3=1+1/r'$. Note that $r$ could be equal to infinity in this case.
Boundedness for the  remaining
vertices follows by symmetry.

The case $m\ge 4$ is proved via the same method. Let $E^{(m)}$ be the intersection of the boundary     of the cube $[0,1]^m$  in $\mathbb R^m$ with the hyperplane $x_1+\cdots + x_m=1+1/r'$. Let us first consider the case where $r'>1$. The set $E^{(m)}$ is a closed convex polygon with $m(m-1) $ vertices which we call $V_i^{(m)}$. These are the   points
 in $\mathbb R^m$ whose coordinates are all zero except for two of them which are    $1$  and $1/r'$.
We claim that $E^{(m)}$ is the closed convex hull of the vertices $V_i^{(m)}$. In fact, we   show this by induction, as follows:
The assertion is true when  $m=3$. If this assertion is valid when $m=n-1$ for some $n\ge 4$, then
the faces of $E^{(n)}$ are the sets
$$
\{(1/q_1, \dots , 1/q_n):\,\, \sum_{i=1}^n q_i^{-1} = 1+ 1/r' \, , \quad \textup{$1\le q_i\le \infty$,  and   $q_i=\infty$ for only one $i$} \}.
$$
 Consider for instance the face
$$
F_1^{(n)}=\{(1/q_1, \dots , 1/q_n):\,\, \sum_{i=1}^n q_i^{-1} = 1+ 1/r' \, , \quad \textup{$1\le q_i\le \infty$ and   $q_1=\infty$} \}.
$$
 Then $F_1^{(n)}$ can be identified with the set
 $$
 \{(1/q_2, \dots , 1/q_n):\,\, \sum_{i=2}^n q_i^{-1} = 1+ 1/r' \, ,  \quad \textup{$1\le q_i\le \infty$} \}
 $$
which is equal to the intersection of the cube $[0,1]^{m-1}$ with the hyperplane $\sum_{i=2}^n x_i  = 1+ 1/r'$ in
$\mathbb R^{m-1}$
and  corresponds to the case $m=n-1$. Applying the    induction
hypothesis, the vertices of $F_1^{(n)}$ are the points $(1/q_1, \dots , 1/q_n)$  where $1/q_1=0$ and the
$1/q_i$ with $i\ge 2$ are   zero except for two of them which are either $1$ or $1/r'$. Thus, since $1/q_1=0$,  the vertices of
$F_1^{(n)}$ are the points $(1/q_1, \dots , 1/q_n)$  where $1/q_1=0$ and
$1/q_i$ with $i\ge 2$ are   zero except for two of them which are either
$1$ or $1/r'$. The same holds for the remaining faces and thus
the assertion that the points $V_i^{(m)}$ are the vertices of  $E^{(m)}$ holds when $m=n$. The same conclusion holds when
$r'=1$, but in this case the  intersection of the hyperplane $x_1+\cdots + x_m=2$ with the cube $[0,1]^m$  is a closed convex set with  $ m(m-1)/2$ vertices $W_i^{(m)}$ instead of $m(m-1)$ vertices.  Observe that the set of
$W_i^{(m)}$ is the set of all vectors with in $\mathbb R^m$ whose coordinates are zero except for two of them which are equal to $1$.

\begin{comment}
since the set $E$ is a closed convex polygon, there exists a line $l$ such that two vertices of set $E$ lie on $l$, i.e., one edge of $E$ is on $l$, whereas other vertices of set $E$ lie on the same side of $l$. Next, we pick one of two directions of $l$ as its positive direction. Assume all the vertices of set $E$ that not on $l$ are on the left side of $l$ according to this positive direction. Find the latter vertex, which we assume is $V_1$, of $E$ that lies on $l$ according to the positive direction. Then we rotate $l$ around $V_1$ anti-clockwise until $l$ meets other vertex of $E$. Assume this vertex is $V_2$, and we again rotate $l$ around $V_2$ anti-clockwise. Continuing these steps repeatedly until $l$ meets $V_1$ when rotating $l$ around $V_{m(m-1)}$ anti-clockwise. Notice that the step could be ended since $E$ has finite vertices. Thus, convex polygon formed by points $V_1,\dots, V_{m(m-1)}$ is the closed convex hull of vertices $V_1,\dots, V_{m(m-1)}$, while this convex polygon is the set $E$. Therefore, we prove our claim.
\end{comment}

 Notice that at each vertex $V_i^{(m)}$ (or $W_i^{(m)}$)
 boundedness holds and is proved
exactly in the same way as when $m=3$ with the only exception being that  the $L^\infty$ norm  of one function  is replaced by
the product of $L^\infty$ norms of the remaining $m-2$ functions. Then multilinear complex interpolation yields
boundedness in the closed convex hull of these vertices, i.e., yields
the required conclusion when $p=\frac{r'}{r'+1}$ (equivalently $q=1$).

Finally, the case where $\frac{r'}{1+r'}<  p< 1$ (equivalently $1< q< r$)  is obtained by
interpolation between the endpoint cases where
$p=\frac{r'}{r'+1}$ (equivalently $q=1$) and $p=1$ (equivalently $q=r$).   Notice that in the case where $p=1$ ($q=r$) the
intersection of the boundary of the unit cube $[0,1]^m$ with the hyperplane $x_1+\cdots + x_m=1$ is the
closed convex hull of the $m$ points $U_j^{(m)}=(0,\dots , 1, \dots , 0)$ with $1$ in the $j$th coordinate and zero in every other coordinate.
Let $Z_k^{(m)}$ be the vertices of the closed convex set formed by the intersection of the boundary of the unit cube $[0,1]^m$ with the hyperplane $x_1+\cdots + x_m=1/q+1/r'$ in $\mathbb R^m$ when $1<q<r$. Then there is a $\theta \in (0,1)$
[in fact $\theta = r'(\frac 1q -\frac 1r)$] such that
the set of all $Z_k^{(m)}$ is contained in the set
$$
\Big\{(1-\theta) U_j^{(m)} + \theta V_i^{(m)}{:} \qquad  j\in \{1,\dots ,m\}, \ \ i \in \{1,\dots , m(m-1)\}\Big\},
$$
where $V_i^{(m)}$ is replaced by $W_i^{(m)}$ for $i\in \{1,\dots , m(m-1)/2\}$,  if $r=\infty$.  Thus interpolation is possible in this case.

\begin{comment}
Then we assume the boundedness estimates holds for $m-1$, that is,
$$\|\vec{f}\otimes g\|_{L^1}\leq C\prod\limits_{i=1}^{m-1}\|f_i\|_{L^{p_i}}\|g\|_{L^r}.$$
Now we consider the $m$-linear case.
By assumption we have that
\begin{align*}
&\|\vec{f}\otimes g\|_{L^1}\leq \|f_m\|_{L^{\infty}}\int_{\mathbb{R}^d}\int_{\mathbb{R}^d}\prod\limits_{i=1}^{m-1}|f_i(x-\theta_i y)||g(y)|dydx\\
&\ \ \ \ \ \ \ \ \ \ \ \ \ \leq \|f_m\|_{L^\infty}\prod\limits_{i=1}^{m-1}\|f_i\|_{L^{p_i}}\|g\|_{L^r}.
\end{align*}
Thus, we obtain the boundedness estimates on a vertex of convex hull in $m$ dimension. The same is true for all other vertices. Then by interpolation, we get $\vec{f}\otimes g$ is bounded from $L^{p_i}\times\cdots\times L^{p_m} \times L^r$ to $L^1$ if $(1/p_1,\cdots,1/p_m)$ lies in this convex hull. Therefore, we have obtained required estimates for the case $q=1$. At last applying interpolation between the case $q=r$ and $q=1$ yields $\vec{f}\otimes g$ maps $L^{p_1}\times\cdots\times L^{p_m} \times L^r$ into $L^q$ when $1\leq q\leq r$.\\
Therefore, we complete the proof of Theorem \ref{lem11}.

\section{Requisite lemmas}
The following lemma will be used as a key ingredient in our argument.

\begin{lemma}\label{lem1}
Assume that $1<r<\infty$, $g\in L^{r,\infty}(\mathbb{R}^d)$, $1/p=1/{p_1}+\cdots+1/{p_m}$, $p_i>1$, $i=1,\dots,m$ and $1/q+1=1/p+1/r$. If $1/m< p<r/(r-1)$ (equivalently $r/(mr-r+1)<q<\infty$). Then, for $\vec{f}\in L^{p_1}(\mathbb{R}^d)\times\cdots\times L^{p_m}(\mathbb{R}^d)$, we have $\vec{f}\otimes g\in L^q(\mathbb{R}^d)$ and the norm inequality as follows,
$$\|\vec{f}\otimes g\|_{L^q(\mathbb{R}^d)}\leq C(\theta,p,q,r)\prod\limits_{i=1}^m\|f_i\|_{L^{p_i}(\mathbb{R}^d)}\|g\|_{L^{r,\infty}(\mathbb{R}^d)}.$$
where $\vec{f}\otimes g(x)=\int_{\mathbb{R}^d}f_1(x-\theta_1y)\cdots f_m(x-\theta_my)g(y)dy$ and $\theta_i$ ($i=1,\cdots,m$) are defined as   above.
\end{lemma}

We notice that Lemma \ref{lem1} extends the range of $p$ and $q$ in \cite{GN}, where $p$ and $q$ satisfy $r/(2r-1)\leq p<r/(r-1)$ and $1\leq q<\infty$. We will take the similar method as in \cite{GN} to prove it. In the proof of \cite{GN}, $p$ is divided into two cases $1<p<r/(r-1)$ and $r/(2r-1)\leq p\leq 1$. Since the upper bound of $p$ in Lemma \ref{lem1} remains the same as in \cite{GN}, we only prove that Lemma \ref{lem1} holds in the case $1/m<p\leq 1$, which implies $q$ could be less than $1$.

We will need the following definitions and lemmas in the proof of Lemma \ref{lem1}.
\begin{definition}
The distribution function of a measurable function $f$ on $\mathbb{R}^d$ is defined by
$$
\lambda_f(s)=|\{x\in\mathbb{R}^d:|f(x)|>s\}|, s\geq 0.
$$
If $f$ is finite almost everywhere on $\mathbb{R}^d$ and $\lambda_f(s)<\infty$ for all $s>0$, then the rearrangement $f^*$ of $f$ and its vector $\vec{f}^*$ are defined as
$$f^*(t)=\inf\{s>0:\lambda_f(s)\leq t\}, t\geq 0,$$ and
$$\vec{f}^*(t)=f_1^*(t)\cdots f_m^*(t).$$ Moreover, we denote the maximal function of $\vec{f}^*$ by
$$\vec{f}^{**}(t)=\frac{1}{t}\int_0^t f_1^*(x)\cdots f_m^*(x)dx,\quad t>0.$$
\end{definition}
We have the following alternative description of the $L^p$ norm in terms of the distribution function and the non-increasing rearrangement, see \cite{BS1}: if $0<p<\infty$, then
$$\int_{\mathbb{R}^d}|f(x)|^pdx=\int_0^\infty f^*(t)^pdt.$$\\
We also need the following lemmas.
\begin{lemma}\cite{GN}\label{lem2}
Let $0<p\leq 1$, $p<q<\infty$, $1<r<\infty$. Then the inequality
$$
\Big(\int_0^\infty\Big(\int_t^\infty\tau^{-\frac{1}{r}}\phi(\tau)d\tau\Big)^qdt\Big)^{1/q}\leq C(r,q)\Big(\int_0^\infty\phi(t)^pdt\Big)^{1/p}$$
holds for all non-negative and non-increasing functions $\phi$ if and only if $1/q+1=1/p+1/r$.
\end{lemma}

\begin{lemma}\cite{GN}\label{lem3}
Let $0<p\leq 1$, then for all non-negative and non-increasing functions $\phi$ the following inequality
$$\int_0^\infty \Big(\frac{1}{t}\int_0^t\phi(\tau)d\tau\Big)^pdt\leq p/(p-1)\int_0^\infty \phi(t)^pdt$$
holds.
\end{lemma}

\begin{lemma}\cite{GN}\label{lem4}
Let $f_i$, $i=1,\dots,m$ be functions that are finite almost everywhere in $\mathbb{R}^d$ and $\lambda_{f_i}(s)<\infty$, $i=1,\dots,m$ for all $s>0$, $g\in L^{r,\infty}(\mathbb{R}^d)$, $1<r<\infty$. Then
$$(\vec{f}\otimes g)^*(t)\leq (\vec{f}\otimes g)^{**}(t)\leq r/(r-1)C(\theta) \|g\|_{L^{r,\infty}}\int_t^\infty \tau^{-\frac{1}{r}}\vec{f}^{**}(\tau)d\tau,$$
where $C(\theta)=|\theta_1\cdots \theta_m|^{-n}$.
\end{lemma}
 Now we are in the position to give the proof of Lemma \ref{lem1} for the case $1/m<p\leq 1$.
By making use of Lemma \ref{lem4}, Lemma \ref{lem2}, Lemma \ref{lem3} and H\"{o}lder's inequality, we obtain
\begin{align*}
&\|\vec{f}\otimes g\|_{L^q(\mathbb{R}^d)}=\|(\vec{f}\otimes g)^*\|_{L^q(0,\infty)}\\
&\leq \|(\vec{f}\otimes g)^{**}\|_{L^q(0,\infty)}\\
&\leq r/(r-1)C(\theta)\Big(\int_0^\infty\Big|\int_t^\infty \tau^{-\frac{1}{r}}\vec{f}^{**}(\tau)d\tau\Big|^qdt\Big)^{1/q}\|g\|_{L^{r,\infty}}\\
&\leq r/(r-1)C(\theta)C(r,q)\Big(\int_0^\infty \vec{f}^{**}(t)^pdt\Big)^{1/p}\|g\|_{L^{r,\infty}}\\
&\leq (p/(p-1))^{1/p}r/(r-1)C(\theta)C(r,q)\Big(\int_0^\infty \vec{f}^{*}(t)^pdt\Big)^{1/p}\|g\|_{L^{r,\infty}}\\
&\leq (p/(p-1))^{1/p}r/(r-1)C(\theta)C(r,q)\prod\limits_{i=1}^m\|f_i^*\|_{L^{p_i}(0,\infty)}\|g\|_{L^{r,\infty}}\\
&=(p/(p-1))^{1/p}r/(r-1)C(\theta)C(r,q)\prod\limits_{i=1}^m\|f_i\|_{L^{p_i}(\mathbb{R}^d)}\|g\|_{L^{r,\infty}}.
\end{align*}
Thus, we prove Lemma \ref{lem1} is true for the case $1/m<p\leq 1$. For the proof of the case $1<p<r/(r-1)$, it is the same as in \cite{GN}. Therefore, we extend the range of $p$ and $q$ in \cite{GN}.

\end{comment}

\section{Proof of Theorem \ref{thm1} }

\begin{proof} Using the inverse Fourier transform we write
\begin{align*}
&T_{\ell}(f_1,\dots,f_m)(x)=\int_{\mathbb{R}^d}f_1(x-\theta_1y)\cdots f_m(x-\theta_my)K_{\ell}(y)dy\\
&\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ =\int_{\mathbb{R}^d}f_1(x-\theta_1y)\cdots f_m(x-\theta_my)e^{2\pi i\ell\cdot y}K(y)dy.
\end{align*}

We claim that the following inequality holds for $T$:
\begin{equation}\aligned \label{2.1}
T(f_1,\dots,f_m)(x)&=\Big(\sum\limits_{\ell\in \mathbb{Z}^d}|T_{\ell}(f_1,\dots,f_m)(x)|^2\Big)^{\frac{1}{2}}\\
&\leq \sum\limits_{u\in \mathbb{Z}^d}\bigg(\int_{Q_u}|f_1(x-\theta_1y)\cdots f_m(x-\theta_my)K(y)|^2dy\bigg)^{\frac{1}{2}}.
\endaligned\end{equation}
To verify the validity of this claim, we   take a sequence $a=\{a_{\ell}\}_{\ell\in\mathbb{Z}^d}\in \ell_2(\mathbb{Z}^d)$ with the property $\|a\|_{\ell_2(\mathbb{Z}^d)}=1$. Then by duality,  we need to prove   that
$$
\Big|\sum\limits_{\ell\in\mathbb{Z}^d}a_{\ell}T_{\ell}(f_1,\dots,f_m)(x)\Big|\leq \sum\limits_{u\in \mathbb{Z}^d}\bigg(\int_{Q_u}|f_1(x-\theta_1y)\cdots f_m(x-\theta_my)K(y)|^2dy\bigg)^{\frac{1}{2}}.
$$
Note that
\begin{align*}
&\sum\limits_{\ell\in\mathbb{Z}^d}a_{\ell}T_{\ell}(f_1,\dots,f_m)(x)=\int_{\mathbb{R}^d}f_1(x-\theta_1y)\cdots f_m(x-\theta_my)K(y)\sum\limits_{\ell\in\mathbb{Z}^d}a_{\ell}
e^{2\pi i\ell\cdot y}dy\\
&\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ =\int_{\mathbb{R}^d}f_1(x-\theta_1y)\cdots f_m(x-\theta_my)K(y)\hat{a}(y)dy{ ,}
\end{align*}
where $\hat{a}$ is the Fourier transform of sequence $a$. It is easy to see that $\hat{a}$ is a periodic function with $\|\hat{a}\|_{L^2([0,1]^d)}=1$. Hence,
\begin{align*}
&\Big|\sum\limits_{\ell\in\mathbb{Z}^d}a_{\ell}T_{\ell}(f_1,\dots,f_m)(x)\Big|\\
&\leq \int_{\mathbb{R}^d}|f_1(x-\theta_1y)\cdots f_m(x-\theta_my)K(y)\hat{a}(y)|dy\\
& =\sum\limits_{u\in\mathbb{Z}^d}\int_{Q_u}|f_1(x-\theta_1y)\cdots f_m(x-\theta_my)K(y)\hat{a}(y)|dy\\
&\leq \sum\limits_{u\in\mathbb{Z}^d}\Big(\int_{Q_u}|f_1(x-\theta_1y)\cdots f_m(x-\theta_my)K(y)|^2dy\Big)^{1/2}\bigg(\int_{Q_u}|\hat{a}(y)|^2dy\bigg)^{1/2}\\
& =\sum\limits_{u\in\mathbb{Z}^d}\bigg(\int_{Q_u}|f_1(x-\theta_1y)\cdots f_m(x-\theta_my)K(y)|^2dy\bigg)^{1/2}.
\end{align*}
Therefore, we proved the claimed inequality (\ref{2.1}).

To show inequality (\ref{1.2}), we consider two cases:

\noindent\textbf{Case 1: $  2\le p<\infty$}.
Making use of Minkowski's inequality  we obtain
\begin{align*}
&\|T(f_1,\dots,f_m)\|_{L^p(\mathbb{R}^d)}\\
&\leq\sum\limits_{u\in\mathbb{Z}^d}\Big(\int_{\mathbb{R}^d}\Big(\int_{Q_u}|f_1(x-\theta_1y)\cdots f_m(x-\theta_my)K(y)|^2dy\Big)^{p/2}dx\Big)^{1/p}\\
&\leq\sum\limits_{u\in\mathbb{Z}^d}\Big(\int_{Q_u}\Big(\int_{\mathbb{R}^d}|f_1(x-\theta_1y)\cdots f_m(x-\theta_my)K(y)|^pdx\Big)^{2/p}dy\Big)^{1/2}\\
&=\sum\limits_{u\in\mathbb{Z}^d}\Big(\int_{Q_u}|K(y)|^2\Big(\int_{\mathbb{R}^d}|f_1(x-\theta_1y)\cdots f_m(x-\theta_my)|^pdx\Big)^{2/p}dy\Big)^{1/2}\\
&\leq\sum\limits_{u\in\mathbb{Z}^d}\Big\{\int_{Q_u}|K(y)|^2\Big[\Big(\int_{\mathbb{R}^d}|f_1|^{p\cdot\frac{p_1}{p}}dx\Big)^{\frac{p}{p_1}}
\cdots\Big(\int_{\mathbb{R}^d}|f_m|^{p\cdot\frac{p_m}{p}}dx\Big)^{\frac{p}{p_m}}\Big]^{\frac{2}{p}}dy\Big\}^{\frac{1}{2}}\\
&=B_2\prod\limits_{i=1}^m\|f_i\|_{L^{p_i}(\mathbb{R}^d)}.
\end{align*}

\noindent\textbf{Case 2: $1\le  p<2$}. By Minkowski's inequality, we obtain
\begin{align*}
&\|T(f_1,\dots,f_m)\|_{L^p(\mathbb{R}^d)}\\
&\leq \Big(\int_{\mathbb{R}^d}\Big|\sum\limits_{u\in\mathbb{Z}^d}\Big(\int_{Q_u}|f_1(x-\theta_1y)\cdots f_m(x-\theta_my)K(y)|^2dy\Big)^{1/2}\Big|^pdx\Big)^{1/p}\\
&\leq \sum\limits_{u\in\mathbb{Z}^d}\Big(\int_{\mathbb{R}^d}\Big(\int_{Q_u}|f_1(x-\theta_1y)\cdots f_m(x-\theta_my)K(y)|^2dy\Big)^{p/2}dx\Big)^{1/p}\\
&\leq\sum\limits_{u\in\mathbb{Z}^d}\Big\{\sum\limits_{n\in\mathbb{Z}^d}\Big(\int_{P_n}\int_{Q_u}|f_1(x-\theta_1y)\cdots f_m(x-\theta_my)K(y)|^2dydx\Big)^{p/2}\Big\}^{1/p}{ ,}
\end{align*}
where $P_n= n+[0, 1)^d$ and we used H\"older's inequality for exponents $2/p$ and $(2/p)'$ over the space $(P_n, dx)$ in the last inequality.

Let
$$
A_{n,u}=\int_{P_n}\int_{Q_u}|f_1(x-\theta_1y)\cdots f_m(x-\theta_my)K(y)|^2dydx
$$
and set   $f_{i}^{n,u} =f_i \chi_{P_n-\theta_i Q_u} $, $i=1,2,\dots,m$, $k_u =K \chi_{Q_u} $ and $\vec{f^{n,u}}=  (f_{i}^{n,u})_{1\le i\le m}$.
Using Theorem \ref{lem11} (part (b)) with $q=1$, $r=p'/2$, and $p_i/2$ in place of $p_i$, we obtain
\begin{align*}
A_{n,u}&\leq\int_{\mathbb{R}^d}\int_{\mathbb{R}^d}|f_{1}^{n,u}(x-\theta_1y)\cdots f_{m}^{n,u}(x-\theta_my)k_u(y)|^2dydx\\
&= \int_{\mathbb{R}^d}|\vec{f^{n,u}}|^2\otimes|k_u|^2(x)dx \\
&\leq C\prod\limits_{i=1}^m\| |f_{i}^{n,u}|^2\|_{L^{p_i/2}(\mathbb{R}^d)}\| |k_u|^2\|_{L^{p'/2  }(\mathbb{R}^d)}\\
&=C\prod\limits_{i=1}^m\|f_{i}^{n,u}\|_{L^{p_i}(\mathbb{R}^d)}^2\|  k_u \|_{L^{p'   }(\mathbb{R}^d)}^2.
\end{align*}
Therefore, it follows from H\"older's inequality that
\begin{align*}
&\|T(f_1,\dots,f_m)\|_{L^p(\mathbb{R}^d)}\\
&\leq C\sum\limits_{u\in\mathbb{Z}^d}\Big(\sum\limits_{n\in\mathbb{Z}^d}\Big(\prod\limits_{i=1}^m
\|f_{i}^{n,u}\|_{L^{p_i}(\mathbb{R}^d)}^2\|k_u\|^2_{L^{p'  }(\mathbb{R}^d)}\Big)^{p/2}\Big)^{1/p}\\
&= C\sum\limits_{u\in\mathbb{Z}^d}\|k_u\|_{L^{p'   }(\mathbb{R}^d)}\Big(\sum\limits_{n\in\mathbb{Z}^d}
\Big(\int_{\mathbb{R}^d}|f_{1}^{n,u}|^{p_1}dx\Big)^{\frac p{p_1}}\cdots\Big(\int_{\mathbb{R}^d}|f_{m}^{n,u}|^{p_m}dx
\Big)^{\frac p{p_m}}\Big)^{\frac 1p}\\
&\leq C\sum\limits_{u\in\mathbb{Z}^d}\|K \chi_{Q_u}\|_{L^{p'  }(\mathbb{R}^d)} \Big\{\prod_{j=1}^m \Big(\sum\limits_{n\in\mathbb{Z}^d}
\Big(\int_{P_n-\theta_j Q_u}|f_j|^{p_j}dx\Big)^{\frac{p}{p_j}\frac{p_j}{p}}\Big)^{\frac{p}{p_j}} \Big\}^{\frac{1}{p}}\\
&\leq C\, B_p   \prod\limits_{i=1}^m\|f_i\|_{L^{p_i}(\mathbb{R}^d)} .
\end{align*}

Thus, we complete the proof of Theorem \ref{thm1}.

\par\bigskip
\end{proof}