\section{Introduction}
\label{Intro}

Quantum correlations play a central role as a resource in quantum information processing and quantum communication tasks.  Traditionally, entanglement was thought of as the unique form of quantum correlation and the reason why quantum computers can outperform classical computers.
Yet there are tasks that are believed to be exponentially hard classically, which can be done efficiently using quantum computational models with little or no entanglement~\cite{KL98,DFCaves05}.  Quantum discord was introduced as a more general measure of quantum correlation for bipartite systems, with no classical analogue~\cite{Zurek01}, and it was suggested as a resource for certain quantum computation models~\cite{DSCaves08}, quantum state merging~\cite{madhok11,cavalcanti}, and for encoding information onto a quantum state~\cite{gu}.  Discord has been generalized to continuous-variable systems to study quantum correlations in Gaussian states~\cite{Datta10,Paris10} and certain nonGaussian states~\cite{TMAdessoK12}.

Recently, schemes have been proposed to test for nonvanishing quantum discord of discrete-variable quantum states~\cite{terno2010,rahimi2010,BC10,DVBrukner10,modi2011,ZYCO11,sarandy2011,long2011,GAdesso12,serra2012}, and some of these have been implemented in nuclear-magnetic-resonance systems~\cite{Serra11,Laflamme11} and in an optical system~\cite{walborn2012}.
Of particular practical interest, however, is a general method for detecting nonvanishing discord in the joint state of both discrete and continuous-variable systems.

In this paper we introduce a measurement-based method for verifying quantum discord of any bipartite quantum state, without requiring any prior knowledge of the joint state.
We consider the post-measurement states of one subsystem, $B$, conditioned to all the outcomes of an informationally complete POVM (IC-POVM) performed on the other subsystem, $A$.  We show that if the post-measurement states of $B$ commute with one another, then the quantum discord from $B$ to $A$ is zero. Conversely, if they do not commute, the quantum discord from $B$ to $A$ is necessarily nonzero. A POVM is informationally complete if its outcome probabilities are sufficient to determine uniquely the quantum state, i.e., to perform quantum state tomography~\cite{Prug77,Busch91}.  In other words, a bipartite state has zero discord from $B$ to $A$ if tomography on $A$ leaves the eigenstates of the density operator of $B$ unchanged.

Our method for verifying quantum discord is an improvement on the existing method~\cite{modi2011}, as it only requires measurement of one IC-POVM on $A$.
Hence, this method can be readily applied to continuous-variable systems where an IC-POVM is available from either heterodyne or homodyne measurements. We discuss in Sec.~\ref{Method} that the commutativity of the conditional states of $B$ can be efficiently tested by checking the commutation relations between one nondegenerate conditional state and all other states.

Quantum discord is defined as the difference between two classically equivalent measures for mutual information~\cite{Zurek01}.  According to Bayes's rule for classically correlated probability distributions, the quantities $I(A:B)=H(A)+H(B)-H(A,B)$, $J(A|B)=H(A)-H(A|B)$, and $J(B|A)=H(B)-H(B|A)$, where $H$ denotes the Shannon entropy and $H(A|B)=H(A,B)-H(B)$ is the conditional entropy, are all equal; they are called the classical mutual information.  For a bipartite quantum system, the quantum mutual information is defined, in analogy to $I(A:B)$, by $I(\rho_{AB})=S(\rho_A)+S(\rho_B)-S(\rho_{AB})$, where $S(\rho)=-\text{Tr}[\rho\log(\rho)]$ is the von Neumann entropy.  A measurement-based, quantum version of the conditional entropy is $S_{\{\Pi_j\}}(A|B)=\sum_{j}p_{j}S(\rho_{A|j})$, where $p_j=\text{Tr}[\rho_{AB}\Pi_{j}]$, $\rho_{A|j}=\text{Tr}_B[\rho_{AB}\Pi_j]/p_j$, and the set $\{\Pi_j\}$, with $\sum_j\Pi_j=\mathbb{I}$, makes up a POVM measurement on subsystem $B$.  This conditional entropy depends on the choice of measurement; hence, the quantum analogue of $J(A|B)$ is defined by minimizing over all possible measurements: $J^{\leftarrow}(\rho_{AB})=S(\rho_{A})-\text{inf}_{\{\Pi_j\}}S_{\{\Pi_j\}}(A|B)$.

The quantum discord from $B$ to $A$ is then defined as the difference between these two ways of defining a quantum mutual information:
\begin{align}
D^{\leftarrow}(\rho_{AB})&=I(\rho_{AB})-J^{\leftarrow}(\rho_{AB}) \nonumber \\
&=S(\rho_{B})-S(\rho_{AB})+\text{inf}_{\{\Pi_j\}}S_{\{\Pi_j\}}(A|B)\;.
\end{align}
The quantum discord is zero if and only if the quantum state can be expressed in the form
\begin{equation}
\rho_{AB}=\sum_{j} p_{j} \rho_{j}\otimes\ket{j}\bra{j\,}\;,
\label{form}
\end{equation}
where $\{\ket{j}\}$ are orthogonal states and $0\leq p_{j}\leq 1$~\cite{Zurek01,MLang11}.  For a quantum state with this form, local measurements on $B$ in the basis $\{\ket{j\,}\bra{j\,}\}$ leave the system unperturbed, and all the state information can be extracted without joint measurements.
Notice that the state~\eqref{form} is diagonal in a conditional product basis pointing from $B$ to $A$, i.e., an orthogonal basis of the form $\{|f_{jk}\rangle\otimes|e_j\rangle\}$, where the states $|f_{jk}\rangle$ are the eigenstates of $\rho_j$. Hence, for a given quantum state, quantum discord can be directly verified by diagonalizing the joint density operator~\cite{MLang11}.

Although there is no general method for minimizing over all possible measurements in order to calculate the conditional entropy of a state, this can sometimes be done when there are restrictions to certain classes of states and POVMs.  Thus the Gaussian quantum discord is defined as the quantum discord for two-mode Gaussian states where the evaluation of the conditional entropy is restricted to generalized Gaussian measurements~\cite{Datta10,Paris10}.

This paper is structured as follows.  In the next section, we introduce the method for verifying quantum discord.  In Sec.~\ref{CVsystems}, we discuss the application of the method to continuous-variable systems.  In Sec.~\ref{Gaussian}, we propose a technique for verifying quantum discord of Gaussian states. Based on that, we show that for Gaussian states, only product states have zero discord; hence states with nonzero Gaussian discord have nonvanishing quantum discord.

\section{Method}
\label{Method}

A previously proposed measurement-based method~\cite{modi2011} for verifying quantum discord is based on testing whether a quantum state can be expressed in the form~\eqref{form} of states with zero discord. We improve this previous proposal by showing that the quantum discord can be verified with only one IC-POVM.

