Longitudinal Fluctuations

\usepackage{nath}
\delimgrowth=1
\def\d{{\rm d}}
\def\calS{{\cal S}}
\def\IC{{\rm IC}}
== Introduction ==
Monte Carlo Glauber has been used in order to generate the transverse energy density distribution at the initial time refs.
In this approach the positions of the nucleons are randomly distributed according to the nucleus charge distribution.
With this the event by event analysis can be carried out in the phenomenological side and a finer comparison can be done with the experimental results.
In this paper, we propose and extension of the MC Glauber approach including the longitudinal direction.
This extension is performed in order to consistently preserve energy and momentum during the collision.
With this we generate initial conditions with longitudinal fluctuations.
== Initial condition coordinates ==
In order to convert the energy and momentum information of the nucleons into the hydrodynamic initial condition, one must compute the change of coordinates from the Cartesian coordinates before the collision into the hyperbolic coordinates where the hydrodynamical initial conditions are defined. In this section we calculate the relation between these two quantities in order to distribute the energy and momentum consistently in the hyperbolic hydrodynamic initial condition.
From the energy-momentum tensor T^{\mu\nu} we compute the momentum flux by integrating the energy-momentum tensor on a surface as
:
P^\mu = \int_\calS \d\sigma_\nu T^{\mu\nu}
where \calS is the surface and \d\sigma_\nu is the surface element.
From the conservation of energy and momentum in general coordinates in differential form,
:
T^{\mu\nu}{}_{;\nu} = 0,
we compute the 4-dimensional integral in general coordinates
:
\int \sqrt{-g} \d^4x\ T^{\mu\nu}{}_{;\nu} = 0.
and reexpress that in terms of the direvative and Christoffel symbols
:
\int \sqrt{-g} \d^4x (\partial_\nu T^{\mu\nu} + \Gamma^\mu_{\nu\alpha} T^{\nu\alpha} + \Gamma^\nu_{\nu\alpha} T^{\alpha \mu}) = 0
which can be further simplified into
:
\int \sqrt{-g} d^4x (\frac{1}{\sqrt{-g}}\partial_\nu( \sqrt{-g} T^{\mu\nu} ) + \Gamma^\mu_{\nu\alpha} T^{\nu\alpha} ) = 0
By applying the divergence theorem in the term with the derivative
:
\oint \d^3\sigma_\nu \sqrt{-g} T^{\mu\nu}
= - \int \sqrt{-g} d^4x \Gamma^\mu_{\nu\alpha} T^{\nu\alpha}
we get the integral equation for the energy-momentum conservation.
In the rhs the terms depend on the derivatives of the metric which means that if we choose a coordinate system where the metric is constant, the rhs term completely vanishes.
General expression
:
\d^3\sigma_\nu = \epsilon_{\nu\alpha\beta\gamma}\frac{\partial x^\alpha }{ \partial u } \frac{\partial x^\beta }{ \partial v } \frac{\partial x^\gamma }{\partial w }\d u \d v \d w
=== Cartesian metric ===
The energy and momentum of the collision are defined in Cartesian coordinates. In this case (with constant metrics) the Christoffel symbols vanish above equation reduces to
:
\oint_{\cal S} \d^3\sigma_\nu T^{\mu\nu} = 0.
and now we can see that, in Cartesian coordinates, the momentum flux P^\mu = (P^t, P^x, P^y, P^z) is a conserved quantity.
This means that for any closed surface, the amount of momentum flux inwards is exactly balanced by the amount of momentum flux outwards.
Therefore, we can open the surface in two regions where one is suitably chosen to be the initial condition and the momentum flux on the rest of the surface is constrained to be the same as the one in the initial condition
:
P^{\mu}_{\IC} = \int_{\IC} \d\sigma_\nu T^{\mu\nu} = - \int_{{\cal S} - \IC} \d\sigma_\nu T^{\mu\nu}.
For convenience we first calculate the case where the IC surface is a time constant one. In this academic exercise, the surface elements assume the form of
:
\d\sigma_\nu = \d x\d y \d z\ (1, 0, 0, 0)
and therefore the momentum-flux expressed as
:
P^\mu_{t=t_0} = \int \d x\d y \d z\ T^{\mu t}.
By using ideal hydrodynamics to express the energy momentum tensor we get that the momentum flux density for a time-constant surface is
\frac{\d P^\mu_{t=t_0}}{\d x\d y \d z} = \omega \gamma u^\mu - Pg^{t\mu}
So that the energy flow density is
\frac{\d E_{t=t_0}}{\d x\d y \d z} = \omega \gamma^2 - P = \varepsilon \gamma^2 + P( \gamma^2 - 1 )
and the longitudinal momentum flux is
\frac{\d P^z_{t=t_0}}{\d x\d y \d z} = \omega \gamma^2 v^z
=== \tau-constant IC ===
In the case of heavy-ion collisions, we assume that the initial condition is given by a surface of equally distant points (in proper time) from the collision point. Since the collision is taken to be between two sheets of Lorentz contracted nuclei, only the longitudinal coordinate enters the calculation. Therefore the proper time is defined as
\tau = \sqrt{t^2 + z^2}
and the IC surface is defined on a \tau-constant curve.
The surface elements in this surface are written (in Cartesian coordinates) as
:
\d\sigma_\nu = \tau \d x\d y \d\eta\ (\cosh \eta, 0, 0, -\sinh\eta)
therefore,
:
P^\mu_{\tau=\tau_0} = \int \tau \d x\d y \d\eta\ ( \cosh \eta T^{\mu t} - \sinh \eta T^{\mu z} )
By using the ideal hydrodynamics energy-momentum tensor
:
T^{\mu\nu} = \omega u^\mu u^\nu - g^{\mu\nu} P,
we can compute the total energy and momentum of the surface as
:
P^\mu_{\tau=\tau_0} 
= \int \tau \d x\d y \d\eta( (\omega \gamma u^\mu - g^{\mu t} P)\cosh \eta - (\omega u^z u^\mu - g^{\mu z} P) \sinh \eta )
:
= \int \tau \d x\d y \d\eta ( \omega \gamma u^\mu \cosh\eta ( 1 - v^z \tanh\eta ) - P\cosh\eta (g^{\mu t} - g^{\mu z} \tanh \eta ) )
So that the (conserved) Cartesian energy and momentum densities can be expressed on a \tau-constant IC
:
\frac{\d P^\mu_{\tau=\tau_0}}{\tau \d x\d y \d \eta}
= \cosh\eta ( \omega \gamma u^\mu ( 1 - v^z \tanh\eta ) - P (g^{\mu t} - g^{\mu z} \tanh \eta ))
Or explicitly as
:
\frac{\d P^\mu_{\tau=\tau_0}}{\tau \d x\d y \d \eta}
= \cosh\eta ( \omega \gamma^2 ( 1 - v^z \tanh\eta ) - P )
\\
= \cosh\eta ( \omega \gamma^2 u^z ( 1 - v^z \tanh\eta ) - P \tanh \eta )
The Cartesian energy and momentum densities are expressed in terms of the local invariant energy density, the fluid velocity, the pressure and the spatial rapidity coordinate all calculated on a τ-constant surface.
From the MC Glauber one can calculate what is the total energy and momentum involved in the collision. 
This puts a constraint on the integrals, however the actual transverse and longitudinal distributions of the local energy density and the fluid velocity are unconstrained. 
In fact, the determination of both the total energy and the total momentum from the collision do impose relations between these distributions so that same amount of energy and momentum are present in the IC as what is missing after the nuclei collision.
As we will see below, the inclusion of the total momentum constraint allows us to consistently distribute the local energy density in the longitudinal direction.
=== Velocity profile ===
An usual assumption when simulating heavy ion collisions is that the initial longitudinal flow is given by
:
v^z = \frac{z}{t}.
This means that the fluid has propagated in free-streaming lines until it reached the \tau-constant surface.
Using the hyperbolic coordinates,
:
\tau = \sqrt{ t^2 + z^2}
\\
\eta = \tanh^{-1}\frac{z}{t}
One can reexpress the longitudinal velocity as
:
v^z = \tanh\eta
 
