\section{Introduction}

A Bell inequality \cite{Brunner2014} experiment consists of $n$ parties who share an entangled state of $n$ (or more) particles.  Each party chooses to perform one of $m$ measurements, with each measurement producing one of $d$ outcomes.  By repeating this round multiple times, the probabilities for the various joint outcomes can be estimated.  A Bell inequality for this scenario is a relationship that these probabilities must satisfy, if they arise from a local realistic model.  The general structure of these inequalities has been studied intensively since Bell's original argument \cite{Bell65}.

The most famous form of the Bell inequality is the Clauser-Horne-Shimony-Holt (CHSH) inequality \cite{Clauser69} for $n=2$ parties (e.g. Alice and Bob), each party performing one of $m=2$ measurement choices, with each measurement registering one of $d=2$ outcomes.  In terms of the joint probabilities, the CHSH inequality reads
\begin{equation}
p(A_1 = B_1) - p(A_1=B_2) + p(A_2=B_1) + p(A_2=B_2) \le 2,
\label{CHSH}
\end{equation}
where the measurement settings are labeled by $1$ and $2$ for Alice's (or Bob's) choice of measurement, with outcomes $A_1$ and $A_2$ for Alice (and $B_1$ and $B_2$ for Bob), and, in a slight abuse of notation, the joint probability $p(A_0=B_0)$ indicates the probability that Alice and Bob's measurement outcomes are identical.  For certain entangled states and measurement choices, quantum mechanics predicts, and experiments confirm, a violation of this inequality \cite{Brunner2014}.

A common experimental procedure, illustrated in Fig. \ref{bellfig}, replaces the alternative measurement settings by a unitary transformation chosen by the parties and performed just before a fixed measurement.  This unitary must be chosen and performed by one party sufficiently quickly to ensure that no information can propagate to the other party.   If this cannot be done, the experiment is subject to the so-called {\it locality} loophole, for which the derivation of Eq. (\ref{CHSH}) fails.  Another requirement is that the measurements be sufficiently accurate so that the probabilities entering the inequality can be reliably estimated.  If this is not the case, the experiment is subject to the so-called {\it detection} loophole.  Modern experiments with photons typically close the locality loophole \cite{Weihs98,Ursin2007}, but are subject to the detection loophole (in the guise of the fair-sampling assumption, but see recent progress \cite{Christensen2013,Giustina2013}).  Experiments with matter qubits (atoms \cite{Hofmann2012}, ions \cite{Rowe2001,Matsuk2008}, solid-state spin qubits \cite{Pfaff2013,Dehollain2015} and superconducting qubits \cite{Ansmann2009, Vlastakis2015}) typically close the detection loophole, but are subject to the locality loophole.   Achieving a loophole-free Bell inequality experiment would be a landmark test of quantum mechanics.

\begin{figure}
\includegraphics[width=3 in]{bell_inequality_setup}
\caption{A general framework for a Bell inequality experiment with $n=2$ parties.  A source of entanglement is shared to Alice and Bob, who choose to rotate their part of the joint quantum state by unitaries $U_{A,a}$ and $U_{B,b}$, respectively.  The measurement outcomes $M_A$ and $M_B$ are then recorded, compiled, and compared against the inequality.  }
\label{bellfig}
\end{figure}

Towards this goal, and to better understand entanglement and nonlocality in general, there have been many studies of generalized Bell inequalites  \cite{Brunner2014}.  These include the Mermin inequality \cite{Mermin90} for multiple qubits ($n>2$), the Collins-Gisin inequality \cite{Collins2004} for multiple measurements ($m>2$), and the Collins-Gisin-Linden-Massar-Popescu (CGLMP) inequality \cite{Collins2002} for higher-dimensional systems known as qudits ($d>2$).  A recent approach \cite{Vertesi2010} to the Collins-Gisin inequality for entangled qudits (with  $n=2$ and $m=d>2$) exhibits the potential to reduce the requirements to close the detection loophole.

In this paper we focus on understanding the CGLMP inequality (for $n=m=2$ and $d>2$), which takes the form
\begin{equation}
I_d = \sum_{k=0}^{\lfloor d/2\rfloor -1} \left(1 - \frac{2k}{d-1} \right) \left[ \mathcal{P}(k) - \mathcal{P}(-k-1) \right] \le 2,
\end{equation}
where
\begin{eqnarray}
\mathcal{P}(k) &=& P(A_1 = B_1+k) + P(B_1 = A_2 + k + 1)  \nonumber \\
& & + P(A_2 = B_2 + k) + P(B_2 = A_1 + k).
\end{eqnarray}
Here the joint probabilities are defined for outcomes $A_a = 0, 1, \dots, d-1$, and the addition is performed modulo $d$.  These can expressed as
\begin{equation}
P(A_a = B_b + k) = \sum_{j=0}^{d-1} P(A_a = j, B_b = j + k \ \mbox{mod} \ d).
\end{equation}
We have also studied a closely related inequality proposed by Zohren and Gill \cite{Zohren2008}
\begin{eqnarray}
P(A_2 < B_2) + P(B_2 < A_1) + P(A_1 < B_1) & & \nonumber \\
+ P(B_1 \le A_2) &\ge& 1.
\end{eqnarray}
These two inequalities have the remarkable property that the violation increases with increasing $d$.  Alternatively, if the initial state is subject to depolarizing noise, the amount of noise that removes the entanglement, and hence the potential to violate the inequality, increases with the dimension \cite{Kasz2000}.  In this sense, these higher-dimensional Bell inequalities exhibit a surprising robustness to noise, and may be useful when exploring advanced tests of quantum systems.

A recent experiment tested the CGMLP inequality using orbital angular momentum (OAM) states of light \cite{Dada2011}, but did not find the enhanced violation with dimension.  In fact, they found that the violation ceased for dimensions higher than $d=12$.  This fact strongly suggests that the robustness of the inequality can be compromised by the actual noise subject to the system.  Previous work on this topic has focused only on a simple form of depolarizing noise \cite{Kasz2000,Collins2002}, and its impact on the detector efficiency needed to observe a violation of the inequality \cite{Durt2001}.  To account for the behavior seen in the recent experiment, and to predict the possible violations in other experiments, more realistic noise models are necessary.

In general, noise can affect all of the stages of the experiment: entangled-state preparation, state rotation, and state measurement.  The noise from each stage must be analyzed to determine the robustness of the inequality.   In this paper we examine a general framework for testing the CGLMP inequality with qudits, analyze the potential complexity of each stage of the experiment, and evaluate the impact of different types of noise (depolarizing, dephasing, and amplitude damping) on the inequality.  Our work is inspired by theoretical proposals to use the quantum states of superconducting resonators as qudits \cite{Strauch2010, Strauch11,  Krastanov2015}; such systems have the potential for long-distance entanglement through microwave photons \cite{Yin2013,Wenner2014}.  However, our analysis is intended to be general enough to guide experiments with other potential matter qudit systems, and our results may have implications for photonic qudits as well.  Our analysis indicates that the CGLMP inequality does not generally exhibit the special robustness to noise claimed in previous work.

This paper is organized as follows.  Section II describes a general Bell inequality experiment with ideal qudits, analyzing  how state preparation, rotation, and measurement can be implemented for a general qudit system.  Section III analyzes the potential inequality violation subject to depolarizing, amplitude-damping, and dephasing noise on the qudits during the full experimental sequence.  In Section IV we conclude and discuss outstanding questions.

\section{Bell Inequality Experiment with Ideal Qudits}

In this section we will describe an analysis of a bipartite ($n=2$) Bell inequality experiment, focusing on the general structure of entangled-state preparation, single-qudit rotations, and qudit-state measurement.  A quantum circuit for this process is shown in Fig. \ref{bellcircuit}.  We will analyze each component of this circuit, drawing inspiration from recent theoretical work on the control of superconducting resonators.  However, we expect our results can be applied to alternative qudit implementations.

\begin{figure}
\includegraphics[width=3 in]{bell_circuit}
\caption{A quantum circuit for a Bell inequality experiment, consisting of entangled-state preparation, single-qudit rotations, and qudit-state measurement.  In the first stage, two qudits, initially in the $|00\rangle$ state, are acted on by generalized Hadamard gates $H$, a controlled-phase gate (with $\theta = 2\pi /d$) to prepared a maximally entangled two-qudit state.  The single-qudit rotations $U_a$ and $U_b$ determine the two measurement bases (depending on the choices of Alice and Bob), while the actual measurement is in a fixed qudit basis. }
\label{bellcircuit}
\end{figure}

\subsection{Entangled-State Preparation}
The first important stage of a Bell inequality experiment is the production of entangled states of two systems.  As shown in Fig. \ref{bellcircuit}, this can be accomplished by using a generalized Hadamard or DFT (for Discrete-Fourier-Transform) operation on each qudit and a controlled-phase gate between the qudits.

The generalized Hadamard gate is the discrete Fourier transform (DFT) matrix with matrix elements $H_{j,k} = \omega^{jk}/\sqrt{d}$, where $\omega = e^{2\pi i /d}$.  For $d=4$, this has the matrix form
\begin{equation}
H(d=4) = DFT_4 = \frac{1}{\sqrt{4}} \left( \begin{array}{rrrr} 1 & 1 & 1 & 1 \\
1 & i & -1 & -i \\
1 & -1 & 1 & -1 \\
1 & -i & -1 & i
\end{array} \right).
\end{equation}
In general this matrix has $d$ independent elements (and $d^2$ total elements), so an analysis of control on qudits suggests that implementing a general qudit Hadamard will require a number of elementary operations that is a polynomial in $d$ \cite{Brennen05}.

The controlled-phase gate is a two-qudit gate that implements the transformation
\begin{equation}
\mathcal{CR}(\theta)|j,k\rangle = e^{-i (j k) \theta} |j,k\rangle.
\end{equation}
This is a natural generalization of the two-qubit controlled phase gate, and can also be implemented in a time polynomial in $d$ \cite{Brennen05}.

The combination of these operations in the preparation stage of Fig. \ref{bellcircuit} performs the transformation
\begin{eqnarray}
|\Psi\rangle &=& (H \otimes I) \mathcal{CR}(2\pi/d) (H \otimes H) |0,0\rangle \nonumber \\
&=& (I \otimes H) \mathcal{CR}(2\pi/d) \frac{1}{d} \sum_{j,k=0}^{d-1} |j,k\rangle \nonumber \\
&=& (I \otimes H) \frac{1}{d} \sum_{j,k=0}^{d-1} e^{-i 2 \pi (j k)/d} |j,k\rangle \nonumber \\
&=& \frac{1}{d^{3/2}} \sum_{j,k,\ell=0}^{d-1} e^{-i 2 \pi k (j-\ell)/d} |j,\ell\rangle \nonumber \\
&=& \frac{1}{\sqrt{d}} \sum_{j} |j,j\rangle.
\label{maxentstate}
\end{eqnarray}
This maximally entangled two-qudit state will be used in the Bell inequality measurements for Alice and Bob.