\begin{theorem}
For a bipartite system $\rho_{AB}$, the necessary and sufficient condition for having zero discord
from $B$ to $A$, $D^{\leftarrow }(\rho_{AB})=0$, is that the states of subsystem $B$, $\rho_{B|k}=\text{Tr}_{A}[M_{k}\rho_{AB}]/\text{Tr}_{AB}[M_{k}\rho_{AB}]
=\text{Tr}_{A}[M_{k}\rho_{AB}]/p_k$, conditioned to the outcomes $k$ of an IC-POVM on $A$ (POVM elements $\{M_{k}\}$), commute with one another, i.e.,
\begin{equation}
\left[\rho_{B|k},\rho_{B|k'}\right]=0\,,\quad\text{for any}\ k\ \text{and}\ k'\,.
\label{condition}
\end{equation}
\end{theorem}
\begin{proof}
For a state having the zero-discord form~\eqref{form},
\begin{equation}
\rho_{B|k}=
\frac{\sum_jp_j\text{Tr}_A[M_k\rho_j]|j\,\rangle\langle j\,|}
{\sum_jp_j\text{Tr}_A[M_k\rho_j]}\;,
\end{equation}
immediately demonstrating that the states $\rho_{B|k}$ are all diagonal in the basis $\{|j\rangle\}$ and thus commute.

For the converse, we assume the condition~\eqref{condition}.  For such a set of commuting conditional states $\{\rho_{B|k}\}$, there exists an orthonormal basis, $\{\ket{j}\}$, that diagonalizes all the conditional states, $\rho_{B|k}=\sum_{j}\lambda_{kj}\ket{j\,}\bra{j\,}$.  To say that the POVM elements $\{M_k\}$ make up an IC-POVM is to say that they span the space of operators and thus there exist operators $\{N_k\}$ such that $\rho=\sum_k N_k\text{Tr}[M_k\rho]$ for any density operator $\rho$.  Applying this identity to the joint state gives
\begin{equation}
\rho_{AB}=\sum_k N_k\text{Tr}_A[M_k\rho_{AB}]=\sum_{k,j} \lambda_{kj}p_kN_k\otimes\ket{j}\bra{j}\;.
\end{equation}
Hence, $\rho_{AB}$ has the form~\eqref{form}, with $\rho_j=\sum_k\lambda{kj}p_kN_k$.
\end{proof}

Physically, what the proof says is that for a state of zero discord from $B$ to $A$, the measurement of the POVM $\{M_k\otimes\ket{j\,}\bra{j\,}\}$ extracts all information about the state $\rho_{AB}$.  From the perspective of the original definition of discord~\cite{Zurek01}, one imagines extracting this information by first measuring $B$ in the basis $\{|j\rangle\}$ and then measuring an IC-POVM on $A$.  Our criterion for zero discord works from the opposite perspective by reversing the order of the measurements on $A$ and $B$.

In order to test whether an unknown quantum state has nonzero discord experimentally, based on this theorem, one needs to measure an IC-POVM on subsystem $A$ and determine, by state tomography for each outcome, the corresponding states of the subsystem~$B$.
This procedure continues until one of the commutation relations between conditional states of subsystem $B$ is nonzero.
If subsystem $A$ has a $d$-dimensional Hilbert space, one can always find an IC-POVM that has $d^2$ POVM elements~\cite{Prug77,Caves02}.  Hence, there are $d^2$ conditional states of subsystem~$B$ and $d^2(d^2-1)/2$ commutation relations between all pairwise states.
However, as the conditional states are Hermitian operators, the most efficient way to check commutativity is to calculate the commutation relations between one of the states with no degeneracy and all other states. In this case, there are at most $d^2-1$ commutation relations to be checked.
Also, if some prior knowledge about the state in question is available, as is often the case in practice, quantum discord can be tested by considering only a few IC-POVM elements.
 Consider, for example, the maximally entangled state, $\ket{\psi}=\sum_{j=1}^d\ket{j\,}\ket{j\,}/\sqrt d$.  Any two rank-one outcomes, $\ket{n}$ and $\ket{\eta}$, on one of the subsystems, provided $0<|\braket{n}{\eta}|<1$ (these could be outcomes from two distinct, nonorthogonal projective measurements), yield conditional states of the other subsystem that do not commute.  Also, as we show below, for Gaussian states only two different heterodyne outcomes are sufficient to verify quantum discord.

\section{Continuous-variable systems}
\label{CVsystems}

An interesting feature of this method is that it can be readily applied to continuous-variable systems, as complete sets of IC-POVMs are available from heterodyne or homodyne detection. Two sets of measurements are required, one on each of the subsystems.  In general, one needs to do state tomography to construct the quasiprobability distributions of subsystem $B$ for all the states conditioned to outcomes of the measurement performed on subsystem~$A$.  Then the commutativity of the states $\rho_{B|k}$, which are represented in terms of quasiprobabilities, must be checked in order to verify discord. This can be efficiently done by finding one nondegenerate state and calculating the commutation relations between that state and all other states using an appropriate relation in terms of the reconstructed quasiprobabilities.  For instance, if the Wigner functions  $W_{B|k}(\alpha)$ of conditioned states of subsystem $B$ are available, the commutation relations between corresponding density operators can be calculated by using the Moyal Bracket~\cite{Moyal}:
\begin{align}
W_{kk'}(\alpha)=\frac{1}{2\pi}\int\text{d}^2\!\beta\,&\text{d}^2\!\beta'\,
W_{B|k}\bigl(\alpha+\textstyle{\frac{1}{2}}\beta\bigr)
W_{B|k'}\bigl(\alpha+\textstyle{\frac{1}{2}}\beta'\bigr) \nonumber \\
&\times \sin\!\left(i\frac{\beta\beta^{\prime*}-\beta'\beta^*}{2}\right)\;.
\end{align}
Here $W_{kk'}(\alpha)$ is the Wigner-like function for the operator $-i[\rho_{B|k},\rho_{B|k'}]$. If the states commute with each other then $W_{kk'}(\alpha)=0$ for all $\alpha$. Alternatively, the commutation relations can be calculated using characteristic functions,
\begin{align}
\mathlarger{\chi}_{kk'}(\xi)
=\frac{2}{\pi}\int\text{d}^2\zeta\,&
\mathlarger{\chi}_{B|k}\bigl({\textstyle\frac{1}{2}}\xi+\zeta\bigr)
\mathlarger{\chi}_{B|k'}\bigl({\textstyle\frac{1}{2}}\xi-\zeta\bigr)\nonumber \\ &\times\sin\!\left(i\frac{\xi\zeta^*-\xi^*\zeta}{2}\right)\;,
\end{align}
or in terms of any other quasiprobability distributions~\cite{Agarwal-Wolf}.

For states with zero discord, the eigenstates of the conditional density operator of $B$ do not change while $A$ is being fully determined from measurements of an IC-POVM. For continuous-variable systems, defined on an infinite-dimensional Hilbert space, the IC-POVM will have an infinite number of outcomes.  In practice, only a finite number of measurement outcomes can be explored.  For instance, in homodyne detection only a finite number of phases are considered, and the phase space is subdivided into a finite number of bins.  This introduces errors in the state estimation and uncertainties for the reconstructed quasiprobabilities of the conditional states, which propagate to the distributions representing the operators $-i[\rho_{B|k},\rho_{B|k'}]$. If one of these commutator distributions takes on a nonzero value at some point, which is larger than its associated uncertainty, then quantum discord is necessarily nonzero; otherwise, it is not clear whether the discord is nonzero.
However, by having some prior knowledge about the state, such as being Gaussian, the error can be estimated, and it can be made arbitrarily small using a sufficiently large number of measurements.

\section{Gaussian states}
\label{Gaussian}

A special class of continuous-variable states consists of the Gaussian states, i.e., those states whose Wigner function is a Gaussian function.  Such states are uniquely characterized by the means and covariance matrix of their quadrature components, $x$ and $p$.  For two systems, with modal annihilation operators $\hat a=x_1+ip_1$ and $\hat b=x_2+ip_2$, we define quadrature vectors for each system, $\mathbf{x}_1=(x_1,p_1)$ and $\mathbf{x}_2=(x_2,p_2)$, and we define an overall quadrature vector $\mathbf{x}=(\mathbf{x}_1,\mathbf{x_2})=(x_1,p_1,x_2,p_2)$.

The means of the quadrature components can be set to zero by locally displacing the two systems.  Then the state is specified by its covariance matrix~\cite{Adesso-Illuminati}
\begin{equation}
\bm{\sigma}=\langle\mathbf{x}^T\mathbf{x}\rangle=
\begin{pmatrix}
\mathbf{A} & \mathbf{C} \\
\mathbf{C}^{T} & \mathbf{B}
\end{pmatrix}
\;.
\label{sigmagen}
\end{equation}
Using local unitary operations that preserve the Gaussian form of the states, the covariance matrix of a bipartite Gaussian state can be brought to a standard form in which
$\mathbf{A}=\text{diag}(a,a)$, $\mathbf{B}=\text{diag}(b,b)$, and $\mathbf{C}=\text{diag}(c,d)$, where $a\ge0$ and $b\ge0$.  This can be accomplishing by first applying local unitary rotations that diagonalize $\mathbf{A}$ and $\mathbf{B}$, then using local squeezing operations to transform these diagonal blocks to $\mathbf{A}=\text{diag}(a,a)$ and $\mathbf{B}=\text{diag}(b,b)$, and finally applying further local unitary rotations to diagonalize $\mathbf{C}$. Notice that positivity of the density operator imposes the uncertainty-principle constraint~\cite{RSimon00},
\begin{align}
&\bm{\sigma}+\frac{i}{4}\bm{\Omega}\ge0\;,\quad
\Omega=\begin{pmatrix}
\mathbf{J} & \mathbf{0} \\
\mathbf{0} & \mathbf{J}
\end{pmatrix}\;,\quad
&\mathbf{J}=\begin{pmatrix}
0&1\\-1&0
\end{pmatrix}\;.
\end{align}
For a covariance matrix in standard form, this implies that $a^2\ge1/16$, $b^2\ge1/16$, $ab\ge c^2$, and $ab\ge d^2$, plus cubic and quartic constraints on $a$, $b$, $c$, and $d$.

It has been shown that Gaussian discord for a two-mode Gaussian state is zero if and only if $\mathbf{C}=0$~\cite{Datta10,Paris10}.  Here we show this condition is also necessary and sufficient for having zero discord.

A zero-mean Gaussian state with the standard form of the covariance matrix has characteristic function $\text{\large${\chi}$}(\mathbf{k})=\langle e^{i\mathbf{k}\mathbf{x}^T}\rangle=e^{-\mathbf{x}\bm{\sigma}\mathbf{x}^T/2}$
and Wigner function
\begin{align}
&W(\mathbf{x}_1,\mathbf{x}_2)=W(\mathbf{x})\nonumber\\
&=\frac{1}{4\pi^2\sqrt{\det\bm{\sigma}}}
\,\exp\!\left(-\frac{\mathbf{x}\bm{\sigma}^{-1}\mathbf{x}}{2}\right)\nonumber\\
&=\frac{1}{4\pi^2\sqrt{(ab-c^2)(ab-d^2)}}\nonumber\\
&\times\exp\!\!\left(-\frac{bx_1^2+ax_2^2-2cx_1x_2}{2(ab-c^2)}-\frac{bp_1^2+ap_2^2-2dp_1p_2}{2(ab-d^2)}\right)\;.
\end{align}

Suppose Alice makes a heterodyne measurement on subsystem~$A$; i.e., she uses the IC-POVM whose POVM elements are the coherent states $\ket{\beta}\!\bra{\beta}$.  Let $\beta=x'_1+ip'_1$ specify the outcomes of here measurement.  Then the state $\rho_{B|\mathbf{x}'_1}$, conditioned on these outcomes, has Wigner function
\begin{align}
W_{B|\mathbf{x'_1}}(\mathbf{x}_2)=\frac{1}{N}\int\text{d}^2\mathbf{x}_1\,
W(\mathbf{x}_1,\mathbf{x}_2)W_{\mathbf{x}'_1}(\mathbf{x}_1)\;,
\end{align}
where $W_{\mathbf{x}'_1}(\mathbf{x}_1)=2\exp[-2(x_1-x'_1)^2-2(p_1-p'_1)^2]/\pi$ is the Wigner function of coherent state $\ket{\beta}=\ket{x'_1+ip'_1}$ and $N$ is a normalization factor. Integration yields
\begin{align}
W_{B|\mathbf{x}'_1}(\mathbf{x}_2)=&\frac{1}{N'} \exp\!\left(-\frac{1}{2}\,x_2^2f(a,b,c)+x_2x'_1g(a,b,c)\right)\nonumber\\
&\times \exp\!\left(-\frac{1}{2}\,p_2^2f(a,b,d)+p_2p'_1g(a,b,d)\right)\;,
\label{WB}
\end{align}
where
\begin{align}
f(a,b,z)&=\frac{1}{ab-z^2}\left(a-\frac{z^2}{b+4(ab-z^2)}\right)\;, \nonumber \\
g(a,b,z)&=\frac{4z}{b+4(ab-z^2)}\;, \nonumber
\end{align}
and $N'$ is a normalization factor.

The peak of system~$B$'s conditional Wigner function~\eqref{WB}, located at
\begin{equation}
\gamma=\frac{g(a,b,c)}{f(a,b,c)}x'_1+i\frac{g(a,b,d)}{f(a,b,d)}p'_1\;,
\label{peak}
\end{equation}
depends on the outcomes of the measurement on $A$, $x'_1$ and $p'_1$, unless $c=0$ and $d=0$. Consequently, the eigenvectors of the conditional state $\rho_{B|\mathbf{x}'_1}$, which are generally displaced, squeezed number states, $\{D(\gamma)S(\zeta)\ket{n}\}$, displaced to the Wigner-function peak $\gamma$, change based on the outcomes of the heterodyne measurement performed on subsystem~$A$. This indicates nonzero discord, since the eigenvectors do not commute.  Therefore, without explicitly calculating any commutation relations, we can see that the bipartite Gaussian state has nonvanishing discord unless $c=0$ and $d=0$. Transforming back from the standard form to the general convariance matrix~\eqref{sigmagen}, one can say that a bipartite Gaussian state has zero discord if and only if $\mathbf{C}=0$, i.e., if and only if the state is a product state. This also implies that states with Gaussian discord ($\mathbf{C}\neq0$) have nonzero quantum discord.

These results show that quantum discord of Gaussian states can be verified using only two different heterodyne outcomes on one subsystem and finding (by tomography) the points in the phase space at which the corresponding conditional Wigner functions attain their maximum values. If those points do not coincide, the quantum discord is nonzero, since having different peaks guarantees that the corresponding eigenstates, which are displaced, squeezed number states, do not commute.
This argument can also be applied to nonGaussian states: if there are two conditional Wigner functions with the same shape, but located at different points in phase space, they correspond to states $\rho$ and $D(\nu)\rho D^{\dagger}(\nu)$, which have two different sets of eigenvectors $\{\ket{\psi_i}\}$ and $\{D(\nu)\ket{\psi_i}\}$, which is sufficient evidence that the quantum discord is nonzero.  Note that discord exists even if only one of $c$ and $d$ is nonzero, so to uncover discord of Gaussian states with only two heterodyne outcomes, one should choose the outcomes to be different for both quadratures.

\section{Conclusion}
We have introduced a method for verifying quantum discord of any bipartite quantum system.  The method is based on the fact that all the information in states with zero discord from subsystem $B$ to subsystem $A$ can be fully extracted by measurements that are diagonal in a single basis of $B$. In order to verify discord, one needs to perform an IC-POVM on subsystem~$A$ and check whether the conditional density operators of the subsystem~$B$ commute, i.e., whether they share the same eigenstates.

It is worth mentioning that, in practice, one would check commutativity of the conditional states as an IC-POVM is being performed on $A$. This can be efficiently done by finding a nondegenerate conditional state and calculate the commutation relations between that state and other states. In this case, the maximum number of commutation relations to be checked scales linearly with the number of IC-POVM elements.  Once one of these commutators is found to be nonzero, that confirms nonzero discord.  This method can be simply applied on continuous-variable systems by using homodyne or heterodyne detection and calculating the commutation relations in terms of quasiprobability distributions.  We have shown that a bipartite Gaussian state has nonzero quantum discord if and only if it is not a product state, which is the same as the condition for having nonzero Gaussian discord. Moreover, we show that with only two heterodyne outcomes and without calculating any commutation relations, quantum discord of Gaussian states can be verified.

\vspace{18pt}

\acknowledgments

We thank M.~Gu and N.~Walk for useful discussions.  This research was conducted by the Australian Research Council Centre of Excellence for Quantum Computation and Communication
Technology (Project number CE11000102) and was partially supported by U.S.\ National Science Foundation Grant Nos.~PHY-1005540, PHY-0903953, and PHY-1212445.


\section{Introduction}

The current paradigm,  which is the basis of the standard cosmological model, states that the large-scale universe is homogeneous and isotropic, i.e., there is no privileged place or direction in the universe  \citep{2002SPIE.4847...86M} in a way that the observed inhomogeneities are local in character and they should vanish at sufficiently large scales  \citep{peacock1999cosmological, longair2008galaxy}. This principle is based mainly on the observations of the cosmic microwave radiation whose isotropy is one part in a hundred thousand \citep{2005astro.ph..5185B},  and philosophical considerations: the observed universe should be the same for any observer, the Copernican Principle.

Despite the success of the physical models constructed from this cosmological principle, some researchers now are debating heavily whether or not the observations from more recent galaxy catalogs corroborate the assumption of homogeneity and isotropy at scales greater than 60 Mpc/h \citep{2009MNRAS.399L.128S, 2009MPLA...24.1743C},  or if instead the galaxies are grouped into highly structured hierarchical pattern, with properties of fractality \citep{citeulike:2502783, 2011ARep...55..324V}.

In the last decade, the interest in the fractal analysis of the galaxy distribution has yielded promising results, becoming a fertile field for cosmological physics research \citep{2005astro.ph..5185B}. In particular, we do not know if the cold dark matter distribution have a fractal behaviour in the same way that is observed in the galaxy clustering. Some authors propose the same fractal dimension for the galaxy distribution and for the dark matter, where galaxies act as tracers of dark matter \citep{2008pc2..conf...60B}, while other researchers assign different fractal dimension; D=3 for dark matter, i.e. homogeneity in accord with cosmological principle and D=2 for galaxy clustering \citep{1998A&A...339L..85D}. Authors like \citet{0295-5075-71-2-332} suggest the use of multi-fractal models in order to analyse the dark matter haloes in numerical N--body simulations. He found that haloes of similar mass have a fractal distribution with a given dimension that grows as the mass diminishes.

The homogeneous self-gravity matter can develop fractal structure \citep{1998CeMDA..72...91C}. For this reason the N--body cosmological simulations can be used like a numerical laboratory where the fractal analysis of matter under gravitational interaction can be performed. Besides, in simulations with sufficient number of particles it is possible to detect the fractal behaviour with the scale, including the transition to homogeneity. The statistical algorithms used in the analyses of the simulation data can be applied subsequently to galaxies surveys \citep{citeulike:2502783}.

In this paper we present the calculations of mass-radius fractal dimension, mono-fractal lacunarity, multi-fractality and lacunarity spectrum of dark matter haloes distribution from the Millennium Simulation Project. We propose the use of multi-fractal dimension with the complementary measure of the lacunarity spectrum as a method for determining the homogeneity scale, i.e., the distance beyond which the transition to homogeneous regime should occur in a fractal matter clustering.

The paper is divided as follows: Section 2 presents the mathematical foundations involved in the fractal approach to the clustering of large scale matter; section 3 summarises the characteristics of the Millennium simulation, focusing on the cosmological parameters and scale used as the basis of the simulation; section 4 discusses the calculation of mass-radius fractal dimension and lacunarity as functions of radial distance in the clustering of dark matter haloes from the Millennium Simulation; section 5 deals with the calculation of  the multi-fractal spectrum and the spectrum of lacunarity on the same data from the Millennium Simulation; in Section 6 the discussion of the fractal results is presented and finally the conclusions are given in section 7.

\section{Theoretical development}

\subsection{Dimension Concept}

The most general definition for the dimension measure was developed by Felix Hausdorff in 1918 \citep{springerlink:10.1007/BF01457179},  looking to define the dimension concept for any metric space. He first considered the number $N(\epsilon )$ of balls of radius $\epsilon$ needed to cover a subset $A$ completely. The  d-dimensional measure is proportional to the limit:

\begin{equation}
 h^{d}(A)= \lim_{\epsilon \to 0}{\mathit{}N(\epsilon)\epsilon^{d}}
\end{equation}

Hausdorff proved that there is only one value $d$ for which this is different from zero or infinity. For this value is satisfied that:
\begin{equation}
 h^{d}(A)=\left\{\begin{matrix}\infty & \mathrm{ if }\ d<D_{H}(A){}\\
                0\ & \mathrm{ if }\ d>D_{H}(A){}\end{matrix}\right.
\end{equation} where $D_H$ is by definition the Hausdorff dimension for the subset $A$,  with $D_H$  the value for which $N(\epsilon )\propto 1/\epsilon ^{d}$.  In order to decide whether a set is fractal or not it is necessary to compare the above definition with the topological dimension,  defined in a simple way as the number of independent directions in which one can move around a given point of the set. In this form Mandelbrot define a fractal as a set for which the Hausdorff dimension strictly exceeds its topological dimension. \citep{mandelbrot1983fractal}

It is not easy to determine the Hausdorff dimension for galaxy clustering, although there are methods developed to apply the Hausdorff dimension, like the Minimal Spanning Tree \citep{martinez2002statistics}. For this reason it is necessary to use alternative concepts of fractal dimension for the analysis of the matter distribution on a large scale in the universe. One of the most used concepts is the mass-radius fractal dimension. According to Blumenfeld and Mandelbrot  \citep{PhysRevE.56.112},  given a sphere of radius $r$ in a space with Euclidean dimension $d$, which encloses a self-similar fractal structure, the total mass measure $M(r)$ enclosed by the sphere take the form:

\begin{equation}\label{Mass-radius}
M(r) = F{r}^{D_{m}}
\end{equation} here $M(r)$ represents the number of galaxies (or in our case dark matter haloes)  in a sphere of radius $r$, $D_{m}$ is the fractal mass-radius dimension, and $F$ is a constant related to the average distance between nearest neighbours. The exponent $D_{m}$, the fractal dimension is smaller than the embedding Euclidean dimension $d$ for inhomogeneity condition in the set.

Although the mass-radius dimension has the advantage of extending the measure of fractal dimension of the set on large scales, its accuracy is limited as it does not move the centre of observation, skewing the statistical analysis. Therefore it is necessary to refine the determination of fractality including average looking at several points in the sample as centres. According to \citet{citeulike:2502783} it is possible to average the mass-radius fractal dimension moving on different occupied points  of the galaxy distribution. The average mass-radius relation is a fractal measure which reflects the conditional cosmological principle: The cosmos appears statistically the same to all observers situated on a galaxy (point of a fractal) but not in a region of void \citep{mandelbrot1983fractal}.

To fully describe a fractal it is necessary to specify how the pre-factor $F$ behaves. For this reason it is necessary to use a complementary fractal quantity, the lacunarity, which is a general measure for the analysis of spatial patterns.

\subsection{Lacunarity Definition}

The majority of the studies aiming to determine the fractality of large-scale structures in the universe begin with the calculation of the fractal dimension, but fractals with similar dimensions can occupy the space where they are embedded in different ways. Therefore, it is necessary to recognise that the fractal dimension is not enough to uniquely characterise fractal sets,  and for this reason \citet{PhysRevE.56.112} proposed a complementary fractal measure, the lacunarity.

The lacunarity appears in the context of the mass-radius fractal dimension, although we can also extend this concept for other fractal dimension definitions like the correlation dimension \citep{2002SPIE.4847...86M}. Lacunarity high values suggest the presence of large empty regions inside the clustering and therefore an increased heterogeneity, while fractals with low values indicate approximation to homogeneity.

The starting point for calculating the lacunarity is the pre-factor $F$  from the mass-radius relation,  Eq. (\ref{Mass-radius}). There is not a unique definition to describe the properties of the factor $F$,  but Mandelbrot  proposed a description based in series of factors of variability \citep{PhysRevE.56.112},

\begin{equation}
S_{k}=\frac{C_{k}(D)}{[C_{1}(D)]^{k}},
\end{equation} where $C_k(D)$ is the $kth$ cumulant, a statistical measure of the correlation designed to go to zero whenever any one or more quantities under study become statistically independent. The cumulants are related with the moments, in particular $k_1 = \mu_1$ is the mean, $k_2$ is the variance, and $k_3= \bigl\langle(X - \mu_1)^3\bigr\rangle$ for a variable $X$ under study. The simplest of this variability factors, is the second order variability factor $S_{2}$ through which the lacunarity $\Phi$ is defined as:

\begin{equation}
\Phi =\frac{\bigl\langle(F - \bigl\langle F \bigr\rangle)^{2}\bigr\rangle}{\ \bigl\langle F \bigr\rangle^{2}}=\frac{\bigl\langle F^{2}\bigr\rangle}{\bigl\langle F\bigr\rangle^{2}}-1.
\end{equation}

The Lacunarity analysis complements the studies of fractal sets enabling us to distinguish between similar fractal patterns which occupy the space in a different form. \citet{martinez2002statistics} applied the above definition of lacunarity to the Las Campanas survey data and the corresponding fractal Levy flight model. They found great difference in lacunarity behaviour with similar fractal mass-radius dimension.

The average mass-radius dimension and the lacunarity are sufficient  to characterise sets where the fractal calculation of matter clustering at different scales and levels of mass density are the same. In the case of scale dependency,  we must go beyond in the characterising of the large scale matter distribution. This kind of analysis can be performed by the multi-fractal formalism.

\subsection{Generalised Dimension}

The dimension definitions presented above, represent particular cases of the multi-fractal spectrum of generalised dimension, method applied by the first time to the large scale structure in the universe by \citet{1992PhyA..185...45C}.  The large scale matter distribution can be characterised with a fractal spectrum defined from a generalised correlation. First, is necessary to introduce the correlation integral $C_{2}$, a function capable to measure the number of neighbours that on average a chosen centre has within a distance  $r$. We outline the development in accordance with the notation described in \citet{2008MNRAS.390..829B}:

\begin{equation}\label{C2}
 C_{2}(r)=\frac{1}{\mathit{NM}}\sum_{i=1}^M n_{i}(r),
\end{equation} where $N$ is the total number of particles inside the distribution, $M$ is the number of particles used as centres and the summation is performed over the set of chosen centres. $n_{i}(r)$ is the number of particles within a radial distance $r$ from a particle at the point $i$, defined as:
\begin{equation}\label{ni}
 n_{i}(r)=\sum_{j=1}^N \Theta(r-\mid \boldsymbol{x_i-x_j} \mid),
\end{equation}
where the summation is performed over all the particles in the sample. The coordinates of each particle in our space of three dimensions are denoted by $\boldsymbol{x_j}$, and $\Theta$ is the Heaviside function, defined such that $\Theta(x) = 0$ for $x < 0$ and $\Theta(x)=1$ for $x \geqslant 0$. The number of particles around each centre $n_{i}(r)$ is determined by counting the number of particles around that centre that lie inside a comoving sphere of radius $r$ from it.

From Equation (\ref {C2}) the correlation dimension is defined similarly to the definition of mass-radius dimension:
\begin{equation}\label{D2}
 C_{2}(r) \approx  r^{D_{2}}
\end{equation}

So that the correlation dimension is calculated as a derivative:

\begin{equation}
  D_{2}=\frac{d\log C_{2}(r)}{d\log r}
\end{equation}

From the correlation integral $C_{2}(r)$, Equation (\ref{C2}), the generalised correlation integral can be defined as:

 \begin{equation}
 C_{q}(r)= \frac{1}{\mathit{NM}}\sum_{i=1}^M [n_{i}(r)]^{q-1}
\end{equation}
where $M$ is the number of centres, $N$ is the total number of particles included in the sample, $n_{i}(r)$ is the same expression defined in the Equation (\ref {ni}) and $q$ is called the structure parameter, which corresponds to an arbitrary real number. From this generalised correlation it is possible to do an expansion in powers of $log(r)$ as described by \citet{1997Chaos...7...82P} and to calculate directly the multi-fractal dimension and the lacunarity spectrum:

 \begin{equation}
log \left[ C_{q}(r)^{1/(q-1)}\right] =D_{q}log(r)+log(F_{q})+O \left( \frac{1}{log(r)}\right).
\end{equation}
Keeping only the first two terms on the right side, we have the relation between the generalised correlation integral and the generalised fractal dimension:

\begin{equation}
C_{q}(r)^{1/(q-1)}=F_{q} r^ {D_{q}}
\end{equation}

In this manner the generalised dimension and the generalised lacunarity can be defined in the same way as for the mass-radius fractal dimension, the generalised dimension $D_{q}$:

\begin{equation}
  D_{q}=\frac {1}{(q-1)} \frac{d\log C_{q}(r)}{d\log r},
\end{equation} and the corresponding generalised lacunarity from the pre-factor  $F_q$ for every structure parameter $q$:

\begin{equation}
\Phi _{q}=\frac{\bigl\langle(F_{q} - \bigl\langle F_{q} \bigr\rangle)^{2}\bigr\rangle}{\ \bigl\langle F_{q} \bigr\rangle^{2}}=\frac{\bigl\langle F_{q}^{2}\bigr\rangle}{\bigl\langle F_{q} \bigr\rangle^{2}}-1.
\end{equation}

If for any $q_1\neq q_2$,  $D_{q_1} = D_{q_2}$ is verified, it is said that the distribution is a homogeneous fractal (mono-fractal). For $q \geqslant 1$, $D_q$  explores the scaling behaviour  in high density environments (clusters and superclusters) and  for $q < 1$ values, $D_q$ explores the scaling behaviour in low-density environments, i.e., voids \citep{2009MNRAS.399L.128S}. In the event that the distribution of dark matter haloes undergo the transition to homogeneity, all values of  the fractal dimension must close to the physical space dimension $(D_q \to 3) $ and the spectrum of lacunarity must tend to zero $(\Phi _{q} \to 0)$, at the same radial distance $r$.

Next we will apply the concepts of average mass-radius dimension,  generalised dimension and generalised lacunarity to the dark matter haloes spatial distribution generated by the Millennium dark matter simulation.

\section{Millennium Simulation}

The Millennium Simulation \citep{2005Natur.435..629S} is one of the most important computational efforts in contemporary cosmology. It uses the most accepted parameters inside the standard cosmology, with particles evolving since redshift $z=127$ until the present $z=0$. Although the standard cosmological model contains two ingredients which have not yet been verified by laboratory experiments (the components of dark matter and dark energy), the $\Lambda$CDM model is almost universally accepted by cosmologists as the best description of the present observational data. The spatial geometry is very close to flat and the initial perturbations which should have given origin to galaxies are Gaussian, adiabatic, and nearly scale-invariant. The cosmological parameters assumed inside this computational algorithm are summarised in Table \ref {table: Millparameters}.

\begin{table}
\centering
\begin{tabular}{llllll}
 \hline
$\Omega_m$ & $\Omega_b$ & $\Omega_\Lambda$ & $h$ & $\sigma_8$ & $n$ \\
\hline
$0.25$ & $0.045$ & $0.75$ & $0.73$ & $0.9$ & $1$
\end{tabular}
\caption{Cosmological parameters for the Millennium Simulation}
\label{table: Millparameters}
\end{table}

Here the density parameters $\Omega_i = \rho_i /  \rho_{crit}$ are defined from the critical density $ \rho_{crit}=3H_0^2/8\pi G $. The total present matter density composed by dark matter and baryonic matter $\Omega_m= \Omega_b+\Omega_{dm}$, the dark energy density $\Omega_\Lambda$ with $\Omega_m+ \Omega_\Lambda=1$, the Hubble parameter $h=H_0/(100km s^{-1}Mpc^{-1})$, $\sigma_8$ the linear-theory power spectrum variance in spheres of radius 8 Mpc/h and the spectral index $n$ used to describe the density perturbations in Fourier space, where $n=1$ corresponds to Harrisson-Zeldovich scale invariant spectrum.

The simulation is implemented with the length of the side of the simulation cube $L_{box} = 500 Mpc/h$, the number of particles $N_{part}=2160^3\approx  10^{10}$ , with the resulting particle mass $m_p=8.6 \times 10^8 M_{\odot}/h$  and the gravitational softening length $\eta=5 kpc/h$. The Virgo consortium, provides free access to the database, which authorised users can entry using the Structured Query Language (SQL). It is expected that different researchers employ the information to explore how the structures of dark matter and galaxies evolve in the standard cosmology.

We use the information provided by this simulation in order to perform a fractal analysis of the dark matter haloes clustering, aiming to find the way in which the dark matter haloes are grouped as a function of  the radial distance and to detect the scale of transition to homogeneity.

\section{Fractal Mass-Radius Dimension and Lacunarity of Millennium Simulation data}

The first step in fractal analysis of dark matter haloes clustering, for $z = 0$ from the millennium simulation, is the choice of a random sample of centres in the entire simulation ($\approx 1.5 \times 10^7$ dark matter haloes).  We chose randomly one hundred centres with the purpose to determine the functional relationship between the logarithm of the number of  dark matter haloes and the logarithm of the radial distance. We want to know if there is a linear function that leads to the determination of the fractal correlation dimension $D_m$ and the pre-factor $F$ as shown in Figure \ref{fig: lognlogr}.

\begin{figure}
\includegraphics[width=1 \linewidth,clip]{lognvslogr.pdf}
 \caption{Millennium Simulation Fractal Analysis:  Number of dark matter
 haloes \textit{vs} Radial distance  over 100 centres randomly chosen, in logarithmic scale. The initial data dispersion and the slope growth is shown. }
	\label{fig: lognlogr}
\end{figure}

In the case of a mono-fractal set, the expected relationship between the logarithm of the number of particles and the logarithm of radial distance should be close to a linear behaviour. By observing Figure \ref{fig: lognlogr} a large scatter in the relationship for distances $r < 20$ Mpc/h is evident. Also, it would be inappropriate to assign a single slope for the entire curve when a gradual increase of slope with the radial distance is observed, suggesting a fractal dimension growing with scale. It is therefore necessary to calculate the fractal dimension and the pre-factor $F$ in small sections of the curve for each centre, in order to determine if there is a radial distance dependency of the fractal quantities.

Dividing the curve in segments using the sliding window technique  \citep{ 2002SPIE.4847...86M, 2005PhRvE..72a6707R} and making a linear least-square fit through each set of  successive points in the log-log plot, it is possible to calculate the slope and intercept of the line that approximates the curve section over every point for each of the curves generated from each centre, and then determine the average over a thousand centres (more representative sample) with their respective standard deviation. We also analyse the behaviour of this average mass-radius fractal dimension as the sample size approaches to the limits of the simulation, as shown in Figure \ref{dmvsr}.

\begin{figure*}
	\includegraphics[width=0.45\linewidth,clip]{FDAE140.pdf}
	\includegraphics[width=0.45\linewidth,clip]{FDAE160.pdf}
         \includegraphics[width=0.45\linewidth,clip]{FDAE180.pdf}
         \includegraphics[width=0.45\linewidth,clip]{FDAE200.pdf}

         \caption{The average of mass-radius fractal dimension \textit{vs} radial distance for the dark matter haloes distribution in the Millennium simulation. Spurious homogenisation effects in the curves around 20 Mpc/h are observed to fractal calculations for depths larger than 180 Mpc/h.  $\pm1\sigma$ error bars are shown.}
	\label{dmvsr}
\end{figure*}

In Figure \ref{dmvsr}, the graphs for different depths have a similar behaviour; a rapid increase in the fractal dimension to radial distances about 20 Mpc/h followed by slower growth until it reaches the physical dimension of space to depths beyond the 100 Mpc/h. Those graphics with radial distances from the centres exceeding 180 Mpc/h, where the dark matter haloes are close to the edges of the simulation,  show that the average mass-radius fractal dimension at short distances is larger than the physical dimension of the space where the whole fractal is embedded, a physically unacceptable situation.

The behaviour of the pre-factor $F$ is shown in Figure \ref{fvsr}. The values of the pre-factor change with scale, starting with high values near the centres and diminishing to reach a low constant value at distances larger than 100 Mpc/h. Therefore, the average distance between dark matter neighbours decreases with radial distance reaching a constant value for radial distances larger than 100 Mpc/h. A change at short distances of the behaviour of the curves for depths from the centres larger than 180 Mpc/h is detected.

\begin{figure*}
	\includegraphics[width=0.45\linewidth,clip]{F140.pdf}
	\includegraphics[width=0.45\linewidth,clip]{F160.pdf}
         \includegraphics[width=0.45\linewidth,clip]{F180.pdf}
         \includegraphics[width=0.45\linewidth,clip]{F200.pdf}

         \caption{Fractal Pre-factor $F$ \textit{vs} radial distance for the dark matter haloes distribution in the Millennium simulation. There is a tendency to achieve a constant value with increasing radial distance. $\pm1\sigma$ error bars are shown.}
	\label{fvsr}
\end{figure*}

This results confirm that there are spurious homogenisation effects in the calculations of the fractal dimension. Thus, it is necessary to calculate the fractal dimension over radial distances smaller than 180 Mpc/h, i.e., radial distances where the fractal dimension does not exceed $3\pm1\sigma$, as recommended by \citet{1989A&A...219....1W}. In our case the maximum radial distance that satisfies the above requirement is 160 Mpc/h.

To complement the vision of the dark matter haloes fractal distribution we calculate the lacunarity up to 160 Mpc/h radial distances, taking into account the pre-factor of the sliding window least-square fit analysis to the same depth within the simulation. The lacunarity of the dark matter haloes distribution, Figure \ref{resume160}, shows a maximum value near each chosen centre, followed by oscillatory lacunarity behaviour at radial distances smaller than 60 Mpc/h. In addition, above 60 Mpc/h we see a smoothly decreasing lacunarity function with tendency to homogeneity at radial distances $\approx 120$ Mpc/h where lacunarity values approach to zero.

\begin{figure*}
	\includegraphics[width=0.45\linewidth,clip]{Dimension160.pdf}
	\includegraphics[width=0.45\linewidth,clip]{Chisquare160.pdf}
	\includegraphics[width=0.45\linewidth,clip]{Prefactor160.pdf}
	\includegraphics[width=0.45\linewidth,clip]{Lacunarity160.pdf}

 \caption{Mass-radius quantities calculated with the sliding window technique: The average of the mass-radius dimension $D_m$, the pre-factor $F$, the lacunarity $\Phi$ and the $\chi^2$ test \textit{vs} the radial distance, for the dark matter haloes distribution. $\pm1\sigma$ error bars are shown. The average $<\chi^2> = 7.9\times10^{-3}$, with $\chi^2 _{max}=5.1 \times10^{-2}$.}
	\label{resume160}
\end{figure*}

The fractal quantities show a strong reliance on the scale. The decrease in the lacunarity curve and the scale dependence of the average mass-radius fractal dimension with the radial distance suggest us that the clustering of dark matter haloes do not behave like a mono-fractal. Therefore it is necessary to study the dark matter haloes distribution from the multi-fractal viewpoint.

\section{Multi-fractality and Multi-lacunarity of Millennium Simulation data}

In the last section we have shown that the Millennium Simulation dark matter haloes average mass-radius fractal dimension have a radial distance dependency with a lacunarity descending to zero around 120 Mpc/h. Now, the results of the generalised correlation integral for twelve different values of $q$ are presented in Figure \ref{Cqr}. This calculations include overdense regions $(q \geqslant 1)$ and low density regions $(q < 1)$. For $q =1$ the numerical limit of the generalised correlation integral is shown. To avoid the problems of spurious homogenisation found in the mass fractal dimension,  the multi-fractal calculus is limited to depths still far from the borders of the simulation. Based on the average mass-radius fractal behaviour, the maximum depth for the calculations in the Millennium data is 160 Mpc/h.

\begin{figure*}
	\includegraphics[width=0.45\linewidth,clip]{Cq_clusters.pdf}
	\includegraphics[width=0.45\linewidth,clip]{Cq_voids.pdf}

 \caption{Generalised Correlation Integral $C_q(r)$ for all the values of $q$ studied in this paper for the Millennium simulation dark matter haloes. Over densities  $q \geqslant 1 $ left, low densities $q < 1$ right. All graphics in log-log scale}
	\label{Cqr}
\end{figure*}

From the generalised correlation integral, the sliding window technique is used to determine the fractal dimension spectrum. We find the multi-fractal dimension behaviour as a function of the radial distance for the structure parameter values $q < 1$ as shown in the Figure \ref{Dqr1a}. Besides, we calculate the corresponding $\chi^2$ test in each multi-fractal relation over all distance range where we made the calculations; the low $\chi^2$ values provide us confidence on the used method. The $\chi^2$ average and maximum values for each structure parameter are shown in Table \ref{table:chi2v}.

\begin{figure*}
	\includegraphics[width=0.45\linewidth,clip]{D-6.pdf}
	\includegraphics[width=0.45\linewidth,clip]{D-5.pdf}
	\includegraphics[width=0.45\linewidth,clip]{D-4.pdf}
         	\includegraphics[width=0.45\linewidth,clip]{D-3.pdf}
	\includegraphics[width=0.45\linewidth,clip]{D-2.pdf}
 	\includegraphics[width=0.45\linewidth,clip]{D-1.pdf}
         \includegraphics[width=0.45\linewidth,clip]{D0.pdf}

	\caption{Multi-fractal Dimension Spectrum  $D_q(r)$ \textit{vs} radial distance $r$ for low density environments $-6 \leqslant q  \leqslant 0$ from the Millennium dark matter haloes. $\pm1\sigma$ error bars in the graphics, with 0.99 confidence level in t-student test.}
	\label{Dqr1a}
\end{figure*}

For structure parameter values $q \geqslant 1$ the radial distance dependency is presented in Figure \ref{Dqr2a}. In this manner we covered the range of the structure parameter $-6 \leqslant q \leqslant 6$ in steps of one. The  $\chi^2$ test was calculated in the same manner as for the multi-fractal relations with $q < 1$. There are low $\chi^2$ values over all radial distances, See Table \ref{table:chi2v}.

\begin{figure*}
	\includegraphics[width=0.45\linewidth,clip]{D1.pdf}
	\includegraphics[width=0.45\linewidth,clip]{D2.pdf}
	\includegraphics[width=0.45\linewidth,clip]{D3.pdf}
	\includegraphics[width=0.45\linewidth,clip]{D4.pdf}
	\includegraphics[width=0.45\linewidth,clip]{D5.pdf}
	\includegraphics[width=0.45\linewidth,clip]{D6.pdf}

	\caption{Multi-fractal dimension spectrum  $D_q(r)$ \textit{vs} radial distance $r$ for high density environments $ 1\leqslant q \leqslant 6$  from the Millennium dark matter haloes. $\pm1\sigma$ error bars in the graphics, with 0.99 confidence level in t-student test.}

\label{Dqr2a}
\end{figure*}

\begin{table*}
         	\centering
	\begin{tabular}{ccccccccccccccccc}
 	\hline
	$q$ & $-6$ & $-5$ & $-4$ & $-3$ & $-2$ & $-1$ & $0$\\
	\hline
	$<\chi^2>$ & $1.9\times10^{-3}$ & $8.8\times10^{-4}$ & $2.5\times10^{-3}$ & $1.2\times10^{-3}$ & $1.5\times10^{-3}$ & $3.1\times10^{-4}$ & $2.5\times10^{-6}$\\
	\hline
	$\chi^2 _{max}$ & $4.5\times10^{-1}$ & $1.0\times10^{-1}$ & $2.9\times10^{-1}$ & $6.7\times10^{-2}$ & $2.1\times10^{-1}$ & $7.1\times10^{-2}$ & $1.3\times10^{-4}$\\
	\hline
	$q$ & $1$ & $2$ & $3$ & $4$ & $5$ & $6$ \\
	\hline
	$<\chi^2>$ & $1.6\times10^{-6}$ & $1.6\times10^{-6}$ & $2.6\times10^{-6}$ & $3.1\times10^{-6}$ & $3.5\times10^{-6}$ & $4.7\times10^{-6}$\\
	\hline
	$\chi^2 _{max}$ & $9.0\times10^{-5}$ & $2.8\times10^{-5}$ & $2.8\times10^{-5}$ & $4.7\times10^{-5}$ & $6.9\times10^{-5}$ & $9.4\times10^{-5}$\\
	\hline
	\end{tabular}

\caption{$\chi^2$ mean and maximum values obtained for each multi-fractal dimension $D_q$. Top: low density environments $-6 \leqslant q  \leqslant 0$, bottom:  high density environments $ 1\leqslant q \leqslant 6$. For all the $D_q$ analysed  the $\chi^2$ values show a peak at scales smaller than $\approx 15$ Mpc/h and very low, almost constant, values at larger scales.}
\label{table:chi2v}
\end{table*}

In order to detect the transition to homogeneity in terms of the fractal dimension, we use the relationship:
\begin{equation}
 \mid D_q-3 \mid / 3 < \epsilon,
\end{equation}
where $\epsilon$ is a fiducial value which is calculated from the average of the relative dispersion $\sigma/D_q$.  In Table \ref{table:Rh} we present the shortest radial distances for which this condition for the transition to homogeneity is fulfilled.

\begin{table*}
         	\centering
	\begin{tabular}{ccccccccccccccccc}
	\hline
	$q$ & $1$ & $2$ & $3$ & $4$ & $5$ & $6$ \\
	\hline
	$\epsilon$ & $1.8 \times10^{-3} $ & $8.4 \times10^{-4}$ & $1.2 \times10^{-3}$ & $1.3 \times10^{-3} $ & $1.8 \times10^{-3}$ & $ 1.7 \times10^{-3}$ \\
	\hline
	$R_{h}[Mpc/h]$ & $101.1$ & $108.9$ & $107.2$ & $108.7$ & $117.7$ & $122.5$ \\
	\hline
	\end{tabular}

\caption{Scale of homogeneity transition $R_{h}$  for every fractal dimension in overdenses environments $1\leqslant q \leqslant 6$. The $\epsilon$ value is calculated from the relative dispersion of the fractal dimension  }
\label{table:Rh}
\end{table*}

This results are confirmed by the relationship between the multi-fractal dimension and the spectral parameter $q$, as shown in Figure \ref{Dqrq}.

\begin{figure*}
	\includegraphics[width=0.45\linewidth,clip]{Dq_70r.pdf}
	\includegraphics[width=0.45\linewidth,clip]{Dq_80r.pdf}
	\includegraphics[width=0.45\linewidth,clip]{Dq_90r.pdf}
         \includegraphics[width=0.45\linewidth,clip]{Dq_100r.pdf}
          \includegraphics[width=0.45\linewidth,clip]{Dq_110r.pdf}
         \includegraphics[width=0.45\linewidth,clip]{Dq_120r.pdf}
         \includegraphics[width=0.45\linewidth,clip]{Dq_130r.pdf}
         \includegraphics[width=0.45\linewidth,clip]{Dq_140r.pdf}
         \caption{Spectrum of the Multi-fractal dimension $D_q(r)$ as a function of structure parameter $q$ in the range $[-6,6]$ for Millennium simulation dark matter haloes. Homogeneity transition around 120 Mpc/h is determined. The solid line is an interpolation by the B\'ezier method.}
	\label{Dqrq}
\end{figure*}

Another possible definition for transition to homogeneity is:

\begin{equation}
{2{[C_2(r_1)/V(r_1)-C_2(r_2)/V(r_2)]} \over {[C_2(r_1)/V(r_1)+C_2(r_2)/V(r_2)]} }< \epsilon,
\end{equation}
$\forall  r_1,r_2 > R_h$, where $R_h$ is the homogeneity scale, and $\epsilon$ is a fiducial value. If  the same value for epsilon as above is used, then the transition to homogeneity appears at scales near 20 Mpc/h, which do not agree with the scales of transition to homogeneity we found from the multi-fractal spectrum or from the minimal values of the lacunarity spectrum.  Therefore, either this definition for transition to homogeneity is less appropriate than the former one or the value of epsilon should be much more smaller than $10^{-3}$ . In the latter case, it is not clear what fiducial value should be assigned to $\epsilon$ and how to estimate it.

To corroborate this result and because the fractal dimension does not indicate in what form the dark matter haloes set is filling the whole space, the analysis is completed with the calculation of the spectrum of lacunarity for overdenses regions, as shown below in Figure \ref{Phiqr}.

\begin{figure*}
	\includegraphics[width=0.45\linewidth,clip]{L1.pdf}
	\includegraphics[width=0.45\linewidth,clip]{L2.pdf}
	\includegraphics[width=0.45\linewidth,clip]{L3.pdf}
	\includegraphics[width=0.45\linewidth,clip]{L4.pdf}
         \includegraphics[width=0.45\linewidth,clip]{L5.pdf}
         \includegraphics[width=0.45\linewidth,clip]{L6.pdf}
         \caption{Multi-lacunarity spectrum $\Phi_q(r)$ in logscale \textit{vs}  radial distance $r$ for high density environments $q \geqslant 1$  from the Millennium dark matter haloes distribution. Oscillatory behaviour modulated by  a decreasing function. $ \pm1\sigma$ error bars are shown, with 0.99 confidence level in t-student test}
	\label{Phiqr}
\end{figure*}

\section{Discussion}

We confirm the dependence of fractal dimension with the radial distance; The most of the previous works reported only a single fractal dimension value to frame the clustering behaviour of large-scale mass \citep{1998A&A...339L..85D, 1999dmap.conf..336S, 2005A&A...443...11J}. The majority of recent works focus on the radial distance dependency of the fractal dimension and the multi-fractal analysis of large scale clustering in the Universe, for example: \citet{2002PhyD..168..404R} describe the scaling properties of large-scale structures in the Universe and the transition to homogeneity around 500 Mpc/h,  \citet{2005BASI...33....1S} found that the universe is homogeneous over scales larger than about 80 to 100Mpc/h,  \citet{2008MNRAS.390..829B} observed that the fractal dimension makes a rapid transition to values close to the physical space dimension at scales between 40 and 100 Mpc, while \citet{2009MNRAS.399L.128S} determinate that the Millennium simulation exhibit a transition to homogeneity at around 70 Mpc/h. In our case, the analysis of The Millennium Simulation confirms fractal behaviour (in the average mass-radius relation) at scales smaller than at least 100 Mpc/h with homogeneity transition around 110 Mpc/h for the dark matter haloes distribution.

On the other hand, we observe limitations in the application of  the average mass-radius fractal dimension and the multi-fractal dimension spectrum to radial distances larger than 180 Mpc/h. We detect spurious homogenisation effects from these calculations as the radial distance from the centres approaches the edges of the simulation. This effect is reported in the literature \citep{1989A&A...219....1W, 1996PhyA..226..195S}. As stated above, due to  this edge effect it is necessary to limit the maximum radius to which the fractal dimension can be calculated within the simulation. In our case the average mass-radius dimension, the lacunarity, the multi-fractal dimension and the lacunarity spectrum have to be limited to maximum radial distances of 160 Mpc/h, where the fractal dimension does not exceed the value of $3 \pm1\sigma$.

The behaviour of mass-radius fractal dimension can be divided in two regions clearly observed: a fast-growing region from each centre up to $\approx$ 20 Mpc/h in radial distance, followed by a slow-growing region to achieve homogeneity at radial distances around 100 Mpc/h. Above this distance the search for the scale of transition to homogeneity using the lacunarity and the multi-fractal treatment should focus.

The calculation of lacunarity as complement in the determination of mass-radius fractal dimension, in addition to confirm the process of increasing the homogeneity (Decreasing lacunarity values), is capable to detect heterogeneity oscillations with local increase of lacunarity. This behaviour indicates alternation between voids (lacunarity high values) and concentrations of matter (lacunarity low values) and was reported earlier by \citet{1997Chaos...7...82P}, \citet{1997PhLA..228..351S}, \citet{1998NYASA.867..258M}. Furthermore, the high values of lacunarity at radial distances smaller than 60 Mpc/h confirm the existence of voids near each centre, followed by oscillatory lacunarity behaviour (gaps and high concentration of dark matter alternatively).

The multi-fractal spectrum calculation exhibits less dispersion than that found on the average mass-radius fractal dimension and its corresponding lacunarity, as seen in the behaviour of the $\chi^2$ test. This gives more confidence on the subsequent analyses. The results for structure parameter values $q < 1$, show an excessive increase of the fractal dimension for radial distances smaller than 40 Mpc/h, reaching values that overcome the physical dimension of space (unphysical values), followed by a tendency to homogeneity for voids at larger radial distances. In the work of  \citet{2009MNRAS.399L.128S}, on transition to homogeneity in the Millennium Simulation data,  it is shown that for $q = - 2$ fractal dimension values larger than $3$ appear not only in the region below 40 Mpc/h, but also a radial distances between 70 Mpc/h and 90 Mpc/h.

For $q \geqslant 1$ we observe a slower growth of the multi-fractal dimension as the parameter structure increases, indicating a heterogeneity growth as hierarchical clustering increases. We find the possible homogeneity transition located in radial distances between 100 Mpc/h and 120 Mpc/h . This result has be corroborated by evaluating the multi-fractal dimension as a function of structure parameter $q$, where a homogeneity transition around 120 Mpc/h is detected.  This homogeneity transition scale is larger than those reported in similar works such as  \citet{2007ApJ...658...11G}, \citet{2009MNRAS.399L.128S}, but lower than the reported by \citet{2010MNRAS.405.2009Y}. Other works on galaxy surveys did not report any homogeneity transition (e.g. \citet{2008EPJB...64..615S} and  \citet{2011ARep...55..324V}). For this reason we propose to use the spectrum of lacunarity as another homogeneity transition indicator.

The lacunarity spectrum presents a similar behaviour to that found in the lacunarity calculations based in the average mass-radius dimension.  An initial increase in the heterogeneity, i.e., great presence of voids in regions close to the centres (below 20 Mpc/h) followed by a decreasing function which module oscillations with lacunarity values alternation. This shows that the clustering of dark matter haloes in the Millennium Simulation have low dense regions followed by regions of high concentration of dark matter haloes. Furthermore, the lacunarity spectrum points to a homogeneity transition (minimum lacunarity values) for radial distances between 100 Mpc/h and 130 Mpc/h,  confirming the region where the cross to homogeneity was found in the multi-fractal spectrum.

\section{Conclusions}

In this paper we have studied the distribution of  the Millennium Simulation dark matter haloes from a multi-fractal viewpoint. Our results provide a complementary method on the detection of the homogeneity transition scale based on the lacunarity concept. Our main conclusions are as follows.

The dark matter haloes clustering in the $\Lambda$CDM Millennium Simulation have a radial distance dependency divided in two regions. First a region with a fast dimensional growth which starts with great variability in regions close to the centres,  around 20 Mpc/h, a symptom of inhomogeneity at short distances. A second region where the clustering of dark matter shows evidence of slower dimensional growth arriving to homogeneity with a progressive decrease in the dispersion. Similar results are reported by \citet{0295-5075-71-2-332, 2007ApJ...658...11G} in his simulation fractal analysis. Therefore, the use of multi-fractal analysis appears indispensable because of the complexity in the clustering of the dark matter haloes.

In the determination of the fractal mass-radius dimension and multi-fractal spectrum it is necessary to consider the spurious homogenisation effects introduced as the radial distance approaches the edges of the volume to analyse. In our case when it comes near the edges, the number of possible centres are being limited (fewer centres from which to estimate the relationship); Besides, the mass-radius fractal dimension starts to exceed the spatial dimension that the clustering of dark matter haloes is embedding, this can lead to confusing behaviour over the cross to homogeneity. \citep{1989A&A...219....1W, 2002SPIE.4847...86M, citeulike:2502783}. For this reason the maximum radial distance for the fractal analysis should be limited in the calculation, according to the sample spatial size.

For the case of low-density environments, $q < 1$, it is remarkable the coherence between high dimension values for voids below $r \approx 20$ Mpc/h, the first region of fast growing in the multi-fractal dimension spectrum for  $q \geqslant 1$, and the lacunarity spectrum for the same spatial region. This is an evidence of complexity due to self-gravitation at short distances.

In the calculation of the multi-fractal spectrum, we show the homogeneity cross at depths between 100 Mpc/h and 120 Mpc/h, without exceeding the physical dimension of space. The cross to homogeneity is more precisely located by the relation between the multi-fractal dimension and the structure parameter $q$, which shows the cross to homogeneity at 120 Mpc/h radial distance. The spectrum of lacunarity confirms this result, every lacunarity function points the homogeneity scale in the same spatial region.

The lacunarity spectrum for every structure parameter in overdenses environments $q \geqslant 1$ reveals regions with relative maxima, allowing us to detect spatial regions where voids are forming inside the clustering of dark matter haloes. The transition of homogeneity for every structure parameter is founded beyond the 100 Mpc/h. Insofar that the structure parameter increases, the minimum value of lacunarity is located farther from the centre, in the same manner that the scale of homogeneity transition is increasing in the multi-fractal spectrum.

\section*{Acknowledgements}
 C. A. Chac\'on acknowledges doctoral financial support from the Universidad Distrital Francisco Jos\'e de Caldas. The authors thank the anonymous referee for the careful reading and the valuable comments and suggestions to improve this paper.


\section{Introduction}

Writing 'good' formal specifications is an essential factor for efficient program verification. Moreover, these specifications, written in the form of contracts of the software components, provide useful documentation with precisely defined requirements. They should be elegant, readable, independent from the code implementation, and should fully express the required code behaviour.

However, specifying concurrent software is not straightforward. Currently, different techniques are being developed to reason about concurrent programs and specific concurrency features~\cite{OwickiG75,SchellhornTER11,Haack76303}, but specifying the functional behaviour of such a program is still a challenge. To address this challenge, this paper introduces a history-based approach to develop behavioural specifications for concurrent data structures. We present the idea on a collection of realistic Java classes, using a combination of JML (the \emph{Java Modeling Language})~\cite{LeavensPCCRCMKC07} and \emph{separation logic}~\cite{OHearnRY01,Reynolds02separationlogic} as a specification language. This results in clear, readable and abstract specifications.

The major difficulty in specifying the behaviour of concurrent code is caused by thread interference. If threads operate on the same object, a method's postcondition can be \emph{unstable}, i.e., it can be invalidated by other threads. The following code fragment is part of the \texttt{put(E e)} method from a shared \texttt{Queue} data structure. The method enqueues an element to the end of the list that internally represents the queue:
\begin{lstlisting}[language=java, basicstyle=\footnotesize]
public void put(E e){
    lock.lock();
    enqueue(e);
    lock.unlock();   . . .
}
\end{lstlisting}
A straightforward way to define the postcondition of this method would be to state that \texttt{e} is the last element in the queue. However, after the lock release another thread may interfere, remove \texttt{e}, and invalidate the postcondition before the caller continues. Therefore, this postcondition would be unstable.

Instead, we propose to maintain a \emph{history} that records all elements that have been added to the queue. In this way, we can reason about the past. For example, for the \texttt{put} method we can use a history to state that: \emph{There was a moment in the past when \texttt{e} was added to the queue}.

The main contributions of this work are formal history-based specifications for the classes from the \texttt{java.util.concurrent.BlockingQueue} hierarchy. To make the specifications abstract, modular and naturally organized, we use abstract predicates, specification inheritance and JML model variables. This way we provide an efficient template for easily specifying the classes that implement the \texttt{BlockingQueue} interface, even if they maintain a different order on the elements. Furthermore, we claim that the history pattern is also useful for specifying other concurrent data structures; in particular we have also specified the \texttt{ConcurrentHashMap} using the same idea. Support to reason about history-based specifications will be integrated into the tool set that is currently being developed as part of the VerCors project~\cite{AmighiBHZ12}  (Verification of Concurrent Data Structures).

\section{Functional Behaviour Specifications}

In this section we explain how we use histories to overcome the problem of \emph{unstable postconditions}.
Mainly, we discuss the specifications for the \texttt{BlockingQueue} interface, and then we shortly describe how these specifications are inherited by the concrete \texttt{BlockingQueue} instantiation classes.

\paragraph{History-based Specifications}

A history in our approach is a linked list of nodes that keep the references of all queue elements. The history is defined as a specification ghost variable:
\texttt{/*@ public ghost JMLValueSequence history = new JMLValueSequence();\@*/}. Therefore, each queue element \texttt{e} has an appropriate history node \texttt{e'} that holds the reference \texttt{e}. We say that \texttt{e'} is \emph{a matching node} of \texttt{e}. History updates are done by adding JML \emph{set statements} at the locations where the queue is changed (i.e., an element is added or removed).
Since the history should remember all \emph{'old'} elements, when a queue element is deleted, its matching history node is not removed, but marked with a negative flag. We define this flag in the node class as a \texttt{boolean exists} variable.
The nodes in the history are represented by a model class \texttt{ListElem<E>}, where \texttt{E} is the type of the \texttt{BlockingQueue} items:
\begin{lstlisting}[language=java , basicstyle=\footnotesize]
/*@ model class ListElem<E> implements JMLType {
  @   boolean exist = true;
  @   E item;
  @   int orderNum = 0;
  @   public int compareOrderTo(ListElem<E> o){
  @       return orderNum - o.orderNum;
  @   }    . . .
  @ } @*/
\end{lstlisting}

\paragraph{The Order of the Elements} Our approach requires preserving the same order of the elements in the history and the queue. Therefore, it is mandatory that each node is added to the correct position in the history list. Nevertheless, different \texttt{BlockingQueue} implementations follow a different order. To keep our approach abstract, we add a \texttt{compareOrderTo(ListElem<E>)} method in the \texttt{ListElem} class, making the history nodes comparable to each other with respect to their own specific order.
We use this method to add a new node to the history at the proper position.

The common order of the elements in a \texttt{BlockingQueue} is FIFO. Thus, we define a default implementation of the \texttt{compareOrderTo} method with respect to FIFO. Since this order (as well as LIFO) depends on the moment when an element is added to the history, we add the \texttt{orderNum} variable in the \texttt{ListElem} class to represent a counter for the element insertion to the history. Hence, the \texttt{ListElem} is a default history node that respects the FIFO order. The \texttt{BlockingQueue} implementations that follow different ordering of the elements should extend this class and override the \texttt{compareOrderTo} method. An example is the \texttt{PriorityBlockingQueue} where the order of the element depends on the item value.

\paragraph{Methods' Contracts}
We define the methods' postconditions to express the behaviour of the queue in terms of its history.
The postcondition of the \texttt{add(E e)} method states that the history contains the matching node of \texttt{e}. We cannot express anything about the value of its \texttt{exist} flag because it is possible that another thread has removed the element before the end of the method:
\begin{lstlisting}[language=java , basicstyle=\footnotesize]
/*@ public behavior
  @ requires e != null;
  @ ensures (\exists ListElem n; history.has(n) * n.item == e); @*/
void put(E e) throws InterruptedException;
\end{lstlisting}
The postcondition of the \texttt{take()} method states that the matching node of the \texttt{\char92result} (the object returned as a result of this method) exists in the history with a \texttt{false} value of the \texttt{exist} flag. It is not sure that the same node existed before method start, but if it existed, its flag was \texttt{true}.
\begin{lstlisting}[language=java , basicstyle=\footnotesize]
/*@ public behavior
   @ requires true;
   @ ensures \exists ListElem n; history.has(n) * n.item == \result * !n.exist *
   @ (\old(history).contains(n) ==> \old(n).exist); @*/
E take() throws InterruptedException;
\end{lstlisting}

\paragraph{Class Invariant}

The crucial part of the specifications is to define a property to state that the queue is \emph{compatible} with the history. This ensures that the queue behaves correctly. Informally, we define a \emph{compatibility} property stating that: if all nodes with a \texttt{false} flag are removed from the history, the remaining nodes \emph{match} with the queue elements and are ordered in the same manner (see Fig.~\ref{fig: comp}).

\begin{wrapfigure}{r}{0.5\textwidth}
\vspace*{-0.5cm}
\centering
\includegraphics[scale=0.5]{full-compatible.pdf}
\caption{Compatibility between the history and the queue}
\label {fig: comp}
\end{wrapfigure}

To express this property, we define a recursive predicate \texttt{compatible}, where both parameters (\texttt{queue} and \texttt{history}) are of \texttt{JMLValueSequence} type. We could also have chosen a different type for the \texttt{queue} variable (for example, an array), but we found that this was a simple and elegant way to express the predicate.

With the approach presented so far, the order of the history nodes \texttt{e1} and \texttt{e2} is not precisely defined if \texttt{e1.compareOrderTo(e2)=0} (in this case we call \texttt{e1} and \texttt{e2} \emph{equally ordered nodes}). Note that equally ordered nodes can exist for example in the \texttt{PriorityBlockingQueue} history. If the nodes \texttt{A} and \texttt{B} in the history (Fig.~\ref{fig: comp}) were \emph{equally ordered}, we cannot guarantee that \texttt{A} precedes or follows \texttt{B}. Anyhow, the queue and the history should be compatible and therefore, they should follow the same order. The following formal specification of the \texttt{compatible} predicate deals with this case:

\begin{lstlisting}[language=java, basicstyle=\footnotesize]
/*@ pred compatible(JMLValueSequence queue, JMLValueSequence history) =
  @ (queue.isEmpty() && history.isEmpty()) ||                              //case 1
  @ (!queue.isEmpty() &&
  @  (\exists int idx; 0<=idx<history.length;
  @    (((ListElem)queue.first()).item == ((ListElem)history.itemAt(idx)).item) &&
  @    ((ListElem)history.itemAt(idx)).exist &&
  @    (\forall int i; 0<=i<idx;
  @     !((ListElem)history.itemAt(i)).exist ||
  @     (((ListElem)history.itemAt(i)).compareOrderTo((ListElem)history.itemAt(idx))
  @	    == 0))
  @    && compatible(queue.removeItemAt(0), history.removeItemAt(idx))))   //case 2
  @*/
\end{lstlisting}

With this definition we even allow expressing the compatibility between queue and a history where the order of elements is not important at all (a set data structure for example). In that case, the method \texttt{compareOrderTo} should be overridden to return always 0, and the \texttt{compatible} predicate will express that the queue contains the correct elements without considering their order.

Furthermore, a class invariant is specified to state that the queue is compatible with the history. To state the invariant, we use a model variable
\texttt{queue} to denote the current contents of the queue.

\begin{lstlisting}[language=java , basicstyle=\footnotesize]
/*@ protected model instance JMLValueSequence queue;
  @ invariant compatible(queue,  history) @*/
\end{lstlisting}

\paragraph{Specifications for the Descendant Classes}

Specifying the descendant classes is simple: the following steps are required:
 \begin{inparaenum}[(\itshape i\upshape)]
\item if the order of the elements in the queue is different then FIFO, extend the \texttt{ListElem} class and override the \texttt{compareOrderTo(ListElem<E>)} method;
\item update the \texttt{history} ghost variable at the locations where the queue is changed i.e., an element is added or removed).
\end{inparaenum}

The complete specifications together with the Java source code are available via \url{http://www.ewi.utwente.nl/~zaharievam/specs_classes/}.

\section{Conclusions}
This paper presents a history-based approach for specifying the functional behaviour of concurrent data structures. It allows one to write \emph{stable} method specifications, i.e., specifications
that cannot be invalidated by other threads.  We apply this approach on the \texttt{java.util.concurrent.BlockingQueue} hierarchy, but it is also suitable for specifying other data structures. Furthermore, our approach introduces a new dimension: abstraction. The result is a specification framework, where specifying a new \texttt{BlockingQueue} implementation requires only a small effort.

\paragraph{Related Work}

Using a history for specifying communication-based parallel programs dates back to the 70's and 80's. A history, also called a \emph{communication trace}, is a sequence that records the communication events sent between the processes. The history helps one to reason about the events that happened up to a certain point in time and their order. This approach has been used to specify the interactions in CSP~\cite{Soundararajan83}. Also in recent work, a similar idea has been used: Dovland \emph{et al.}~\cite{Dovland:2008} use histories to record the communication events between the components in a distributed environment.

The specifications we provide are closely related to the linearizability-based approach where proving that a data structure is linearizable guarantees its correctness. In this approach the specification contains sequential code that updates a ghost variable, which is an abstract representation of the data structure itself. A library is linearizable with respect to a specification if it behaves as if it was sequentially executed. Linearizability is automatically proved by searching a \emph{linearization point} for each method ~\cite{Vafeiadis10}.
 Gotsman \emph{et al.}~\cite{linown11} explain the importance of linearizability stating that: when proving properties for a client that uses a concurrent library, this library can be soundly replaced by its abstract implementation.

Our work is more practically oriented. We are focused on making clear and understandable specifications, written in a standard form of method contracts and class invariants. The specifications are readable even without the code, and therefore can be used as formally written documentation.

\paragraph{Future Work}

The work described here are only the first steps in the VerCors
project~\cite{AmighiBHZ12}, and there is still much work to be done in
the future. A first important step is to develop tool support, so that
the history-based specifications can be validated w.r.t. the API
reference implementation.

Another point we are currently working on is providing full support for class invariants in a concurrent environment.
These class invariants must be \emph{strong}, i.e., invariants that hold in all possible states.
We allow an invariant to explicitly break under certain conditions, while guaranteeing that no other thread is able to observe that it is broken.

It still has to be investigated how the history-based specifications
can be exploited to verify client's code. One particular case where
this will be useful is if the behaviour of the client depends on the
elements that have been in the queue at some point. For example, in a
game application, the queue could hold the \(k\) highest scores in
descending order. A client can use the history-based specifications to
derive properties as: the queue holds the best scores, and it knows
the correct number of total scores.

Finally, we will also investigate whether we can automatically
generate updates to the history, instead of depending on the user
manually inserting these updates in the descendant classes.

\bibliographystyle{eptcs}