The relation between the Cartesian velocities and the hyperbolic ones can be derived directly from the definition as
:
u^t = \cosh\eta\ u^\tau + \tau \sinh\eta\ u^\eta
\\
u^z = \sinh\eta\ u^\tau + \tau \cosh\eta\ u^\eta
With this we can express the general form of the conserved energy-momentum density in terms of the hyperbolic velocities as
:
\d E = T^{tt} \cosh\eta - T^{tz} \sinh\eta = (\omega u^t{}^2 - P)\cosh\eta - \omega u^t u^z \sinh\eta = (\omega u^\tau{}^2 - P)\cosh\eta + \tau \omega u^t u^\eta \sinh\eta = T^{\tau\tau} \cosh\eta + \tau T^{\tau\eta} \sinh\eta
\\
\d P_x = T^{xt} \cosh\eta - T^{xz} \sinh\eta = \omega u^t u^x\cosh\eta - \omega u^x u^z \sinh\eta = \omega u^\tau u^x = T^{\tau x}
\\
\d P_z = T^{zt} \cosh\eta - T^{zz} \sinh\eta = \omega u^t u^z \cosh\eta - (\omega u^z{}^2 - P) \sinh\eta = \omega u^\tau u^\eta \cosh\eta + \tau (\omega u^\tau{}^2 - P) \sinh\eta = \tau T^{\tau\eta} \cosh\eta + T^{\tau\tau} \sinh\eta
Therefore the above imposition translates into
:
u^\eta = 0.
Using this assumption we revisit the energy and momentum densities expressions
:
\frac{\d P^\mu_{\tau=\tau_0}}{\tau \d x\d y \d \eta}
= \cosh\eta ( \omega \gamma^2 ( 1 - (v^z)^2 ) v^\mu - P (g^{\mu t} - v^z g^{\mu z} ))
where we substituted the hyperbolic tangent terms by the longitudinal velocity. 
With this we observe that the term
:
\gamma^2 ( 1 − ( v^z )^2)
vanishes if one further assumes that v^x = v^y = 0 (so that \gamma^2 = 1 + (u^z)^2 = \frac{1}{1 - (v^z)^2}).
This means that there also no initial transverse flow.
By applying this second approximation we further simplify the relation between the Cartesian energy and momentum densities and the hydrodynamic distributions
:
\frac{\d P^\mu_{\tau=\tau_0}}{\tau \d x\d y \d \eta}
= \cosh\eta\ ( \omega v^\mu - P (g^{\mu t} - g^{\mu z} \tanh \eta ))
\\
= \varepsilon \delta^{\mu z} \sinh \eta + \varepsilon \delta^{\mu t} \cosh\eta
The above expression depends only on the local energy density profile and the \eta coordinate.
Explicitly as
:
E = \int \d x\d y \tau\d \eta\ \varepsilon(x,y,\eta) \cosh\eta
\\
P_z = \int \d x\d y \tau\d \eta\ \varepsilon(x,y,\eta) \sinh\eta
In this section we have applied the assumption that the velocity profile is zero in the x, y and \eta directions in the IC.
With this we simplify the relation between the conserved energy and momentum densities and the local densities.
In particular, the dependence of the pressure drops and the conserved quantities depend only on the local energy density and the \eta coordinate.
With this assumption, the additive quantities of the conserved energy and momentum are completely determined by integrals over the additive quantity of the local energy density.
This allows us to apply the same partitioning on both the lhs and rhs of the equation above as
:
E = \sum_i E_i = \int \d x\d y \tau\d \eta\ \cosh\eta\ \sum_i \varepsilon_i(x,y,\eta) 
\\
P_z = \sum_i P^z_i = \int \d x\d y \tau\d \eta\ \sinh\eta\ \sum_i \varepsilon_i(x,y,\eta) 
which is very convenient in order to convert the energy and momentum of the MC Glauber picture into the hydrodynamic IC as we shall see below.
Furthermore, the longitudinal profile of the local energy density is constraint by both the conserved energy and the conserved momentum.
In the next section we shall discuss an Ansatz for the local energy distribution that satisfies these constraints.
If one ignores the conserved momentum (as often done in the literature), the longitudinal profile of the local energy density is not constrained.
=== Energy density profile ===
It is clear from the above expressions that in the case that \varepsilon is an even function of \eta, there is no longitudinal momentum.
Therefore, in order to generate an initial condition with the appropriate conserved energy and momentum, \varepsilon must be such that it satisfies the above constraints.
Since the difference in the computation between the energy and momentum is due to the \eta-dependence, we factorize the energy density profile into \varepsilon(x,y,\eta) = \varepsilon_\perp(x,y) N(\eta_0)\varepsilon_\eta(\eta - \eta_0 ) so that
:
E = I_\perp N(\eta_0) \int \tau\d \eta\ \varepsilon_\eta(\eta - \eta_0) \cosh\eta
\\
P_z = I_\perp N(\eta_0) \int \tau\d \eta\ \varepsilon_\eta(\eta - \eta_0) \sinh\eta
where
:
I_\perp = \int \d x\d y\ \varepsilon_\perp(x,y)
is the integral of the transverse part of \varepsilon.
Without loss of generality we assume 
:
I_\perp = 1
which means that \varepsilon_\perp(x,y) is a normalized distribution with finite support (generaly a Gaussian).
In the above expressions, \eta_0 are longitudinal shifts that account for the non-zero momentum and N(\eta_0) are the normalizations that ensure that the conserved energy is not affected by the shifts.
We can simplify the above expressions by displacing the integral so that shifts appear in the hyperbolic functions.
After that we use the relations of the hyperbolic functions to write
:
E = N(\eta_0) (\cosh\eta_0\ \int \tau\d \eta\ \varepsilon_\eta(\eta) \cosh\eta + \sinh\eta_0\ \int \tau\d \eta\ \varepsilon_\eta(\eta) \sinh\eta )
\\
P_z = N(\eta_0) (\cosh\eta_0\ \int \tau\d \eta\ \varepsilon_\eta(\eta) \sinh\eta + \sinh\eta_0\ \int \tau\d \eta\ \varepsilon_\eta(\eta) \cosh\eta )
where the shifts are out of the integrals.
The remaining integrals are over the longitudinal part of \varepsilon_\eta centered at \eta_0=0.
As we discussed above, this is the case where there is no longitudinal momentum in the system, and therefore, because of symmetry, \varepsilon_\eta is taken to be an even function in \eta.
This choice means that the integrals containing the \sinh functions vanish and the above expressions reduce to
:
E = N(\eta_0) \cosh\eta_0\ I_\eta
\\
P_z = N(\eta_0) \sinh\eta_0\ I_\eta
where
:
I_\eta = \int \tau\d \eta\ \varepsilon_\eta(\eta) \cosh\eta
which is independent of the shifts and can be computed once the actual shape of the distribution is given.
In order to determine the shifts, one can divide both expressions
:
\tanh \eta_0 = \frac{P_z}{E}
which is consistent with our previous Ansatz that in the case of no shifts, no momentum is generated.
The normalization is easily determined by the conserved energy
:
N(\eta_0) = \frac{E}{I_\eta \cosh \eta_0 }
which means that more longitudinal momentum translates into larger shifts and therefore smaller densities.
This is due to the already present kinetic energy given by the velocity profile.
Once the local energy density is shifted, it probes larger velocities and contribute with more kinetic energy so that the local rest energy density (in eta coordinates) drops to compensate.
The final for of the local energy density is given by
:
\varepsilon(x,y,\eta) = \frac{E}{I_\eta \cosh \tanh^{-1}(\frac{P_z}{E}) }\varepsilon_\perp(x,y) \varepsilon_\eta(\eta - \tanh^{-1}(\frac{P_z}{E}) )
where the \varepsilon_\perp and \varepsilon_\eta are taken to be normalized finite support distributions (generally 2d and 1d Gaussians, respectively).
The normalized transverse energy distribution is given by
:
\frac{1}{E}\frac{\d E}{\cosh\eta\ \tau\d\eta} = \frac{1}{E}\frac{\d P_z}{\sinh\eta\ \tau\d\eta} = \frac{1}{E}\int\d x\d y\ \varepsilon(x,y,\eta) = \frac{1}{I_\eta \cosh \tanh^{-1}(\frac{P_z}{E}) }\varepsilon_\eta(\eta - \tanh^{-1}(\frac{P_z}{E}) )
This quantity is interesting to compare different distribution choices as we shall see below.
=== Plots ===
We plot the normalized transverse energy distribution below for a choice of a Gaussian profile for the longitudinal distribution.
The only free parameter is the Gaussian width \sigma_\eta (which determines I_\eta.
In the plots Eeta denotes
:
\frac{1}{E}\frac{\d E}{\cosh\eta_0 \tau \d\eta}(\eta; \sigma_\eta; \frac{P_z}{E})
The first plot shows the conserved momentum distribution for different fractions of momentum and different Gaussian widths.
This clearly show the influence of the velocity profile.
A vanishing total momentum means that the distribution is centered at \eta=0 so that both sides perfectly cancel.
The second plot shows the conserved energy distribution for different fractions of momentum and different Gaussian widths.
The two peaks are also a consequence of the choice of the velocity profile and the hyperbolic nature of the initial condition surface.
The peaks show the maximum kinetic energy values.
The third plot shows the local energy density.
As discussed above, larger shifts suppress the amplitude of the distribution.
We have proposed an Ansatz for the local energy density profile that together with the flow scaling Ansatz completely determined the hydrodynamic initial condition in a \tau-constant surface consistent with the conserved energy and momentum given.
In addition, we can perform the partitioning of the energy and momentum and connect the MC Glauber picture with the hydrodynamic IC consistently.
In the next section we shall discuss how to partition the energy and momentum in the MC Glauber picture.
== Initial condition generator ==
In the MC Glauber picture, the 3-dimensional position of the nucleons is distributed according to the charge distribution (Woods-Saxon)
:
n( r ) = \frac{n_0}{ \exp( \frac{r - R}{\xi} ) + 1 }.
Each nucleon i carries an equal ammount of energy and momentum of the total nucleus
:
p^z_i = \pm \sqrt{s_{\rm NN}}/2
\\
e_i = \sqrt{ m_p^2 + p_z^2 }
where the sign of the momentum is different for both colliding nuclei.
With this we can compute the rapidity of these nucleons as
:
Y_i = \tanh^{-1} \frac{p_z}{\sqrt{ m_p^2 + p_z^2 }} = \tanh^{-1} \frac{1}{\sqrt{ 4\frac{m_p^2}{s_{\rm NN}} + 1 }}
In the case of the LHC (\sqrt{s_{\rm NN}} = 2760GeV) we get Y_i \sim 8
Next we chose an impact parameter b and shift the nucleons in each nuclei for either b/2 or -b/2 in the x-direction.
Because of the Lorentz boost, the two nuclei are represented by thin sheets and amplitude of the collision between the nucleons i and j is computed by the condition
:
\Gamma( d_{ij} ) = \Gamma( \sqrt{ \frac{(x_i - x_j)^2 + (y_i - y_j)^2}{\sigma_{\rm NN}} } )
In the hard sphere approximation this amplitude takes the form of
:
\Gamma( d_{ij} ) = 1\ {\rm if}\ d_{ij}<R\ {\rm else}\ 0
where R is determined by the constraint that
:
\int \d^2 x \Gamma( x ) = 1
therefore
:
\int \d^2 x \Gamma( x ) = \pi R^2
\\
R = \frac{1}{\sqrt{\pi}}
This means that if two nucleons from different nuclei overlap in the transverse plane, they collide with each other.
From this two straightforward quantities can be defined, the number of colliding nucleons (participants, N_{\rm part}) and the number of pairs of colliding nucleons (binary collisions, N_{\rm bin}).
In order to distribute the energy of the incoming nuclei into the hydrodynamic initial condition, it is usual to assume that each collision generates a transverse energy distribution (generally a Gaussian) which takes a fraction of the total energy of the system.
The total energy is estimated by comparing the final results such as the total multiplicity or the total energy measured.
In this simple approach, the momentum is neglected.
Furthermore, there is still an ambiguity in defining the collisions.
They could be defined as the CM of two nucleons involved in a binary collisions, meaning there would be N_{\rm bin} independent energy density distributions which each carries 1/N_{\rm bin} of the total energy.
Another scenario is to place the energy distributions in the center of wounded nucleons, meaning there would be N_{\rm part} independent energy density distributions which each carries 1/N_{\rm part} of the total energy.
Both pictures assume that the colliding nucleons deposit a fixed amount of energy in the collision region and no momentum.
This means that in the participant picture, each wounded nucleon will have lost the same energy regardless of how many collisions it took part.
Also as no momentum is lost by the nucleon, this pictures generates inconsistencies since the momentum of the nucleon might be larger than its energy after the collision.
In the binary picture, there is also the problem that since a single nucleon can undergo multiple collisions, by subtracting a fixed amount of energy from it on each collision, one can end up with a nucleon with negative energy after the collision process.
In order to address these issues, we propose that first both energy and momentum are subtracted from each nucleon during collision.
Second, instead of subtracting a fixed amount, a fraction \alpha of the remaining energy and momentum of the nucleons is taken for each collision it is involved.
In the next sections we explore how this applies to the binary and participant scenarios.
=== Binary scenario ===
The assumptions addressed above can be straightforwardly applied to the binary scenario.
First we count the number of collisions each nucleon participates, N.
With this we can compute its momentum after the collision to be
:
p_z' = (1 - \alpha)^N p_z
The same can be done to the energy, however in this work, for simplicity, we assume that all nucleons are protons and they keep their masses after the collision so that e' = \sqrt{ m_p^2 + p_z'^2 }.
By subtracting a fraction \alpha for each collision we guarantee that the final energy is always non-negative and also by keeping the mass we guarantee that the final energy is always larger than the momentum.
This can be seen by computing the case of a nucleon participating on a very large number of collisions, i.e., N\gg 1.
In this case the ratio between the final quantities and the initial ones are
:
\frac{p_z'}{p_z} \to \exp(-N\alpha) 
\\
\frac{e'}{e} \to \frac{ m_p }{ e } + \exp(-N\alpha) \frac{ p_z^2 }{ m_p e }
On the other hand, in the limit of vanishing \alpha, i.e., \alpha\to 0 we get
:
\frac{p_z'}{p_z} \to 1 - N\alpha,
\\
\frac{e'}{e} = 1 - N\alpha\frac{p_z^2}{2e^2}.
Now that we have computed the total energy and momentum loss for each nucleon, we distribute this energy into the collisions they have each participated.
Because of Lorentz contraction, we assume that all these collisions happened simultaneously and therefore, each collision receives an equal amount of the total energy lost by each nucleon.
In other words, if nucleon A participates in N^A collisions, each of the collision will receive
:
\frac{p_z^A - p_z'^A}{N^A} = \frac{1 - (1 - \alpha)^{N^A} }{N^A} p_z = \Delta_p(N^A)
\\
\frac{e^A - e'^A}{N^A} = p_z \frac{ \sqrt{ \frac{m_p^2}{p_z^2} + 1 } - \sqrt{ \frac{m_p^2}{p_z^2} + (1 - \alpha)^{N^A} } }{N^A} = \Delta_e(N^A) 
here we defined the quantities \Delta_p and \Delta_e as the momentum and energy deposition per collision (which depend on the number of collisions), respectively.
In the binary scenario, we assume that for each binary collision a tube is generated which gets its energy and momentum from both the nucleons involved in that binary collision.
The tube is transversely position in the center of mass of the collision, r' = \frac{r_A + r_B}{2}.
Therefore, the final energy and momentum of the tube generated by the binary collision between nucleon A (involved in N^A collisions) and B (involved in N^B collisions) is
:
p_z^{AB} = \Delta_p(N^A) - \Delta_p(N^B) = (\frac{1 - (1 - \alpha)^{N^A}}{N^A} - \frac{1 - (1 - \alpha)^{N^B}}{N^B}) \frac{\sqrt{s_{\rm NN}}}{2}
\\
e^{AB} 
= \Delta_e(N^A) + \Delta_e(N^B)
We compute the energy and momentum consistently for each tube and then we can compute the its \eta shifts and fully determine the local energy density distribution.
The total hydrodynamic initial condition is then the superposition of all the tubes after the collision.
From the above definition, we note that it is clear that if both the colliding nucleons take part in the same number of collisions the momentum of that single collision is zero.
On the other hand, if they are part of a different number of collisions, the tube gets a residual momentum.
This momentum is maximum when one nucleon takes part only in that collision N^A=1 and the other takes part in many other collisions N^B\gg 1.
This situation happens for instance in pA collisions, where the proton experiences multiple collisions while each nucleon experiences only one.
In this case we get
:
p_z^{AB} = \Delta_p(1) - \Delta_p(N^B) \to \alpha p_z
which is the maximum momentum a tube can have in the binary scenario.
Below we make an example of a DD collision where the one incident proton D_1 hits the two target protons (D'_1, D'_2) while the other (D_2) does not hit (spectator).
In this example, there are two tubes formed since there are two binary collisions.
Schematically
:
D+D = (D_1,D_2) + (D'_1,D'_2) = (D_1 + D'_1), (D_1 + D'_2), D_2 \to (D^{**}_1 + D'^{*}_1 + X), (D^{**}_1 + D'^{*}_2 + X), D_2
where the * denotes the excitation from each collision and X denotes the tubes.
The first proton D_1 participates in both (N=2) while the target ones participate each in one.
The momenta of the generated tubes (from the binary collisions D_1+D'_1 and D_1+D'_2) are
:
p_z^{D_1D'_1} = \Delta_p(2) - \Delta_p(1) = p_z ( \frac{1}{2}(1 - \alpha)^2 - (1 - \alpha) ) = - \frac{p_z}{2} ( 1 - \alpha^2)
\\
p_z^{D_1D'_2} = \Delta_p(2) - \Delta_p(1) = p_z^{pD_1}
And the total longitudinal momentum of the hydrodynamic IC is non zero
:
P_z = 2\Delta_p(2) - 2\Delta_p(1) = - p_z ( 1 - \alpha^2)
The figure below shows the rapidity configuration after the collision.
The participating nucleons (represented by disks) suffer a loss of rapidity and the tubes (represented by ellipses) have are shifted in rapidity.
In the next figures we compare this binary scenario for full PbPb collisions, central and peripheral.
We note that although the participants experience different numbers of collisions (by their rapidity losses) the tubes are mostly centered at zero rapidity.
Also there is a slight tilt of the IC for the peripheral case where the IC tends to tilt as if it was dragged by the nuclei.
Another feature from this model is the very large number of tubes which scale with the number of binary collisions.
=== Participants scenario ===
In the participant scenario, the tubes are generated by each wounded nucleon (participating).
Therefore, one could follow the same procedure as before but instead of adding the two contributions from the colliding pair of nucleons, deposit only that wounded nucleon energy and momentum in each tube.
This brings up the problem that the momentum loss is can be larger than the energy loss.
In fact, for \alpha \to 0
:
\frac{p_z - p'_z}{e - e'} = \frac{\Delta_p(N)}{\Delta_e(N)} \to \frac{e}{p_z} + \frac{m_p^2 N \alpha}{2p_z e}
Therefore, the simple picture where each tube receives only the energy and momentum from the wounded nucleon fails.
To amend that, we must include the contributions of nucleons that were hit by that wounded nucleon.
Therefore, the tube receives the contribution from the incident nucleon itself and the other (N) nucleons in the target nucleus that it hits.
On the other hand, in the perspective of the target the collision changes and each of the target nucleons hits the incident nucleon plus others.
To accommodate for the symmetry of the collision, we sum the two views so that the participant nucleon itself being part of N+1 collisions, where N is the number of binary collisions (as target) and its collision (as projectile).
The momentum and energy loss are completely analogous to the binary case, only incrementing one to the number of binary collisions.
:
p_z' = (1 - \alpha)^{N+1} p_z
And the total energy and momentum loss per collision is completely analogous to the previous case
:
\frac{p_z - p_z'}{N + 1} = \Delta_p(N+1)
\\
\frac{e - e'}{N+1} = \Delta_e(N+1)
The difference comes on how to distribute this lost energy and momentum.
Each wounded nucleon generates a tube in the center of mass of its collision against N nucleons of the other nucleon r' = \frac{ r_A + \sum_n r_n}{ N + 1 } and its momentum content is calculate from the contribution of all the nucleons involved in that collision.
Each tube then takes the contribution of the incident nucleon and the multiple (N^A) targets it hit
:
p_z^{AB_1B_2...} = \Delta_p(N^A+1) - \sum_n^{N^A} \Delta_p(N^{B_n}+1) 
In this case, the maximum amount of momentum deposited in the tube can be much larger when compared to the binary scenario.
In the case of pA with the proton hitting a very large number N^A\gg 1 of targets (N^{B_i}=1) the momentum is
:
p_z^{pB_1B_2...} = \Delta_p(N^A+1) - N^A \Delta_p(2) \to - N^A \frac{1-(1-\alpha)^2}{2} p_z
We compare the previous example of a DD collision with this scenario.
Schematically it gets
:
p + D = p + (D'_1,D'_2) = p + (D'_1,D'_2) \to p^{***} + (D'^{**}_1,D'^{**}_2) + X + 2 X'
where the X stands for the tubes from the proton and X' are the tubes from the deuteron.
There are 3 tubes since there are 3 participants.
The momentum distribution in further detail is
In the case of pD where all nucleons hit we have three tubes
:
p_z^{pD_1 D_2} = \Delta_p(3) - \Delta_p(2) - \Delta_p(2)
\\
=((\frac{1 - (1 - \alpha)^{3}}{3} - 2(\frac{1 - (1 - \alpha)^{2}}{2}) p_z
\\
p_z^{D_1 p} = \Delta_p(3) - \Delta_p(2)
\\
=((\frac{1 - (1 - \alpha)^{3}}{3} - (\frac{1 - (1 - \alpha)^{2}}{2}) p_z
\\
p_z^{D_2 p} = \Delta_p(3) - \Delta_p(2) = p_z^{D_1 p}
and the total momentum is
:
P_z = 3\Delta_p(3) - 4\Delta_p(2)
this can be seen in the picture below where we perform a DD collision where one nucleon is not participating.
In the next figures we compare the participant scenario for full PbPb collisions, central and peripheral. 
We note that the tubes are much more shifted from the center in this scenario compared to the binary case.
Due to the multiple nature of the collisions, the tubes are always shifted to the opposite direction of the participant momentum (or not shifted at all).
This generated two tilted strings of tubes are observed, once coming from the participants of each nucleus.
Also there is a slight tilt of the IC for the peripheral case where the IC tends to tilt as if it was dragged by the nuclei. 
The number of tubes is also reduced in this scenario when compared to the binary case.