We note that there are alternative approaches to generating the initially entangled state for superconducting resonators, which typically scale linearly in the qudit dimension \cite{Sharma2015}.  Qudit operations can also be optimized to scale linearly with qudit dimension \cite{OLeary06}.  We optimistically conclude that the state-preparation stage can be performed in a timescale that is linear in $d$.

\subsection{Single-Qudit Rotations}

After Alice and Bob have chosen their measurement basis ($a$ or $b$), they will adjust their measurement by rotating their half of the joint-qudit state by one of the unitary operators $U_{A,a}$ or $U_{B,b}$.  These have the matrix elements
\begin{equation}
[U_{A,a}]_{j,k} = \frac{1}{\sqrt{d}} e^{ i 2 \pi (jk)/d} e^{i \alpha_a j}
\end{equation}
and
\begin{equation}
[U_{B,b}]_{j,k} = \frac{1}{\sqrt{d}} e^{-i 2 \pi (jk)/d} e^{i \beta_b j},
\end{equation}
where $\alpha_1 = 0$, $\alpha_2 = 1/2$, $\beta_1 = 1/4$, and $\beta_2 = -1/4$.  These operations involve the discrete Fourier transform (or its inverse), along with diagonal phases.  These again can be performed in a time polynomial in $d$, which we again optimistically take as linear in the qudit dimension \cite{OLeary06}.

\subsection{State Measurement}

In order to close both the locality and detection loopholes, one wants a fast and efficient measurement of the qudit states.  For some qudit implementations, such as hyperfine states of atoms or multilevel superconducting devices, one can implement a direct $d$ outcome measurement.   For qubit-resonator systems in trapped ion, cavity, or circuit-QED systems, one often implements an indirect measurement coupling the resonator to one or more qubits.   We consider a state measurement approach that maps the resonator state of $d=2^n$ dimensions onto $n$ qubits, inspired by the discussion in Chapter 6 of \cite{HarocheBook}.  After this mapping, one can measure the state bit by bit.   This allows a single-shot measurement of the qudit state.

\begin{figure*}
\includegraphics[width=6 in]{measurement_circuit}
\caption{A measurement circuit to map the quantum state of a resonator (of dimension $d = 2^n$) onto the state of $n$ qubits.  This circuit is composed of $2n$ single-qubit Hadamard operations,  $n$ qubit-resonator controlled-phase gates (with phases indicated in the circular gates on the resonator), and $n(n+1)/2$ qubit-qubit controlled-phase gates (with phases indicated in the rectangular gates on the qubits).}
\label{measurementfig}
\end{figure*}

The measurement scheme can be represented by a quantum circuit, illustrated in Fig. \ref{measurementfig}, composed of single-qubit Hadamard gates, two-qubit controlled-phase gates, and a special qubit-resonator controlled-phase gate.  This last gate acts on a qubit-resonator state as follows:
\begin{equation}
|x\rangle_{\subs{qubit}} \otimes |y\rangle_{\subs{res}}  \to e^{-i \theta x y} |x\rangle_{\subs{qubit}} \otimes |y\rangle_{\subs{res}}.
\end{equation}
This gate can be implemented using the natural evolution of a qubit-resonator system in the dispersive regime \cite{Strauch11}, resonant qubit-resonator logic operations \cite{Strauch2012}, or other methods \cite{Mischuck2013,Krastanov2015}.  In brief, the measurement circuit performs the mapping by sequentially correlating the resonator state $|y\rangle$ with a register of qubits ($\{x_1, x_2, \dots, x_n\}$).   Using a bit string $\{y_1, \dots, y_n\}$ to label the resonator state $|y\rangle$, with the (big-endian) binary representation
\begin{equation}
y = y_1 2^{n-1} + y_2 2^{n-2} + \cdots + y_{n-1} 2 + y_n,
\end{equation}
the circuit performs the mapping
\begin{equation}
|0 0 \cdots 0 \rangle_{\subs{qubits}} \otimes |y\rangle_{\subs{res}} \to |y_1 y_2 \cdots y_n \rangle_{\subs{qubits}} \otimes |y\rangle_{\subs{res}}.
\end{equation}
Then, measurement of the qubits will provide a bit-by-bit measurement of the resonator state $|y\rangle$.

We now analyze the measurement circuit of Fig. \ref{measurementfig} in more detail.  The qubits are initially prepared in an equal superposition state (by Hadamard gates):
\begin{equation}
|\Psi_0\rangle = \frac{1}{2^{n/2}} \sum_{\{x_1, \dots, x_n\}, y} c_y |x_1 \cdots x_n \rangle \otimes |y \rangle
\end{equation}
The first qubit-resonator controlled-phase gate produces the state
\begin{equation}
\frac{1}{2^{n/2}} \sum_{\{x_1, \dots, x_n\}, y} c_y (-1)^{x_n y} |x_1 \cdots x_n \rangle \otimes |y\rangle.
\end{equation}
Using the binary representation, we have
\begin{eqnarray}
(-1)^{x_n y} &=& (-1)^{x_n y_1 2^{n-1}} (-1)^{x_n y_2 2^{n-2}} \cdots(-1)^{x_n y_n} \nonumber \\
&=& (-1)^{x_n y_n}.
\label{phaseeqn}
\end{eqnarray}
A subsequent Hadamard on qubit $n$ yields
\begin{equation}
\frac{1}{2^{(n+1)/2}} \sum_{\{x_1, \dots, x_n,z\},y} c_y (-1)^{x_n y_n + x_n z} |x_1 \cdots x_{n-1} z\rangle \otimes |y\rangle.
\end{equation}
Performing the sums over $x_n$ and $z$, using the fact that
\begin{equation}
\sum_{x_n=0}^1 (-1)^{x_n y_n + x_n z} = 2 \delta_{y_n,z},
\label{hadamardsum}
\end{equation}
we find
\begin{equation}
|\Psi_1\rangle = \frac{1}{2^{(n-1)/2}} \sum_{\{x_1, \dots, x_{n-1}\}, y} c_y |x_1 \cdots x_{n-1} y_n \rangle \otimes |y\rangle.
\end{equation}

The next qubit-resonator phase gate produces the controlled phase
\begin{equation}
\left(e^{i \pi/2}\right)^{x_{n-1}y} = \left(-1\right)^{x_{n-1} y_{n-1}} (i)^{x_{n-1} y_{n}},
\end{equation}
where we have used an argument similar to Eq. (\ref{phaseeqn}).  The final phase can be eliminated by a qubit-qubit phase gate, while a subsequent Hadamard on qubit $n-1$ yields
\begin{equation}
\frac{1}{2^{n/2}} \sum_{\{x_1, \dots, x_{n-1},z\},y} c_y (-1)^{x_{n-1} y_{n-1} + x_{n-1} z} |x_1 \cdots x_{n-2} z y_n \rangle \otimes |y\rangle.
\end{equation}
Performing the sums over $x_{n-1}$ and $z$, using the corresponding form of Eq. (\ref{hadamardsum}), the second stage yields
\begin{equation}
|\Psi_2\rangle = \frac{1}{2^{(n-2)/2}} \sum_{\{x_1, \dots, x_{n-2}\}, y} c_y |x_1 \cdots x_{n-2} y_{n-1} y_n \rangle \otimes |y\rangle.
\end{equation}

After $k$ stages we will find
\begin{eqnarray}
|\Psi_k\rangle &=& \frac{1}{2^{(n-k)/2}} \times \nonumber \\
& & \sum_{\{x_1, \dots, x_{n-k}\},y} c_y |x_1 \cdots x_{n-k} y_{n-k+1} \cdots y_n \rangle \otimes |y\rangle \nonumber \\
\end{eqnarray}
The next qubit-resonator gate produces the controlled phase
\begin{eqnarray}
\left(e^{i \pi/2^k}\right)^{x_{n-k}y} &=& \left(-1\right)^{x_{n-k} y_{n-k}} \times \nonumber \\
& &  \left(e^{i \pi/2}\right)^{x_{n-k} y_{n-k+1}} \cdots \left( e^{i \pi/2^k}\right)^{x_{n-k} y_n}. \nonumber \\
\end{eqnarray}
Eliminating all of the trailing controlled phases by qubit-qubit phase gates, performing a Hadamard gate on qubit $n-k$, and evaluating the summation over $x_{n-k}$ yields
\begin{eqnarray}
|\Psi_{k+1}\rangle &=& \frac{1}{2^{(n-k-1)/2}} \times \nonumber \\
& & \sum_{\{x_1, \dots, x_{n-k-1}\},y} c_y |x_1 \cdots x_{n-k-1} y_{n-k} \cdots y_n \rangle \otimes |y\rangle. \nonumber \\
\end{eqnarray}

Repeating until $k=n$, we finally obtain
\begin{equation}
|\Psi_n\rangle = \sum_y c_y |y_1 \cdots y_n \rangle \otimes |y\rangle.
\end{equation}
A measurement of the $n$ qubits will produce the outcome $(y_1, y_2, \dots, y_n)$ with probability $|c_y|^2$.  This circuit uses $n$ qubit-resonator gates, $2n$ Hadamard gates, and $n(n-1)/2$ qubit-qubit controlled-phase gates.  Thus, we see that this circuit requires a timescale of order $(\log_2 d)^2$.

\section{Bell Inequality Experiment with Noisy Qudits}

The preceding analysis shows that, in general, one can expect that the state preparation, rotation, and measurement stages of a Bell inequality experiment will each scale with the qudit dimension.   If these stages are composed of fundamental elements, such as quantum gates, that are each subject to noise, then the resulting experiment will be subject to noise that scales with the qudit dimension.  How that scaling affects the inequality is the subject of this section.

While we will analyze idealized models of noisy qudits, our results can be understood in physical terms.  Each fundamental quantum gate will take some time, during which the quantum state can be subject to fluctuating fields, lose energy to the environment, or lose quantum information in some other fashion.  Each of these noise processes will affect the quantum state of the system, in terms of its density matrix.  If the total time for the experiment scales with the qudit dimension, the resulting density matrix will be subject to a correspondingly increased amount of noise.

A detailed analysis of this process would require modeling the qudits' Hamiltonian and its coupling to external fields and the environment.  To understand the scaling with dimension, however, we can simplify our analysis to consider a single parameter figure-of-merit for each step of the experiment, such as the gate fidelity, and look at three types of noise: depolarizing noise, amplitude damping, and dephasing noise \cite{mikeandike}.  Using these models, we allow the strength of the noise to vary with the qudit dimension in the following way.

We take the initial density matrix $\rho_0 = |\Psi\rangle \langle \Psi|$, with $|\Psi\rangle$ given by Eq. (\ref{maxentstate}), and iterate the appropriate trace-preserving map
\begin{equation}
\rho_{n+1} = \sum_{m} E_m \rho_n E_m^{\dagger},
\end{equation}
where the error operators $E_m$ depend on the type of noise, to be defined below, and are parametrized by a single number $p$.   This map is iterated $N$ times, where $N$ is proportional to the number of fundamental steps during the preparation and rotation stages of the Bell circuit of Fig. \ref{bellcircuit} and $p$ is proportional to the fidelity of each step.  The final density matrix, after $N$ applications of the noise map and subject to the measurement choices $a$ and $b$, is given by
\begin{equation}
\rho_{a,b} = (U_{A,a} \otimes U_{B,b}) \rho_N (U_{A,a} \otimes U_{B,b})^{\dagger}.
\end{equation}
The probabilities of the measurement outcomes are then
\begin{equation}
P(A_a=j, B_b=k) = \langle j,k| \rho_{a,b}|j,k\rangle.
\label{jointprob}
\end{equation}

As discussed above, the number of steps required to prepare the initial state and rotate the measurement is taken to be linear in $d$, while the single-shot measurement circuit scales with $(\log d)^2$.  For the dimensions $d \le 16$ considered below, there is little difference between linear and logarithmic scaling.  Thus, to understand how scaling affects the Bell parameter $I_d$, we set $N=d$.  For comparison, we will also consider $N=1$, similar to previous work.

\begin{figure*}[t]
\includegraphics[width=6 in]{depolarfig}
\caption{Generalized Bell parameter $I_d$ as a function of qudit dimension $d$ for two types of depolarizing noise and various noise strengths.  (a) Single application of the depolarizing operator, with probabilities of $1$, $0.99$, $0.95$, and $0.9$, from top to bottom.   (b) Repeated application of the depolarizing operator with the number of applications is linear with the qudit dimension, with probabilities of $1$, $0.99$, $0.95$, and $0.9$, from top to bottom.   }
\label{depolarfig}
\end{figure*}

\subsection{Depolarizing Noise}

Depolarizing noise simulates the interaction of a system with a high temperature environment.  With each iteration, the system has probability $1-p$ of becoming depolarized, i.e. replaced by the completely mixed state, so that
\begin{equation}
\rho_{n+1} = p \rho_n +  (1-p) \frac{I}{d^2},
\label{depolar}
\end{equation}
while the system has probability $p$ of being unaffected.   Note that for this map, smaller values $p$ indicate stronger noise (this convention is chosen to match \cite{Collins2002}).  As described above, we allow this map to be repeated for a variable number of times before calculating the Bell parameter.

The value of the generalized Bell parameter $I_d(p)$ is shown in Fig. \ref{depolarfig}, as a function of $d$ and for several values of the probability $p$.   From these results, we can make several observations.  First, for $p=1$, we see that the Bell parameter increases with dimension, as found previously.  In fact, we find
\begin{equation}
I_d(p=1) \approx 2.97 \left( 1 - \frac{1}{10 d} \right).
\label{bellincrease}
\end{equation}
Second, for $p < 1$ and $N=1$, the Bell parameter again increases with dimension, as seen in Fig. \ref{depolarfig}(a).  However, the overall violation does decrease with decreasing $p$.  In fact, it can be shown that, for this type of depolarizing noise, $I_d(p) = p I_d(1)$ \cite{Collins2002}.  This behavior underlies the reputed robustness to noise.

By constrast, if we apply the map $N=d$ times before calculating $I_d$ [Fig. \ref{depolarfig}(b)], we find that the violation decays with dimension.  This can be understood by noting that repeated iteration of Eq. (\ref{depolar}) has the solution
\begin{equation}
\rho_d = p^d \rho_0 + (1-p^d) \frac{I}{d^2},
\end{equation}
so that, in the presence of depolarizing noise scaling with $d$, the Bell parameter behaves as $p^d I_d(p=1)$.  Thus, the small increase in the Bell parameter (for $p=1$) of Eq. (\ref{bellincrease}) is quickly reduced by the exponentially decreasing factor $p^d$.

This may have some bearing on the experimental results discussed above \cite{Dada2011}.  If one of the stages of the experiment has a fidelity that decays with the dimension of the entangled state, as evidenced in previous experiments \cite{Jack2009}, one can easily reproduce the observed decay of the inequality violation.  For example, setting $p=0.998$ and $N=d$ produces results in rough agreement with Fig. 3 of \cite{Dada2011}.

\subsection{Amplitude Damping Noise}

Amplitude damping describes the effects of energy dissipation on the quantum system.  We use a simplified amplitude damping model, chosen to represent the damping of a quantum resonator over a time $\Delta t$ with decay time $T$.  In this model, the singly excited state $|1\rangle$ decays with probability $p = e^{-\Delta t/T}$, while state $|j\rangle$ decays with probability $p^j$ \cite{Wang2008}.

The specific model of amplitude damping is given by the map
\begin{equation}
\rho_{n+1} = \sum_{\ell,m=0}^1 (E_{\ell} \otimes E_m) \rho_{n} (E_{\ell} \otimes E_m)^\dagger,
\end{equation}
where the single-qudit amplitude-damping operators are given by
\begin{equation}
E_0 = \sum_{j=0}^{d-1} \sqrt{p^j} |j\rangle \langle j|
\end{equation}
and
\begin{equation}
E_1 = \sum_{j=1}^{d-1} \sqrt{1-p^j} |j-1\rangle \langle j|.
\end{equation}
We note that while this model intrinsically scales with dimension, we continue to allow the number of iterations of the map to scale with dimension as well.  It is also relevant to observe that this map is an approximation to the actual decay process, in that single application of this map with $p=0$ corresponds not to complete loss of energy, but to the subtraction of a single photon.

Using this model of amplitude damping, the resulting Bell parameter $I_d$ as a function of dimension $d$ is shown in Fig. \ref{ampdampfig}, for several values of the probability $p$.  For $N=1$, the violation decays slowly with $d$, but remains a violation for up to $d=16$ for $p=0.9$.  Already, however, we see a significant difference from the behavior of $I_d$ under depolarizing noise.

For $N=d$ iterations, the violation decays much faster with $d$.  However, the decay is surprisingly less than the depolarizing map for the same value of $p$.  We will soon return to this issue, but for certain values of $d$ and choices of the probability, the two are similar.  For example, to model the OAM experiment \cite{Dada2011}, one could use amplitude damping with $p=0.992$ and $N=d$.

\begin{figure*}
\includegraphics[width=6 in]{ampdampfig}
\caption{Generalized Bell parameter $I_d$ as a function of qudit dimension $d$ for two types of amplitude damping and various noise strengths.  (a) Single application of the amplitude damping operators, with probabilities of $1$, $0.99$, $0.95$, and $0.9$, from top to bottom.   (b) Repeated application of the amplitude damping operators, with the number of applications is linear with the qudit dimension, with probabilities of $1$, $0.99$, $0.95$, and $0.9$, from top to bottom.   }
\label{ampdampfig}
\end{figure*}

The results above have used a discrete model of energy loss, in which a single quantum is removed from the system with some probability.  Of course, the actual physics involves a continuous loss of energy.  We have compared our discrete amplitude damping model to a more continuous model, by breaking up the decay into a number of intervals with the same overall decay factor $p$.  As one might expect, the difference between the two models is negligible for $p$ near 1, in which the probability of losing two quanta can be neglected.

\subsection{Dephasing Noise}
Dephasing describes loss of quantum information without loss of energy: rather than changing the amplitudes of the states as a function of time, the energy eigenstates of a system accrue random phases with some probability.   This is modeled by
\begin{equation}
\rho_{n+1} = p \rho_{n} + (1-p) \sum_{j,k=0}^{d-1} \left(\langle j, k| \rho_n |j, k\rangle \right) |j, k\rangle \langle j,k|.
\end{equation}
After each iteration, the off-diagonal elements of $\rho_n$ are reduced by a factor $p$.

We note that the off-diagonal elements of the final density matrix in this case are identical to those found for the depolarizing map of Eq. (\ref{depolar}).  Remarkably, when measuring the density matrix using the DFT operators, the joint probabilities of Eq. (\ref{jointprob}) do not depend on the diagonal elements.  We thus arrive at the interesting result the effect of dephasing noise is identical to that of the depolarizing noise considered previously.

To verify that the probabilities do not depend on the diagonal elements, we consider an arbitrary ``diagonal'' density matrix of the form
\begin{equation}
\rho_{\subs{diag}} = \sum_{j,k=1}^d c_{j,k} |j,k\rangle \langle j,k|,
\end{equation}
where
\begin{equation}
\sum_{j,k=1}^d c_{j,k} =1.
\end{equation}
For this density matrix we use Eq. (\ref{jointprob}) to calculate the joint probablitiies
\begin{eqnarray}
P(\ell,m) &=& \langle \ell,m| U_{A,a} \otimes U_{B,b} \rho_{\subs{diag}} \left(U_{A,a} \otimes U_{B,b} \right)^{\dagger} |\ell,m\rangle \nonumber \nonumber \\
&=& \sum_{j,k=1}^{d} c_{j,k} | \langle \ell | U_{A,a} |j \rangle|^2 |\langle m| U_{B,b} | k\rangle|^2 \nonumber \\
&=& \sum_{j,k=1}^{d} c_{j,k} \frac{1}{d^2} = \frac{1}{d^2}.
\end{eqnarray}
Since these probabilities are independent of $c_{j,k}$, the diagonal elements of any density matrix will contribute the constant value of $1/d^2$ to the probability $P(\ell,m)$.  Of course, the off-diagonal elements will also contribute to the probability.  However, any two density matrices with the same off-diagonal elements will lead to the same probabilities, confirming the result claimed above.

We can also use this observation regarding the probabilities to help understand the difference between amplitude damping and depolarizing / dephasing noise.  For a single qubit, amplitude damping reduces the diagonal elements by $p$, but the off-diagonal elements by $\sqrt{p} > p$.  Thus, the effect of amplitude damping on the probabilities will be less than that for depolarizing / dephasing noise, for the same value of $p$.  That this holds true for dimensions $d>2$ is still somewhat surprising.  One possible explanation is that coherence between neighboring qudit states is largely maintained, much like the decay of a coherent oscillator state \cite{HarocheBook}, but we have not explored this conjecture further.

\subsection{Thresholds for Inequality Violation}

A convenient way to summarize the results obtained above is to consider the minimum probability $p_{\subs{min}}$ for which the CGLMP inequality can be violated.  This probability satisfies
\begin{equation}
I_d(p_{\subs{min}}) = 2,
\end{equation}
and depends on the different models of noise.  Recall that a smaller probability $p$ indicates a greater amount of noise. For $p < p_{\subs{min}}$, the inequality is not violated, so we call $p_{\subs{min}}$ the threshold probability.  Previous work found that $p_{\subs{min}}$ decreases with dimensions, with the conclusion that higher dimensional Bell inequalities can be violated even if the system is subject to greater amounts of noise.

We have calculated the minimum probabilities for the various noise models described above and their dependence on dimension, as shown in Fig. \ref{threshold}.  We find that only for depolarizing / dephasing noise with $N=1$ does the threshold increase with dimension.  For all other cases---depolarizing / dephasing noise with $N=d$, and amplitude damping with $N=1$ and $N=d$---the minimum probability increases with dimension.  For these noise models, a violation of the CGLMP inequality with higher-dimensional systems requires higher-fidelity operations.

\begin{figure}
\includegraphics[width=3 in]{CGLMP_threshold}
\caption{The minimum probability $p_{\subs{min}}$ for which the CGLMP inequality is violated, as a function of qudit dimension $d$ and for the noise models (from top to bottom) depolarizing / dephasing noise with $N=d$, amplitude damping with $N=d$, amplitude damping with $N=1$, and depolarizing / dephasing noise with $N=1$.  Systems with noise $p < p_{\subs{min}}$ will not violate the inequality. }
\label{threshold}
\end{figure}

\subsection{Alternative Procedures}

We have studied a number of alternative procedures of the Bell inequality test.  First, we have looked at the inequality proposed by Zohren and Gill \cite{Zohren2008} described in the introduction, and find that the threshold probabilities are identical to those shown in Fig. \ref{threshold}.  Second, we have looked at the alternative entangled states
\begin{equation}
|\Psi\rangle_{\subs{app}} = \frac{1}{\mathcal{N}} \sum_{j=0}^{d-1} \frac{1}{\sqrt{(j+1)(d-j)}} |j,j\rangle,
\end{equation}
where $\mathcal{N}$ is a normalization factor; these states can achieve a higher degree of violation than the maximally entangled states \cite{Chen2006}, and approximate those states with maximal violation.  While the thresholds are slightly different for these states, they exhibit the same general behavior seen in Fig. \ref{threshold}.

Finally, we have considered another variation of the entangled state
\begin{equation}
|\Psi\rangle_{\subs{rev}} = \frac{1}{\sqrt{d}} \sum_{j=0}^{d-1} |j, d-j\rangle,
\end{equation}
in which Bob's qudit states have been reversed (Bob's unitaries $U_B,b$ must also be ``reversed'').  These states are easier to produce in superconducting circuits \cite{Sharma2015}, and one might think they would be less sensitive to amplitude damping.  However, the thresholds are again only slightly different for this alternative procedure.  In short, the behavior seen in Fig. \ref{threshold} appears to be generic for these inequalities and noise models.

\section{Conclusion}

We have studied the effects of noise on Bell inequality experiments using the higher-dimensional inequalites proposed by Collins, Gisin, Linden, Massar, and Popescu \cite{Collins2002} and Zohren and Gill \cite{Zohren2008} and variations thereof.  By modeling the required operations needed in the preparation, rotation, and single-shot measurement stages of these experiments, we have analyzed how the number of operations scale with the qudit dimension.  For most types of noisy operations, namely amplitude damping, depolarizing, and dephasing noise, we find that the higher-dimensional inequalities require increasingly higher-fidelity operations.  This conclusion runs counter to previous work, which considered depolarizing noise with a fixed noise probability, independent of qudit dimension.  However, this conclusion is in agreement with a very recent analysis of higher-dimensional Bell inequalities with random {\em coherent} errors \cite{Weiss2015}.

We conclude with a few observations.  Our results show that these higher-dimensional inequalities are not more robust against noise when taking into account the full experimental procedure.  Thus, these inequalities do not appear to offer a quick route towards closing the detection loophole with photons.  However, there are other multi-setting ($m>2$) inequalities with qudits \cite{Vertesi2010} that allow for reduced detection efficiencies for both atom-photon and photon-photon entanglement, and in the presence of depolarizing noise.  Our results motivate continued analysis of these and other inequalities (including those with $n>2$) in the presence of realistic noise appropriate for matter and photonic qudits.  The continued acquisition of evidence for nonlocality from a wide variety of physical systems remains an intriguing and important goal of modern quantum physics.

\acknowledgements
This work was supported by the NSF under Project Nos. PHY-1005571 and PHY-1212413.


\section{Introduction}

Hydrodynamics is considered  now as the basic  part of a spatiotemporal
picture of the matter evolution in the processes of ultrarelativistic heavy ion collisions (see recent reviews \cite{Kodama,Song}). To complete the
description of \emph{A+A} collision processes, hydrodynamics must
be supplied with a generator of an initial non-equilibrated state, pre-thermal dynamics which forms the initial near locally equilibrated conditions for hydro-evolution, and prescription for particle production during the breakup of the
continuous medium at the final stage of the matter expansion.

As for the initial state, since it fluctuates on an event-by-event basis, Monte Carlo event generators
are widely used to simulate it in relativistic
$A+A$ collisions. The most commonly used models of initial state
are the MC-Glauber (Monte Carlo Glauber) \cite{MCG}, MC-KLN (Monte Carlo
Kharzeev-Levin-Nardi) \cite{MC-KLN}, EPOS (parton-based Gribov-Regge model) \cite{EPOS}, EKRT (perturbative QCD + saturation model) \cite{EKRT, EKRT-2}, and IP-Glasma (impact parameter dependent
glasma) \cite{IPG}.\footnote{Note that the last two models are able to
reproduce correctly the centrality systematics of the $\delta v_2$ \cite{EKRT-2,IPG-2}.
} The last model also includes non-trivial non-equilibrium
dynamics of the gluon fields which, however, does not lead to a proper
equilibration. To apply these models to data description some
thermalization process has to be assumed. Evidently, in order to reduce
uncertainties of results obtained by means of hydrodynamical models,
one needs to convert a far-from-equilibrium initial state of matter in a
nucleus-nucleus collision to a close to locally equilibrated one by means
of a reasonable pre-equilibrium dynamics. A relaxation time method \cite{preth}, initially developed for the post-hydrodynamic stage, is modified for applications to pre-thermal dynamics.

As for the post-hydrodynamic (afterburner) stage, good progress was made in understanding and
modeling the breakup of hadron matter into hadron gas at the late stage of expansion. Typically, particlization, i.e., the transition from the hadron (or quark-gluon) fluid
to hadronic gas expansions,  is described by means of the so-called Cooper-Frye prescription within hybrid
models. In such an approach, the sudden conversion of the fluid to particles happens
 at a hypersurface of hadronization or chemical
freeze-out. It has long been known that such a
matching prescription has problems with the energy-momentum
conservation laws when fluid is converted to particles at the
hypersurface which contains non-space-like parts (see up-to-date review \cite{Kodama} and references therein). These problems can
be avoided by using the HydroKinetic Model (HKM) that was proposed in Ref.
\cite{hydrokin-1} and further developed in Refs. \cite{hydrokin-2, hydrokin-2a}
(see also Ref. \cite{hydrokin-3}). In basic HKM, particlization is considered not as a sudden, as in hybrid models, but as a realistic continuous process.\footnote{For recent
discussions of the particlization procedure see, e.g., Refs.
\cite{hydrokin-3,convert}).} In such an approach the description of continuous particle emission from the fluid  is based not on the distribution function  but on the so-called
escape function \cite{hydrokin-1, hydrokin-3}.

In Refs. \cite{preth, Naboka} the method of escape probabilities, similar to that used in HKM formalism for post-thermal dynamics, is applied to the pre-thermal one. It makes it possible to describe the pre-thermal evolution of the energy-momentum tensor using three free parameters: the time of the initial (non-equilibrium) state formation, mean relaxation time, and time of thermalization, say, 1 fm/c. This phenomenological model allows one to provide  initialization of the  hydrodynamical evolution, and is a more general and simple alternative to the so-called ``anisotropic hydrodynamics'' \cite{RSF, aniz, Mart-1, Mart-2, Mart-3}, as discussed in detail in \cite{Naboka}.

The target energy-momentum tensor of the pre-thermal dynamics of \emph{A+A} collisions, which the initially created  system reaches during evolution, is the tensor of viscous relativistic hydrodynamics. The further evolution is governed by the equations of the Israel-Stewart framework of relativistic viscous hydrodynamics until the lowest possible temperature when the system is still close to local thermal and chemical equilibrium. It defines the hypersurface which is often called the chemical freeze-out isotherm.
The particlization procedure starts typically at this very hypersurface. The hypersurface contains non-space-like parts that formally enforce one to use the HKM continuous particlization procedure, to eliminate the problem with violation of the energy-momentum conservation. However, with almost smooth (averaged) initial conditions, the effects originated from these parts are quite small and do not affect significantly  the observables  \cite{hydrokin-2a}.\footnote{The reason is that at small proper times $\tau$ the three-dimensional areas of non-space-like parts are small, and for large times the hydro-velocities at the corresponding parts are large.} Moreover, in viscous hydrodynamics the effects of violation of local equilibrium, which is produced in HKM for an ideal fluid, are partially taken into account. Therefore, summarizing all that and accounting for very time consuming calculations in HKM, we will use in this paper the sudden transition from hadron matter to hadron gas, similar to the hybrid models. In future studies we plan to include the option of continuous HKM particlization for situations when undesirable  contributions from non-space-like parts are fairly large: it can be for  $4\pi$- geometry observables and event-by-event analysis. Considering also that an important component of the full model is the pre-thermal evolution where the description is built on the same idea as the basic HKM formalism, we will call the complete model, which includes generation of the initial state, its thermalization, sudden (with the possible option of ``continuous'') particlization, subsequent ultrarelativistic QMD (UrQMD) hadronic  cascade \cite{urqmd}, and one- and multi- particle spectra formation, the integrated HydroKinetic Model (iHKM).

The paper is organized as follows. Section 2 includes subsections that describe the ingredients of the complete iHKM. Section 3 is devoted to the results and discussions, while the Section 4 summarizes the physical meaning of obtained results.

\section{The model description}
As discussed in the Introduction the integrated HydroKinetic Model (iHKM) consists of five ingredients which we describe below.

\subsection{Initial state}
Though prethermal dynamics is used to initialize hydrodynamics, it itself needs the initial state of the matter at the starting time $\tau_0$, say $\tau_0=0.1$ or 0.5 fm/c, when one can speak about energy density distribution in non equilibrated matter. Following Ref. \cite{{Krasnitz}}, an appropriate scale
controlling the formation of gluons with a physically well-defined energy is roughly $3/g^2\mu$, where $g^2 = 4\pi\alpha_s$ and $\mu^2$ are dimensionless parameters, which
is the variance of the Gaussian weight over the color charges
$\rho$ of partons. The estimate of this time gives $\tau_0=0.1$ fm/c for the top energy available at the BNL Relativistic heavy Ion Collider (RHIC) and $\tau_0=0.07$ fm/c for the CERN Large Hadron Collider (LHC) full energy \cite{Lappi}. We will utilize  $\tau_0=0.1$ fm/c as the default value for the LHC energy $\sqrt{s_{NN}}= 2.76$ TeV; however, we will use also $\tau_0=0.5$ fm/c to analyze the corresponding variation of the results.

We attribute the initial state to the time $\tau_0$. The generation of the initial state is based on the GLISSANDO 2 \cite{generator, gliss2}  package. This package works in the frame of the semiclassical Glauber model. Within this approach at the very initial stage of the collision, individual interactions between nucleons deposit transverse energy.  Each deposition of the transverse energy at a certain space-time point or region is called a source and each source has its weight which is called relative deposited strength, or RDS. The normalization of RDS can be treated as an insertion of the additional parameter. We choose its value in such a way that all charged-particle multiplicity distribution, obtained in iHKM for 5\% most central collisions, fits its experimental value. The RDS may be different for the wounded nucleons and the binary collisions, and it can fluctuate from source to source. We use a mixed model, amending the wounded nucleon model with some binary collisions. In this model, a wounded nucleon obtains the RDS of $(1-\alpha)/2$ and a binary collision has the RDS of $\alpha$. The total RDS, averaged over events, is then $(1-\alpha)N_W/2+\alpha N_{bin}$. Of course, the result of simulations fluctuates from event to event. In this paper we do not study fluctuations of
observables, and will simulate inclusive observables basing on the approximation of the averaged initial state for each centrality class. This approximation implies that for each selected centrality only the one hydro run is performed.
 The simulation of a single event is made in three stages:
\begin{enumerate}
  \item Generation of the positions of nucleons in the two colliding nuclei according to the fluctuating nuclear density distribution. The form of the mean distribution depends on the mass of nuclei $A$. For sufficiently large nuclei this distribution has a Woods-Saxon form with an account of nuclear deformations, the latter are small in our case of Pb+Pb collisions.
 \item Generation of the transverse positions of the sources and their RDSs.
  \item Calculation of the physical quantities and writing the results in the output file.
\end{enumerate}
In order to obtain the initial conditions from the RDS, one can put the RDS to be proportional to energy or entropy. We choose it to be proportional to energy because entropy is not created yet in a non-equilibrium initial state. Then, for the initial state averaged over many events with a chosen centrality bin, one can write ($b$ is the impact parameter associated with centrality)
\begin{equation}
\epsilon(b;\tau_0, {\bf r}_T) = \epsilon_0 \frac{(1-\alpha)N_W({b,\bf r}_T)/2+\alpha N_{bin}(b,{\bf r}_T)}{(1-\alpha)N_W({b=0,\bf r}_T=0)/2+\alpha N_{bin}(b=0,{\bf r}_T=0)},
\label{RDS}
\end{equation}
where $\epsilon_0$ is a fitting parameter defining the maximal energy density of the initial state that is reached at minimal centrality, say 0--5\% ($b\approx 0$), at ${\bf r}_T=0$. The value $\epsilon_0$ and $\alpha$ are the same for all centrality classes. The values of $\epsilon_0$ and $\alpha$ are fixed by the  reproduction of the charged-particle multiplicity at different centralities.

\subsection{Pre-thermal stage of the matter evolution}
The generation of the initial energy density configuration with the help of GLISSANDO is supposed to be followed by a thermalization process.  As we have already discussed in the Introduction we describe the matter evolution at the pre-thermal stage with the relaxation model \cite{preth, Naboka}. We assume a longitudinal boost invariance, which is a good approximation for the central rapidity region of the fireball at high collision energies. The relaxation model simulations are applied to the system for the early stage of the heavy ion collision, from the proper time $\tau=\tau_0$ to $\tau=\tau_{th}$, where $\tau=\sqrt{t^2-z^2}$. As mentioned above, we use $\tau_0=0.1$ fm/c as the basic value and $\tau_0=0.5$ fm/c for a comparison. Thermalization time is supposed to be $\tau_{th}=1.0$ fm/c, the``conventional'' time of the formation of strongly interaction matter, $t\approx 1/m_{\pi}.$\footnote{Note, that calculations with $\tau_{th}=1.5$ fm/c and with the other parameters to be the same as in our basic scenario (see later) lead to very similar results as calculations with $\tau_{th}=1.0$ fm/c. } The total energy-momentum tensor of the matter is taken in the form \cite{preth, Naboka}
\begin{eqnarray}
T^{\mu \nu}(x)=T^{\mu \nu}_{\text{free}}(x){\cal
P}(\tau)+T_{\text{hydro}}^{\mu \nu}(x)[1-{\cal P}(\tau)], \label{1}
\end{eqnarray}
where $T^{\mu \nu}_{\text{hydro}}(x)$ and $T^{\mu \nu}_{\text{free}}(x)$ are hydrodynamic (local equilibrium) and free (fully or almost free evolving)  components of the energy-momentum tensor, $\cal P(\tau)$ is the weight function satisfying the following conditions: $0 \leq {\cal P}(\tau) \leq 1$, ${\cal P}(\tau_0)=1$, ${\cal P}(\tau_{th})=0$, and $\partial_{\mu}{\cal P}(\tau)_{\tau_{\text{th}}}=0$. Its form will be discussed below. As for the free (or almost free) evolving component, it is well defined in such models as IP-Glasma and EPOS, but not in MC Glauber. To analyze different kinds of anisotropy of the initial state we present the boost-invariant distribution function at the initial hypersurface $\sigma_0$: $\tau=\tau_0$ in the factorized form
\begin{equation}
f(t_{\sigma_0},\textbf{r}_{\sigma_0},\textbf{p}) = \epsilon(b;\tau_0, {\bf r}_T)f_0(p)
\label{f0}
\end{equation}
where $\epsilon(b;\tau_0, {\bf r}_T)$ is defined by Eq. (\ref{RDS}).   At other times $f(x,p)$ is defined by the free-streaming requirement $f(t,\textbf{r},\textbf{p})=f(t_{\sigma_0},{\textbf{r}}-\frac{{\textbf{p}}}{p_{0}} (t-t_{\sigma_0}),\textbf{p})$, so that at all times the free evolution of the energy-momentum tensor is defined by the formula
\begin{eqnarray}
T^{\mu \nu}_{\text{free}}(x)=\int d^{3}p \frac{p^{\mu}p^{\nu}}{p_{0}}f(x,p). \label{2}
\end{eqnarray}
Then the evolution of $T^{\mu \nu}_{\text{free}}(x)$ and corresponding distribution function $f(x,p)$ are defined by the initial conditions, which are generated by GLISSANDO, and the initial momentum anisotropy of the function $f_0(p)$. For ${\cal P}(\tau_0)=1$, we have $T^{\mu \nu}_{\text{total}}(\tau_0,x)=T^{\mu \nu}_{\text{free}}(\tau_0,x)$.

 The hydrodynamic component of the energy-momentum tensor has the form
\begin{eqnarray}
T^{\mu \nu}_{\text{hydro}}(x)= [\epsilon_{\text{hydro}}(x) +
p_{\text{hydro}}(x)+\Pi]u^{\mu}_{\text{hydro}}(x)u^{\nu}_{\text{hydro}}(x)
- [p_{\text{hydro}}(x)+\Pi]g^{\mu \nu} +\pi^{\mu \nu},
 \label{6}
\end{eqnarray}
where $\epsilon_{hydro}$ is the energy density in the local rest frame, $p_{hydro}$ is the pressure, $\pi^{\mu \nu}$ is the shear stress tensor, $\Pi$ is the bulk pressure, and $u^{\mu}_{\text{hydro}}(x)$ is the energy flow four-vector. We neglect the bulk pressure, $\Pi=0$. In curvilinear (hyperbolic) coordinates the equation of motion for the shear stress tensor is taken in the form
\begin{eqnarray}
\langle u^\gamma \partial_{;\gamma} \pi^{\mu\nu}\rangle
=-\frac{\pi^{\mu\nu}-\pi_\text{NS}^{\mu\nu}}{\tau_\pi}-\frac 4 3
\pi^{\mu\nu}\partial_{;\gamma}u^\gamma, \label{evolution-1}
\end{eqnarray}
where the semicolon means the covariant derivation, brackets are defined as $\langle A^{\mu\nu}\rangle=(\frac 1 2 \Delta^\mu_\alpha
\Delta^\nu_\beta+\frac 1 2 \Delta^\nu_\alpha \Delta^\mu_\beta -
\frac 1 3 \Delta^{\mu\nu}\Delta_{\alpha\beta})A^{\alpha\beta}$, $\Delta^{\mu\nu}=g^{\mu\nu}-u^\mu u^\nu$, and $\pi_{NS}$ is the Navier-Stokes shear stress tensor:
\begin{eqnarray}
\pi^{\mu\nu}_\text{NS}&=\eta(\Delta^{\mu\lambda}\partial_{;\lambda}u^\nu+\Delta^{\nu\lambda}\partial_{;\lambda}u^\mu)-\frac
2 3 \eta\Delta^{\mu\nu}\partial_{;\lambda}u^\lambda.
\label{evolution-2}
\end{eqnarray}
For $\cal P(\tau)$ we use the ansatz \cite{preth, Naboka}
\begin{eqnarray}
{\cal P}(\tau)=  \left (
\frac{\tau_{\text{th}}-\tau}{\tau_{\text{th}}-\tau_0}\right
)^{\frac{\tau_{\text{th}}-\tau_0}
 {\tau_{\text{rel}} }},
 \label{3}
\end{eqnarray}
where, as mentioned above, $\tau_0=0.1$ fm/c, $\tau_{th}=1.0$ fm/c, and $\tau_{rel} \leq \tau_{th}-\tau_0$. The value of $\tau_{rel}$ is the model parameter, which characterizes the rate of thermalisation, that is, the speed of conversion of the non-equilibrium state to about the equilibrated one. Writing down the conservation laws for the total energy-momentum tensor in the form $\partial_{;\mu}T^{\mu\nu}_{total}=0$ and accounting that for the free streaming $\partial_{;\mu}T^{\mu\nu}_{free}=0$, we have
\begin{eqnarray}
\partial_{;\mu}\{[1-{\cal P}(\tau)]T^{\mu
\nu}_{\text{hydro}}(x)\}= - T^{\mu
\nu}_{\text{free}}(x)\partial_{;\mu}{\cal P}(\tau). \label{7}
\end{eqnarray}
Let us introduce the re-scaled hydrodynamic tensor $\widetilde{T}^{\mu \nu}_{\text{hydro}}(x)$ with initial conditions $\widetilde{T}^{\mu \nu}_{\text{hydro}}(x)=0$ at $\tau = \tau_{0}$ for all $x$. Then
\begin{eqnarray}
\partial_{;\mu}\widetilde{T}^{\mu
\nu}_{\text{hydro}}(x)= - T^{\mu \nu}_{\text{free}}(x)\partial_{;\mu}
{\cal P}(\tau). \label{8}
\end{eqnarray}
This is the hydrodynamic-type equation of motion with a source. As mentioned above, $T^{\mu \nu}_{free}$ is defined by the the initial state at proper time $\tau_0$, and ${\cal P}(\tau)$ is defined explicitly, then the source in the right-hand side of Eq. (\ref{8}) can be calculated for all $\tau$ and $x$. The introduction of the re-scaled hydrodynamic tensor leads to the rescaling of the shear stress tensor  $\widetilde{\pi}^{\mu\nu} =\pi^{\mu\nu}(1-{\cal P})$. Then, multiplying Eq. (\ref{evolution-1}) by $(1-{\cal P})$, we have the equation of motion for the re-scaled shear stress tensor
\begin{eqnarray}
[1-{\cal P}(\tau)]\left \langle u^\gamma \partial_{;\gamma}
\frac{\widetilde{\pi}^{\mu\nu}}{(1-{\cal P}(\tau))}\right \rangle
=-\frac{\widetilde{\pi}^{\mu\nu}-[1-{\cal
P}(\tau)]\pi_\text{NS}^{\mu\nu}}{\tau_\pi}-\frac {4}{3}
\widetilde{\pi}^{\mu\nu}\partial_{;\gamma}u^\gamma.
\label{evolution-3}
\end{eqnarray}

We also have to specify the equation of state in order to close the set of evolutionary equations. In all calculations the Laine-Schroeder equation of state (EOS) \cite{Laine} is applied.
At $\tau = \tau_{th}$ we switch to the viscous hydrodynamic approach. The equations of motion are nearly the same as in Eq. (\ref{8}), but with the right-hand side equal to zero.

\subsection{Matter evolution in thermal and chemical locally near-equilibrated zone}
At $\tau=\tau_{th}=1$ fm/c the ${\cal P}(\tau_{th})=0$ and the target function is reached: $T^{\mu \nu}(x) = T^{\mu\nu}_{hydro}(x)$. The further evolution is described by the relativistic viscous hydrodynamics according to the equations  (\ref{6}) - (\ref{evolution-2}), (\ref{8}), and (\ref{evolution-3}) with ${\cal P}(\tau)\equiv 0$. A numerical solution of the viscous hydrodynamic equations is constructed with the \texttt{vHLLE} code \cite{vhlle}. Such an evolution describes the expansion of superdense quark-gluon and hadron matter close to local chemical and thermal equilibria with a baryon chemical potential $\mu_B=0$ (which is a good approximation for LHC energies) until the temperature when such an approach breaks down. Then the system has lost the properties of local equilibrium, thermal and chemical as well, and another approximation should be used.

\subsection{Particlization stage}

As discussed in the Introduction, the basic HKM describes particlization as  a continuous process (for review see \cite{hydrokin-3}): particles gradually ``escape'' from the expanding fluid, forming the non-equilibrium Wigner function, which then can be used at some space-like hypersurface as the input for UrQMD hadronic cascade  \cite{hydrokin-2a}. As it was already discussed, at smooth initial conditions the results are very similar to a sudden particlization scenario and we will use the latter in this paper. We assume that the chemically and thermally locally equilibrated evolution takes place until temperature $T=165$ MeV (corresponding to an energy density $\epsilon=0.5$ GeV/fm$^3$ for the Laine-Schroeder EOS) is reached, and switch to particle cascade at the hypersurface defined by this criterion. Such a switching surface is built during the hydrodynamic evolution with the help of the Cornelius \cite{convert} routine. We apply the Cooper-Frye formula to convert the fluid to the cascade of particles:
\begin{eqnarray}
p^0 \frac{d^3 N_i(x)}{d^3 p}=d\sigma_{\mu}p^{\mu}f(p\cdot u(x),T(x),\mu_i(x))
\label{Cooper-Frye1}
\end{eqnarray}
In order to take into account the viscous corrections to the distribution function, Grad's 14-moment ansatz is used. We assume that the corrections are the same for all hadron species. Then formula (\ref{Cooper-Frye1}) transforms to
\begin{eqnarray}
\frac{d^3\Delta N_i}{dp^{*}d(cos\theta)d\phi}=\frac{\Delta \sigma^{*}_{\mu}p^{*\mu}}{p^{*0}}p^{*2}f_{eq}\left(p^{*0};T,\mu_i\right)\left[1+(1\mp f_{eq})\frac{p^*_\mu p^*_\nu \pi^{*\mu \nu}}{2T^2(\epsilon+p)}\right]
\label{Cooper-Frye2}
\end{eqnarray}
This distribution function is used to create an ensemble of particles on the hypersurface. At first, the average number of hadrons of every sort is calculated:
\begin{equation}
\Delta N_i = \Delta \sigma_{\mu}u^{\mu}n_{i,th}=\Delta \sigma ^*_0 n_{i,th}
\label{av1}
\end{equation}
The average total number of particles then equals $\langle N_{tot} \rangle=\Sigma_i N_i$. The exact total number of the particles to be created is sampled according to Poisson distribution with a mean value $<N_{tot}>$. The type of each generated particle is chosen randomly based on probabilities $N_i/N_{tot}$. Then the momentum is assigned to each particle in the local rest frame of the fluid. The direction of momentum is chosen randomly in the $4\pi$ solid angle and its modulus is generated according to the isotropic part of Eq. (\ref{Cooper-Frye2}). After that, the corrections for $W_{residual}\cdot W_{visc}$ are applied via the rejection sampling. The particle position is set to be equal to the position of centroid of the surface element, and the spacetime rapidity is sampled randomly within the longitudinal size of the volume element. Finally, the particle momentum is Lorentz boosted to the center-of-mass frame of the fireball.

\subsection{Hadronic cascade}
The generated hadrons are then fed into the UrQMD cascade. Since the cascade accepts only a list of particles at an equal Cartesian time as an input, the created particles are propagated backward in time to when the first particle was created. The particles are not allowed to interact in the cascade until their trajectories cross the particlization hypersurface. The Laine-Schroeder EoS, that is applied in our analysis corresponds to an equilibrium hadron-resonance gas consisting of about 360 hadrons in the low temperature limit. Many of those heavy hadrons are not included in the UrQMD hadron list. To prevent violation of the energy-momentum conservation even at space-like parts of the isotherm $T=165$ MeV, we decay the heavy resonances, which are not in the UrQMD list, just at the switching hypersurface. The particle propagation is stopped at Cartesian time 400 fm/c, where their coordinates and momenta are recorded. We generate 50000 GLISSANDO events for each centrality class to produce averaged initial energy density profile for the pre-thermal/hydrodynamic evolution. Only the one hydrodynamic calculation is performed for each centrality class with a given set of parameters. For each hydrodynamic evolution 20000 UrQMD events are generated. The generated sets of events are stored in \texttt{ROOT} trees and files, which are further processed with scripts to plot momentum distributions and calculate flow coefficients or correlation functions for physical analysis.

\section{Results and discussion}

To perform calculations of the spectra, anisotropic flow coefficients, and interferometry radii we need to specify the initial states and model parameters. As discussed above, we generate these initial conditions with GLISSANDO 2, a Monte Carlo Glauber generator. The output of this generator is the relative deposited strength (RDS), which defines the initial energy density (\ref{RDS}) through  the number of wounded nucleons $N_W$ and number of binary collisions $N_{bin}$. The relative weight of these quantities $\alpha$ is crucial for the description of all charged-hadron multiplicities at different centralities. It is found that the value $\alpha = 0.24$ at different reasonable values of other parameters  gives the  best fit for the multiplicities at LHC energy $\sqrt{s}=2.76$ TeV. Figure \ref{fig:mult} demonstrates a typical fit at this value of $\alpha$.

\begin{figure}[H]
     \centering
     \includegraphics[width=0.55\textwidth]{mult_bw.eps}
     \caption{\scriptsize
Multiplicity centrality dependence for the iHKM basic scenario with the parameter values $\tau_0=0.1$ fm/c, $\tau_{rel}=0.25$ fm/c, $\eta/s=0.08$, and $\Lambda=100$ for the LHC energy $\sqrt{s_{NN}}=2.76$ TeV. Experimental data are from \cite{Alica_mult}.}
\label{fig:mult}
\end{figure}

To analyze the effects of momentum anisotropy of the initial state, similarly to Ref. \cite{Naboka} we choose the momentum dependence of the distribution function $f_0(p)$ in Eq. (\ref{f0}) in the form \cite{Rub}:
\begin{eqnarray}
f_0(p)=g \exp\left(-\sqrt{\frac{(p\cdot U)^2-(p\cdot
V)^2}{\lambda_{\perp}^2}+\frac{(p\cdot
V)^2}{\lambda_{\parallel}^2}}\right),
\label{anis1}
\end{eqnarray}
where $U^{\mu}=( \cosh\eta, 0, 0, \sinh\eta)$, $V^{\mu}=(\sinh\eta, 0, 0,\cosh\eta)$. Note, that in the rest frame of the fluid element, $\eta=0$, $(p\cdot U)^2-(p\cdot
V)^2=p_{\perp}^2$ and $(p\cdot
V)^2=p_{\parallel}^2$. Then one can see that the parameters $\lambda_{\parallel}^2$ and $\lambda_{\perp}^2$ can be associated with the different temperatures along the beam axis and orthogonal to it correspondingly. The parameter $\Lambda=\lambda_{\perp}/\lambda_{\parallel}$ is the main parameter, that defines  momentum anisotropy of the initial state.  We use two values for it: $\Lambda=1$ (momentum isotropic case) and $\Lambda=100$ (anisotropic case, almost no pressure in the longitudinal direction). We put $\lambda_{\perp}= 1.4$ GeV \cite{Naboka}. At $m^2/\lambda^2_{\perp}<<1$ the last value is just absorbed into a factor $g=g(\lambda_{\perp},\lambda_{\parallel})$ that is the normalization constant in Eq. (\ref{anis1}). This constant is defined from the requirement that the initial energy-momentum tensor zero component, formally associated with the dimensional function $f_0(p)$ only and so calculated according to Eq. (\ref{2}) with substitution $f(x,p) \rightarrow f_0(p)$, is unity:   $T^{00}[f_0(p)]=1$ [see Eqs. (\ref{RDS}) and (\ref{f0})]. As for the anisotropy factor $\Lambda$, it also does not affect the initial  energy density distribution (\ref{RDS}), but reveals itself in the process of the matter evolution at the pre-thermal stage, see Sec. II B.

The generated distribution functions are related to the two initial times  $\tau_0=0.1$ and $\tau_0=0.5$ fm/c for comparison. The thermalization time is chosen to be, as discussed before, $\tau_{th}=1$ fm/c. When other parameters are fixed,  we make calculations for 0--5\%, 5--10\%, 10--20\%, 20--30\%, 30--40\%, 40--50\%, and 50--60\% centralities with the same normalization coefficient $\epsilon_0$ in Eq. (\ref{RDS}) and fixed relative contribution of binary collision $\alpha = 0.24$. The latter, as we discussed above, does not in fact depend, on all the other parameters, contrary to normalization $\epsilon_0$. The $\epsilon_0$ values at different scenarios, that are associated with iHKM parameters, are shown in Table 1.

\begin{center}
  \begin{tabular}{| c | c | c | c | c | c | c |}
   \hline
   model & $\Lambda$ & $\tau_{rel}$ & $\eta/S$ & $\tau_0$ & $\langle \frac{\chi^2}{ndf} \rangle$ & $\epsilon_0$ [GeV/fm$^3$] \\ \hline
   hydro & - & - & 0 & 0.1 & 5.16 & 1076.5 \\  \hline
   hydro & - & - & 0.08 & 0.1 & 6.93 & 738.8 \\  \hline
   iHKM & 1 & 0.25 & 0.08 & 0.1 & 3.35 & 799.5 \\  \hline
   iHKM & 100 & 0.25 & 0.08 & 0.1 & 3.68 & 678.8 \\ \hline
   iHKM & 100 & 0.75 & 0.08 & 0.1 & 3.52 & 616.5 \\ \hline
   iHKM & 100 & 0.25 & 0.2 & 0.1 & 6.61 & 596.9 \\ \hline
   iHKM & 100 & 0.25 & 0.08 & 0.5 & 5.36 & 126.7 \\ \hline
	\end{tabular}
\end{center}
	Table 1. The maximal initial energy densities $\epsilon_0$ and the mean $\chi^2/ndf$ for different scenarios. The values $\tau_0$, $\tau_{rel}$ correspond to fm/$c$. The mean $\chi^2/$ndf is calculated by an averaging $\chi^2/$ndf values taken from model comparison with observed pion, soft pion, kaon, and antiproton spectra, mid-rapidity charged multiplicity dependence on centrality, $v_2$ transverse momentum dependence for charged particles, and the pion interferometry radii as functions of the transverse momentum.

For each set of parameters with corresponding $\epsilon_0$ (Table 1),  the experimental multiplicity dependence on centrality is reproduced well, similar to Fig. \ref{fig:mult}, with the same $\alpha = 0.24$. There is a  correspondence between the ratios of $\epsilon_0$ at the different parameters and inverse  ratios of the energy densities at the thermalization time $\tau_{th}=1$ fm/c, obtained at the such parameters in \cite{Naboka} when the initial energy densities at $\tau_0$ are the same.

After specifying the initial conditions we run the relaxation model, which describes the pre-thermal stage with different relaxation times $\tau_{rel}=0.25$ and 0.75 fm/c.  The grid spacing and time step for the relaxation model and further pure hydrodynamic calculations are $\delta x = \delta y = 0.2$ fm/c and $\delta \tau = 0.05$ fm/c. Here we perform 2+1-dimensional longitudinally boost-invariant calculations in iHKM, which is a good approximation for the central rapidity interval at LHC energies. At the  thermalization (proper) time  $\tau_{th}=1.0$~fm/c, the evolutionary equations at the pre-thermal stage are smoothly switched (automatically) to relativistic viscous hydrodynamic equations in the Israel-Stewart framework.

We use the two values of the shear viscosity to entropy ratio in the hydrodynamic phase: the minimal one $\eta/s=0.08\approx \frac{1}{4\pi}$ and $\eta/s=0.2$ for comparison. Also we compare the basic scenario (BS) with $\tau_0=0.1$, $\tau_{rel}=0.25$ fm/c, $\eta/s=0.08$, anisotropy parameter $\Lambda=100$, which is selected as giving the optimal description of the experimental data (see below),  with results of calculations with other different parameters, including  viscous and ideal pure hydrodynamic scenarios (without the pre-thermal stage but with the subsequent hadronic cascade).

\begin{figure}
     \centering
     \includegraphics[width=0.8\textwidth]{pion_arxiv.eps}
     \caption{\scriptsize
Resulting pion spectra in the $0.1 <p_T<3$ GeV/\emph{c} region for various centrality classes obtained in the iHKM basic scenario (as in Fig. \ref{fig:mult}). The results are compared with those in iHKM at the other parameter $\tau_0=0.5$ fm/c and with  pure viscous hydro at the starting time $\tau_{th}\rightarrow\tau_0=0.1$ fm/c for centrality classes 0--5\%, 10--20\%, and 30--40\%. The experimental data are from \cite{Alice_spectra}. The spectra for different centralities are multiplied by a factor of 2 ($2^6=64$ for 0--5\% centrality). The lines correspond to different centrality classes as marked by letters: a, 0--5\%; b, 5--10\%; c, 10--20\%; d, 20--30\%, e, 30--40\%; f, 40--50\%; g, 50--60\%}
\label{fig:pi1}
\end{figure}

\begin{figure}
     \centering
     \includegraphics[width=0.8\textwidth]{pion2_arxiv.eps}
     \caption{\footnotesize
 Detailed picture of  pion spectra in soft $p_T$ region for 0--5\% centrality in the iHKM basic scenario (as in Fig. \ref{fig:mult}) in comparison with the results obtained with (1) the other relaxation time $\tau_{rel}=0.75$ fm/c instead of 0.25 fm/c, or with (2) the isotropic parameter  $\Lambda=1$ instead of the anisotropy one $\Lambda=100$,  or with (3) the other initial time $\tau_0=0.5$ fm/c instead of 0.1 fm/c. Also the results for pure viscous hydro, starting at $\tau_{th}\rightarrow\tau_0=0.1$ fm/c are presented. The experimental data are from \cite{Alice_spectra}.}
\label{fig:pi2}
\end{figure}

\begin{figure}
     \centering
     \includegraphics[width=0.85\textwidth]{kaon_arxiv.eps}
     \caption{\footnotesize
The resulting kaon spectra under the same conditions as in Fig. \ref{fig:pi1}. }
\label{fig:k1}
\end{figure}

\begin{figure}
     \centering
     \includegraphics[width=0.85\textwidth]{aproton_arxiv.eps}
     \caption{\footnotesize
The resulting antiproton spectra under the same conditions as in Fig.
 \ref{fig:pi1}. }
\label{fig:aproton1}
\end{figure}

In Fig. \ref{fig:pi1}  we compare the results with the pion spectra at $0.1<p_T<3$ GeV/\emph{c}.  As one can see at fairly large $p_T > 1.5$ GeV/\emph{c} the results for pion transverse spectra in the scenario with initial time $\tau_0= 0.5$ fm/c describe poorly the experimental data: $\chi^2_{pion}$/ndf=14.33. The best plots correspond to shear viscosity value $\eta/s=0.2$, $\chi^2_{pion}$/ndf=2.48, in the basic scenario (BS) with $\eta/s=0.08 \rightarrow 0.2$. The BS itself leads to satisfactory results for pion spectra also. According to Table 1 the basic scenario leads to mean value $\langle \frac{\chi^2}{\tt{ndf}}\rangle$=3.68 over all  observables under consideration. while the scenario with $\eta/s=0.08 \rightarrow 0.2$ has $\langle \frac{\chi^2}{\tt{ndf}} \rangle$= 6.61 because it does not describe $v_2$ (see below Fig. \ref{fig:v2}), $\chi^2_{v_2}$/ndf=27.92.

In addition to the basic scenario the relatively good description of the wide class of observables gives, according to Table 1, the similar scenarios with, however, larger relaxation time, $\tau_{rel}=0.25$ fm/c $\rightarrow 0.75$ fm/c, and the scenario with isotropic initial state, $\Lambda = 100 \rightarrow 1$. So, for more detailed analysis, let us consider the interesting phenomenon of soft pion enhancement presented in Fig. \ref{fig:pi2}. As one can see, the iHKM scenarios with isotropic initial state, $\Lambda=1$ or with larger relaxation time, $\tau_{rel}=0.75$, obviously underestimate  $p_T$ spectra in the soft region, namely,  $\chi^2_{\tt {soft pions}}$/ndf=1.03 for the isotropic case and $\chi^2_{\tt {soft pions}}$/ndf=1.27 for $\tau_{rel}=$ 0.75 fm/c; pure (viscous) hydrodynamics scenario with $\tau_{th}=\tau_0= 0.1$ fm/c is also less effective at soft $p_T$ for pions. At the same time,  the basic scenario leads to twice less value $\chi^2_{\tt {soft pions}}$/ndf=0.59. The description of this soft momentum region becomes the best when the initial state formation is ascribed to the later time $\tau_0=0.5$ fm/c. However, it excludes, as we discussed above (see also Table 1), a satisfactory spectra description at $p_T>1.5$ GeV/\emph{c} despite that there is no problem with $v_2$ coefficients and the interferometry radii.  This could mean the necessity of including mini-jets in the hydrokinetic picture in order to  describe the large $p_T$ spectra  simultaneously with very small transverse momenta. One of the other proposals considers the specific mechanism of soft pion radiation, namely, the Bose-Einstein condensation in the thermal model \cite{Begun}. However, since the pion spectra in the basic scenario are within the experimental errors with small $\chi^2_{\tt {soft pions}}$/ndf , we consider now the BS as the most realistic, allowing to describe well the bulk observables in wide $p_T$ region.

As one can see from Fig. \ref{fig:k1}, iHKM also describes well the kaon spectra in the basic scenario, $\chi^2_{kaon}$/ndf=1.9.  As for the antiproton spectra (Fig. \ref{fig:aproton1}), one can see that the basic scenario is not the best one (Fig. \ref{fig:aproton1}). That is why, if one compares the mean values of $\langle \frac{\chi^2}{\tt{ndf}}\rangle$ over all observables in Table 1, this scenario brings  slightly worse results than scenarios with an isotropic initial state and anisotropic one with larger relaxation time. At the same time the last two scenarios do not describe the enhancement of soft pions as we already discussed.

Let us apply iHKM to a description of elliptic flows.  The elliptic flows in iHKM are associated with
the reaction plane that is well-defined in GLISSANDO 2. We compare the results with the experimental $v_2$-coefficients   estimated
with four-particle cumulants measured for unidentified charged particles as a function of
transverse momentum for various centrality classes, $v_2\{4\}$ \cite{Alice_v2}.
Figure \ref{fig:v2}  demonstrates all charged-particle $v_2$-coefficients for the basic scenario in comparison with iHKM results with $\eta/s=0.2$ instead of $\eta/s=0.08$ and also with ideal hydro starting at
$\tau_{th}=\tau_0= 0.1$ fm/c. One can see that the two last scenarios lead to poorly described $v_2$-coefficients. At the same time the basic scenario, as well as the iHKM one with $\tau_0=0.5$ instead of $\tau_0=0.1$ fm/c and also the viscous pure hydrodynamics scenario with $\tau_{th}=\tau_0= 0.1$ fm/c describe the data quite satisfactorily, see Fig. \ref{fig:v2_2}. However, the  two last scenarios describe the spectra worse: the one with $\tau_0=0.5$, at fairly large $p_T$; the other, the pure viscous hydro, at both large and small momenta.

\begin{figure}
     \centering
     \includegraphics[width=0.75\textwidth]{v2_arxiv.eps}
     \caption{\scriptsize
  All charged-particle $v_2$ coefficients for centrality classes 0--5\%, 5--10\%, 10--20\%, 20--30\%, 30--40\%, 40--50\% and 50--60\%, obtained in the iHKM basic scenario (as in Fig. \ref{fig:mult}). The results are compared with those in iHKM at the other parameter dissipation condition,  $\eta/s=0.2$ instead of $0.08$ and with  ideal hydro with the starting time $\tau_{th}=\tau_0=0.1$ fm/c for centrality classes 5--10\%, 10--20\%, and 20--30\%. The experimental data are from \cite{Alice_v2}.}
\label{fig:v2}
\end{figure}

\begin{figure}
     \centering
     \includegraphics[width=0.75\textwidth]{v2_2_arxiv.eps}
     \caption{\scriptsize
 Same as Fig. \ref{fig:v2}, but the results are compared with those in iHKM at the other initial time,  $\tau_0=0.5$ fm/c instead of $0.1$ fm/c and with  viscous hydro at the starting time $\tau_{th}=\tau_0=0.1$ fm/c.}
\label{fig:v2_2}
\end{figure}

An important observation for evolutionary models is the interferometry radii that reflect the space-time structure of the particle emission from an expanding fireball. The detailed study within iHKM shows that if a selection of some set of parameter values, such as $\tau_0$, $\tau_{rel}$, $\Lambda$, $\eta/s$,  is accompanied by the renormalization of the maximal initial energy density $\epsilon_0$, similar to that in Table 1, to keep the iHKM charged-particle multiplicity to be  equal to the experimental one, the interferometry radii change within a few percent only, which is a much smaller change than the one caused by a change in the centrality class of the collision. Published results (see, e.g. review \cite{hydrokin-3}) have shown that the interferometry radii depend on particle species, multiplicity, and initial system size. The other details are not so important.

\begin{figure}
     \centering
     \includegraphics[width=0.8\textwidth]{Rout_fin.eps}
     \caption{\small (Color online)
The $R_{out}$ dependence on transverse momentum for different centralities in the iHKM basic scenario under the same conditions as in Fig. \ref{fig:mult}.}
\label{fig:Rout}
\end{figure}

\begin{figure}
     \centering
     \includegraphics[width=0.8\textwidth]{Rside_fin.eps}
     \caption{\small (Color online)
The $R_{side}$ dependence on transverse momentum for different centralities in the iHKM basic scenario under the same conditions as in Fig. \ref{fig:mult}. The experimental data are from \cite{Alice_radii}.}
\label{fig:Rside}
\end{figure}

\begin{figure}
     \centering
     \includegraphics[width=0.8\textwidth]{Rlong_fin.eps}
     \caption{\small (Color online)
The $R_{long}$ dependence on transverse momentum for different centralities in the iHKM basic scenario under the same conditions as in Fig. \ref{fig:mult}.  The experimental data are from \cite{Alice_radii}.}
\label{fig:Rlong}
\end{figure}

\begin{figure}
     \centering
     \includegraphics[width=0.8\textwidth]{ROS_fin.eps}
     \caption{\small (Color online)
The $R_{out}/R_{side}$ ratio for different centralities in the iHKM basic scenario under the same conditions as in Fig. \ref{fig:mult}.  The experimental data are from \cite{Alice_radii}.}
\label{fig:RoRs}
\end{figure}

In Figs. \ref{fig:Rout}--\ref{fig:Rlong} we present the pion interferometry radii in the iHKM  basic scenario for three different centralities. In general the radii are described well except some small deviation for $R_{side}$ and $R_{out}$ from the central experimental points; these deviations are oppositely directed  for these transverse radii, which results in slight overestimation of the $R_{out}$ to $R_{side}$ ratio, as demonstrated in Fig. \ref{fig:RoRs}.

\section{Summary}
The integrated hydrokinetic model (iHKM) of \emph{A+A} collisions is developed. It includes the generation of the initial (generally momentum anisotropic) state, matter evolution at the pre-thermal stage leading to thermalization, subsequent viscous hydrodynamic expansion, particlization, and hadronic cascade UrQMD. This model is applied to describe at different centralities the mid-rapidity charged-particle multiplicities, the pion, kaon, and antiproton spectra, the $v_2$ coefficients for all charged hadrons, and the pion interferometry radii in Pb-Pb collisons at LHC with $\sqrt{s_{NN}}=2.76$ TeV. It is found that a quite satisfactory description of the bulk observables can be reached when the initial states are attributed to the very small initial
time 0.1 fm/c,  the pre-thermal stage (thermalization process)   lasts at least until  1 fm/c, and the shear viscosity at the hydrodynamic stage of the matter evolution has its minimal value, $\eta/s = \frac{1}{4\pi}$. The requirement of satisfactory description of the soft pion enhancement phenomenon also leads to the further discrimination between the iHKM scenarios, namely, the initial state should be maximally momentum anisotropic and have a small mean relaxation time,  $\tau_{rel}\ll \tau_{th}$. In this basic scenario the bulk observables are described in a wide $p_T$ region, including the pion and antiproton spectra, as well as their yields when  annihilation processes at the afterburner (UrQMD) stage are taken into account \cite{hydrokin-2a}.

 It is worthy noting that while the found tendencies, which direction
the changes of spectra and $v_2$ values go when changing the basic model parameters, are probably stable, an accounting of event-by-event fluctuations, as well as modification of such important factors as EoS and particlization temperature, can change the numerical values of the model parameters, ensuring an optimal description of the observables. Such an investigation is the subject of subsequent work.

Another important point which we would like to stress is the approximate similarity of the results at different reasonable values of the main parameters of the model, even including the pure hydrodynamic scenario starting from $\tau_{th}\rightarrow \tau_0=0.1$ fm/c. The reason for such similarity  is probably that each change of basic parameters is accompanied by a re-normalization of the maximal energy density of the initial state, to keep the agreement with the experimental multiplicity of all charged particles at midrapidity. Then, at the same particlization temperature for all the scenarios, the main characteristics of the bulk observables are approximately preserved. Of course, some details are different, as iHKM demonstrates, but the variations are not so dramatic. This observation explains satisfactory agreement of  different variants of hydrodynamic and hybrid models with the experimental data for A+A collisions, especially when not all bulk observables are considered.

In this paper the Monte Carlo Glauber model is used for the generation of the initial state. The profile from GLISSANDO is associated here with the one for energy density of the non-equilibrium  and non-thermal initial state. Note that the entropy density cannot be used to characterize such a state at $\tau_0=0.1$ fm/c. The interesting point is that with such treatment of the Monte Carlo Glauber model, the share of the binary collisions, contributed to the total energy density, is constant at about 1/4 for LHC energies at {\it different}  iHKM parameters. We are planning to utilize other initial state models, such as EPOS and IP-Glasma, which have other structures of the initial state and already fixed values of the maximal initial energy density, as the generators of the initial states in iHKM. The comparison of the models have to be done on an event-by-event basis and include besides spectra, their anisotropy and femtoscopic scales for different hadron pairs, as well as fluctuations of the observables.

\begin{acknowledgments}
The research was carried out within the scope of the EUREA: European Ultra Relativistic Energies Agreement (European Research Group: ``Heavy ions at ultrarelativistic energies''), and is  supported by the National Academy of Sciences of Ukraine (Agreements F7-2016). I.K. acknowledges financial support by the Helmholtz International Center for FAIR and the Hessian LOEWE initiative.

\end{acknowledgments}