Corresponding contributions of core and corona components in $0$-$10\%$ multiplicity and centrality classes in (c) $p$+$p$ and in (d) $Pb$+$Pb$ collisions, respectively.
Corresponding contribution of core and corona components in $50$-$100\%$ multiplicity and centrality classes in (e) $p$+$p$ and in (f) $Pb$+$Pb$ collisions, respectively.

Figure <ref> (a) shows the mean transverse mass, $\langle m_T \rangle -m$, as a function of the rest mass of hadrons, $m$, in high-multiplicity ($0$-$10\%$) and low-multiplicity ($50$-$100\%$) $p$+$p$ collisions at $\sqrt{s}= 7$ TeV.
To pin down the effect of hadronic rescatterings, the results from full simulations and simulations without hadronic rescatterings are compared with each other.
The result from 8 with the default settings including color reconnection is also plotted as a reference.
Overall, $\langle m_T \rangle-m$ in high-multiplicity events ($0$-$10\%$) tends to exhibit an almost linear increase with increasing $m$ except for phi mesons.
On the other hand, such a clear mass dependence is not seen in low-multiplicity events ($50$-$100\%$), which is consistent with the $m_T$ scaling.
It should also be mentioned that the violation of linear increase of phi mesons in core components appears after resonance decays against direct hadrons (not shown).
An apparent flow-like linear mass dependence is seen in results from 8 in high-multiplicity events as well, which is due to the color reconnection <cit.>.
The almost linear increasing behavior in DCCI2 is caused by both radial flow for core components from hydrodynamic expansion and color reconnection for corona components obtained with 8.
As a result, the results from both DCCI2 and 8 have similar tendencies.
Therefore, it is difficult to discriminate each effect by merely seeing the mean transverse mass.
The effect of hadronic rescatterings is almost absent for pions. This comes from an interplay between small $pdV$ work in the late hadronic rescattering stage and approximate conservation of pion number <cit.>.
The small effects of hadronic rescatterings are seen for phi mesons and omega baryons because they do not form resonances in scattering with pions unlike other hadrons <cit.>.

Figure <ref> (b) shows the mean transverse mass as a function of hadron rest mass for $0$-$10\%$ and $50$-$100\%$ centrality classes in $Pb$+$Pb$ collisions at = $2.76$ TeV.
The almost linear increasing behavior of $\langle m_T \rangle-m$ appears even in $50$-$100\%$ centrality class as one can expect from the centrality dependence of the fraction of core components in Fig. <ref> (b).
The larger enhancement of the mean transverse mass due to hadronic rescatterings, in particular, for protons is seen in high-multiplicity events in comparison with the low-multiplicity events.
This is a manifestation of the famous “pion wind" in the late rescattering stage <cit.>.

Figure <ref> (c) shows each contribution of core and corona components to the final result without hadronic rescatterings in 0-10% multiplicity class in $p$+$p$ collisions. The inclusive result here is identical to the one shown as the result without hadronic rescatterings in Fig. <ref> (a).
The difference between results of core and corona components is seen in protons, lambdas, and omega baryons.
The linear mass ordering of $\langle m_T \rangle-m$ from core components is slightly diluted
for protons and lambdas in the inclusive result
due to the sizable contribution of corona components.
In contrast, the core result and the inclusive result are almost on top of each other since the contribution of corona components for omega baryons is smaller in 0-10% multiplicity class compared to other particle species.

Figure <ref> (d) shows each contribution of core and corona components to the final result without hadronic rescatterings in 0-10% centrality class in $Pb$+$Pb$ collisions. The linear increase except phi meson is seen very clearly for the core component, and the increase rate is more than the one from $p$+$p$ results shown in Fig. <ref> (c).

Figure <ref> (e) shows the same variable with Fig. <ref> (c) but in 50-100% multiplicity class in $p$+$p$ collisions.
Since the fraction of the core components is less than 10% in this range of multiplicity class as shown in Fig. <ref> (a), the final result and the result from corona components are almost top of each other showing no significant dependence on hadron rest mass.

Figure <ref> (f) shows the same variables with Fig. <ref> (d) but in 50-100% centrality class in $Pb$+$Pb$ collisions. According to Fig. <ref> (b), the fraction of core components shows $R_{\mathrm{core}}\approx 0.8$ to $\approx 0.9$ in this centrality range. Eventually the result of core components is found to be slightly diluted by corona components.

§ ANISOTROPIC FLOW

Anisotropic flow is one of the observables that has been often investigated in high energy nuclear collisions to quantify collectivity.
As I mentioned in Sec. <ref>,
the word “collectivity” means that there is a correlation in azimuthal angle of momentum space for produced particles.
The physical origin of such a correlation has been, so far, understood as follows: in $Pb$+$Pb$ collisions, once there is a geometrical anisotropy in initial state of QGP fluids, the anisotropy in geometry is converted into that in momentum space due to hydrodynamic response. Final hadrons from the QGP fluids reflect the anisotropy even after hadronization, thus the momentum anisotropy in final hadrons is a signal of QGP formation.
On the other hand, in $p$+$p$ collisions, it has been understood that most of particles originate from fragmentation from hard/non-equilibrium particles.
Those hard particles called jets usually are produced in back-to-back due to momentum conservation
and have a strong correlation in momentum space.
Thus, it comes to ambiguous which physical origin momentum anisotropy implies in $p$+$p$ collisions.
In this sense, implication of anisotropy in $p$+$p$ collisions and $Pb$+$Pb$ collisions can be different. Therefore I explain results from $p$+$p$ and $Pb$+$Pb$ collisions separately in the following sequential subsections.

It would be worth mentioning the motivation to see anisotropic flows in DCCI2 here.
As I explained, there are mainly the two factors that pronounce
anisotropy in momentum space: hydrodynamic response of QGP fluids and jets from hard scatterings, and the anisotropy originates from these are often referred to as flow and non-flow, respectively
It should be noted that these are sorts of jargon in relativistic heavy-ion community, and there is not, at my best knowledge, solid definitions.
Generally speaking, it is assumed that flows are manifested in small to intermediate transverse momentum regime and have azimuthal correlations in long range towards (pseudo)rapidity directions.
On the other hand, non-flows are manifested in large transverse momentum regime and have azimuthal correlation in short range correlations.
If one wants to investigate the flow in experimental data,
techniques called non-flow subtractions are adopted to, namely, subtract non-flow correlations,
and those techniques are built based on the aforementioned assumptions.
For instance, imposing a certain gap in long range direction between two-particle pair when one calculate two-particle correlation is expected to remove short range correlations of non-flows.
It is often simply expected that in core–corona picture,
the anisotropy from core and corona components corresponds to flow and non-flow.
In the following subsections,
I also investigate how each anisotropic flow changes when one imposes non-flow subtraction within DCCI2.

§.§ Multiplicity dependence in $p$+$p$ collisions

Second and third order of anisotropic flow coefficients obtained from two-particle correlation for charged hadrons without eta gap as a function of the number of produced charged particles in $p$+$p$ collisions at $\sqrt{s} = 13 \ \mathrm{TeV}$ are shown in left and right figures. Results from full simulations (orange squares) and simulations without hadronic rescatterings (yellow diamonds) are shown.  Results only from core (red diamonds) and corona (blue triangles) components as a breakdown of results of simulations without hadronic rescatterings are shown simultaneously.

Let me start from discussion on results from $p$+$p$ collisions.
Figure <ref> shows the second- and third- order anisotropic flow coefficient of charged particles obtained from two-particle cumulants, $\vtwtw$ and $\vthree$, as functions of $N_{\mathrm{ch}}$ in $p$+$p$ collisions at $\sqrt{s}=13$ TeV.
Kinematic cuts for all of $\vtwtw$, $\vthree$, and $N_{\mathrm{ch}}$ are $0.2<p_T<3.0$ GeV and $|\eta|<0.8$ as used in Ref. <cit.>.
One sees that $\vtwtw$ obtained from both core and corona components is larger than that from simulations without hadronic rescatterings for overall $N_{\mathrm{ch}}$.
This suggests that the event plane angle of core components might be different from that of corona components, which dilutes $\vtwtw$ of core and corona components with each other.
In contrast, $\vthree$ shows different tendency: the results from the summation of core and corona give the values of $\vthree$ lying in-between the results from core and corona components.
It might be possible to understand this tendency that the third-order of anisotropic flow is clearly manifested in corona components
One of the possible origins of the third-order anisotropic flow in corona component is a
Bremsstrahlung from a back-to-back jets.
while those in core components vaguely appears and weaken those seen in the corona components.

Comparisons among second- (green crosses), third- (navy x-s), and fourth- (blue stars) order of anisotropic flow coefficients obtained from two-particle correlation for charged hadrons as functions of the number of produced charged hadrons in $p$+$p$ collisions at [proton] = 13 TeV from full simulations.

Figure <ref> shows
comparisons among second-, third-, and fourth- order of anisotropic flow coefficients obtained from two-particle correlation as functions of the number of produced charged hadrons in $p$+$p$ collisions at [proton] = 13 TeV from full simulations.
One sees that the absolute values of anisotropic flow coefficients decrease in ascending order of the order of coefficients.
It should be noted that, within DCCI2, initial momentum anisotropy is embedded in
initial matter profile due to dynamically deposited energy and momentum from initial partons.
Hence, the origin of anisotropy pronounced even in core component is not limited to geometrical effects.
To investigate the origin of anisotropic flows,
it has been proposed that a sign change of correlation of anisotropic flow and mean $p_T$ in event-by-event is a proxy of the transition of origin of anisotropic flow in final hadrons between
initial momentum anisotropy or geometry of initial condition.
Further investigation can be performed through such observables in the future.

Comparison of with and without eta gap ($\Delta \eta > 1.0$) in second-order of anisotropic flow coefficient obtained from two-particle correlation for charged hadrons as a function of the number of produced charged particles in $p$+$p$ collisions at $\sqrt{s} = 13 \ \mathrm{TeV}$ are shown in left and right figures, respectively. Results from simulations without hadronic rescatterings (yellow diamonds) and its breakdown into core (red diamonds) and corona (blue triangles) components are shown.

Figures. <ref> (left) and (right) show comparisons of with and without eta gap in second-order of anisotropic flow coefficients obtained from two-particle correlations as a function of the number of produced charged particles in $p$+$p$ collisions at [proton] = 13 TeV, respectively.
In order to see the effects of imposition of eta gap with abundant statistics,
correlations are obtained with particles produced within $|\eta|<3.0$ instead of $|\eta|<0.8$.
When eta gap is imposed, each one particle used to calculate 2-particle correlation is taken from 2-sub region independently, and eta gap is imposed in-between the 2-sub regions. Detailed explanation are in Appendix. <ref>,
As I aforementioned at the beginning of this section,
a certain eta gap is expected to remove short range correlation of non-flows
where non-flows are often expected to be the correlations in the corona.
In contrast to the expectation,
one sees that both $\vtwtw$ obtained from core and corona are suppressed
simultaneously due to imposition of eta gap of $\Delta \eta > 1.0$.
Thus, within DCCI2, such non-flow subtraction also subtract correlation from core components.
This is because QGP fluids are dynamically generated via four-momentum deposition of initial partons and there is a certain momentum correlation between core and corona components,
which I discuss more details in Sec. <ref> by seeing two-particle correlations in 2-D map of $\Delta \phi$ and $\Delta \eta$.

Second order of cumulant obtained from four-particle correlation for charged hadrons as a function of the number of produced charged hadrons in $p$+$p$ collisions at [proton] = 13 TeV.
Kinematic ranges are the same as used in (left) the ALICE experiment in which 3-sub event method is imposed and (right) the CMS experiment.

I close discussion on $p$+$p$ collisions with a proxy of QGP formation in small systems nu showing second order of cumulant obtained from four-particle correlation, $\ctwofour$.
It is said that $\ctwofour$ is a better observable rather than $c_2\{2\} = (\vtwtw)^2$
because the former can truncate back-to-back 2-particle correlations, which are assumed to be non-flow and inherit in the latter.
Here, it should be noted that negative $\ctwofour$ is the proxy of
hydrodynamic behavior because of its definition in Eq. <ref>.
Figure <ref> shows $\ctwofour$ as a function of the number of produced charged hadrons in $p$+$p$ collisions.
Left and right figures are results obtained by adopting the same kinematic ranges and methods
used in the ALICE experiment and the CMS experiment in which
negative $\ctwofour$ is observed.
The results shown in the left panel using the ALICE kinematic range are ones that have been adopted
in calculations of anisotropic flows in this section so far.
On the other hand, $\ctwofour$ is calculated with, so called, 3-sub event method where
four particles are picked up from 3-sub rapidity ranges, $-0.8<\eta<-0.4$, $-0.4<\eta<0.4$, and $0.4<\eta<0.8$, to reduce non-flows by cutting short range correlations. For details see Appendix. <ref>.
The CMS kinematic range adopted to calibrate the results shown in the right panel is as follows:
correlations are obtained by using particles produced with $0.3<p_T<3.0$ GeV
and within $|\eta|<2.4$ while $\Ntrk$ is the number of particles produced with $p_T>0.4$ GeV and within $|\eta|<2.4$.
The results from DCCI2 with full simulations show positive $\ctwofour$ for the entire $\Nch$ while much more statistics are required to make a solid conclusion.
To resolve this in detail, I discuss what is happening in DCCI2 with two-particle correlations in 2-D map in Sec. <ref>.

§.§ Multiplicity dependence in $Pb$+$Pb$ collisions

Figure <ref> shows
second-, third-, and fourth- order of anisotropic flow coefficients obtained from two-particle correlation as functions of the number of produced charged particles in $Pb$+$Pb$ collisions at = 2.76 TeV are shown for each.
There is a clear difference in behavior of $v_n\{2\}$ between core and corona components for $Pb$+$Pb$ collisions while there is not for $p$+$p$ collisions.
Let me first focus on results of $\vtwtw$.
From a comparison between the results with and without hadronic rescatterings, one can tell that a slight enhancement of $\vtwtw$ comes from generation of elliptic flow in the late hadronic rescattering stage <cit.>.
As an over all tendency, $\vtwtw$ of core components shows maximum around
semi-central events, which is around $\Nch\approx 300-400$, and decreases from semi-central to central, which is similar to the tendency observed in experimental data <cit.>.
This can be understood as follows:
generated QGP fluids are expected to have an almond shape when one sees it in transverse plane around this centrality regime. The almond-shaped
QGP has azimuthally anisotropic geometry, which leads to azimuthal anisotropy
of pressure. Since fluid velocity is generated due to the pressure gradient, the azimuthal momentum anisotropy is inherited by the emitted particles from the pushed-out fluids.
On the other hand, $\vtwtw$ from corona components are quite small compared to those from core components.
The results show slow decrease from low to high multiplicity events.
Here again, one can see the correction from corona components in the comparison between the core result and the inclusive result in the case without hadronic rescatterings.
The correction from corona components is found to be $\approx 15$-$30\%$ below $N_{\mathrm{ch}} \approx 400$, which originates from the small peak seen at very low $p_T$ region in the $p_T$ spectra in Fig. <ref> (b) [The leftmost point of the contribution from core components is slightly shifted to large $N_{\mathrm{ch}}$ since there are some events  in which one cannot calculate two-particle cumulants due to less than two charged particles from the core parts are measured in a given kinematic window in this $N_{\mathrm{ch}}$ bin. Therefore the event average of $N_{\mathrm{ch}}$ for core components is biased to larger $N_{\mathrm{ch}}$.].
This suggests that one would need to incorporate corona components in hydrodynamic frameworks to extract transport coefficients from comparisons with experimental data.

Second, third, and fourth order of anisotropic flow coefficients obtained from two-particle correlation for charged hadrons as functions of the number of produced charged particles in $Pb$+$Pb$ collisions at = 2.76 TeV are shown in upper left, upper right, and bottom figures. Results from full simulations (orange squares) and simulations without hadronic rescatterings (yellow triangles) are shown.  Results only from core (red diamonds) and corona (blue circles) components as a breakdown of results of simulations without hadronic rescatterings are shown simultaneously.

In both $p$+$p$ and $Pb$+$Pb$ results,
there are two factors that would give a finite anisotropy in corona components: color reconnection and feed-down from surviving partons.
The color reconnection effect implemented in default 8 and 8 Angantyr can arise collectivity <cit.>.
With the color reconnection, dense color strings formed due to multiparton interactions interact with each other and eventually induce flow-like behavior of final hadrons.
Its effect can be enhanced due to more multiparton interactions in initial parton generation with DCCI2 compared to default 8 and 8 Angantyr.
The detailed discussion on multiparton interactions in initial parton generation is made in Sec. <ref>.
Under the dynamical core–corona initialization, partons originating from hard scatterings and emitted in the back-to-back directions tend to survive.
In contrast, soft partons, which originate from multiparton interactions and are randomly directed, tend to be converted into fluids.
Since the low $p_T$ charged hadrons come from such surviving partons through string  fragmentation, $\vtwtw$ of corona components could reflect that of their parents.
As a result, the corona components show larger anisotropy compared to results from the default 8 <cit.> and 8 Angantyr.

I would like to mention some points from the comparison of anisotropic  flows among different orders.
In Fig. <ref>,
note that the scale of vertical axes are different in the three figures for each.
The overall tendency is quite similar for all of $\vtwtw$, $\vthree$, and $\vfour$.
From comparisons between full simulations and simulations without hadronic rescatterings, one sees that the effect of hadronic rescatterings are comparatively manifested in $\vtwtw$ rather than $\vthree$ and $\vfour$.
It is also noticeable that the fraction of anisotropy of corona components
against that of core components becomes closer in fourth- compared to second- or third- order.

Comparisons among second (green diamonds), third (navy squares), and fourth- (blue stars) order of anisotropic flow coefficients obtained from two-particle correlation for charged hadrons as functions of the number of produced charged particles in $Pb$+$Pb$ collisions at = 2.76 TeV from full simulations are shown.

To see the overall behavior of anisotropic flows in different orders,
a comparison of full simulation results is made as a function of the number of produced charged particles.
As it is the case in other existing hydrodynamic models,
the magnitude decreases as the order of anistoropic flows ascends
similarly to the $p$+$p$ collision results in Fig. <ref>.
Again, it would be worth mentioning that, within DCCI2, initial momentum anisotropy is embedded in initial matter profile due to dynamically deposited energy and momentum from initial partons.
Hence, it can be possible that the anisotropy pronounced even in $Pb$+$Pb$ collisions
might reflect initially generated anisotropy.

Second order of anisotropic flow coefficients obtained from two-particle correlation for charged hadrons with and without eta gap ($\Delta \eta > 3.0$) as functions of the number of produced charged particles in $Pb$+$Pb$ collisions at = 2.76 TeV are shown in left and right figures. Results from simulations without hadronic rescatterings (yellow squares) and its breakdown into core (red triangles) and corona (blue diamonds) components are shown simultaneously.

Second order of cumulant obtained from four-particle correlation for charged hadrons as a function of the number of produced charged hadrons in $Pb$+$Pb$ collisions at = 2.76 TeV. Results from full simulations (orange squares) and simulations without hadronic rescatterings (yellow triangles) are shown.  Results only from core (red diamonds) and corona (blue circles) components as a breakdown of results of simulations without hadronic rescatterings are shown simultaneously.

As I discuss the eta gap dependence of anisotropic flows in $p$+$p$ collision  in Fig. <ref>,
the corresponding results in $Pb$+$Pb$ collisions are shown in Fig. <ref>.
In order to take a large eta gap, $\vtwtw$ is obtained with particles produced within $|\eta| < 3.0$.
In the left figure, $\vtwtw$ without eta gap is shown, while that with eta gap of $\Delta \eta > 3.0$ is shown in the right.
The results from simulations without hadronic rescatterings,
which include both core and corona components, show a slight decrease with eta gap compared to the one without eta gap especially in low multiplicity events, which is although negligible.
On the other hand, its breakdowns, the results both from core and corona components, show negligible decrease as well
despite of the expectation that anisotropic flow from corona components would be suppressed by imposing eta gap.

Figure <ref> shows $\ctwofour$ as a function of the number of produced charged particles in $Pb$+$Pb$ collisions at = 2.76 TeV.
There is no eta gap imposed because we know that the effect of eta gap is not so significant from results in Fig. <ref>.
Compared to the corresponding $p$+$p$ collision results,
one sees trends of each component clearly.
First, from the comparison of results between full simulations and simulations without hadronic rescatterings,
one sees that the absolute value of $\ctwofour$ is enhanced in full simulations due to the effect of hadronic rescatterings,
which indicates that the momentum anisotropy is grown due to the hadronic interactions.
Second, the core component shows negative $\ctwofour$, which manifests clear signal of hydrodynamic behavior.
On the other hand, one from corona components is zero consistent for the entire multiplicity range: non-flow as back-to-back 2-particle correlations is successfully subtracted in $\ctwofour$.
Here, the point that we should focus on is $\ctwofour$ from results without hadronic rescatterings.
As we have been discussed, comparison of observables between core components and results without hadronic rescatterings extracts the corrections to pure hydrodynamic results
due to the existence of corona components in the system, especially in $Pb$+$Pb$ collisions.
In Fig. <ref>, one sees the corona correction:
the absolute value of $\ctwofour$ of core components is decreased to the one from results without hadronic rescatterings while the result from corona components shows zero-consistent $\ctwofour$.
Simply speaking, multi-particle correlation is the correlation per the number of permutation of $m$-particle.
Even if there are particles that does not contribute to correlation, in this case they are the corona components, the permutation of $m$-particle in the denominator is calculated from all particles including corona components.
From this result, one lesson that we can learn is that
anisotropic flows obtained from pure hydrodynamic models cannot directly be compared with
experimental data if there is a certain portion of non-equilibrium components.

§.§ Transverse momentum dependence

We have been, so far, seeing $p_T$ integrated anisotropic flows.
Especially in $Pb$+$Pb$ collisions, we found that the corona components at low $p_T$ regime
make an non-negligible effects on anisotropic flows from core components.
To investigate characteristics of anisotropic flow from each core and corona components,
I analyze anisotropic flows differentiating by $p_T$.

Second order of anisotropic flow coefficients obtained from two-particle correlation with and without eta gap as functions of transverse momentum from INEL$>0$ $p$+$p$ collisions at [proton]=13 TeV are shown in left and right figures, respectively.
Results from full simulations (orange diamonds) and simulations without hadronic rescatterings (yellow squares) are shown.  Results only from core (red triangles) and corona (blue circles) components as a breakdown of results of simulations without hadronic rescatterings are shown simultaneously.

Figure <ref> shows second order of anisotropic flow coefficients obtained from two-particle correlation as a function of $p_T$ from INEL$>0$ $p$+$p$ collisions at [proton] = 13 TeV.
The left and right figures show results without and with eta gap, respectively, while the eta gap imposed here is $\Delta \eta > 1.0$.
Note that anisotropic flow is obtained within each $p_T$ bin.
As an overall tendency, firstly, there is no significant difference of the results among full simulations, simulations without hadronic rescatterings, and core components while corona components show slightly larger $\vtwtw$.
On the other hand, once eta gap is imposed on $\vtwtw$ to eliminate short range correlation, $\vtwtw$ as a function of $p_T$ show significant suppress for intermediate to high $p_T$ regime,
and which is seen in both core and corona components.
Thus, it can be said that the distance of correlation becomes short in rapidity regardless of whether they are core or corona components, which indicates that there can be a strong anisotropic correlation between core and corona.

Second order of anisotropic flows obtained from two-particle correlation as a function of transverse momentum in $0-10\%$ and $40-60\%$ in centrality without eta gap are shown in left and right of upper panels, respectively.
Corresponding results with eta gap ($\Delta \eta > 1.0$) are shown in lower panels.
Results from full simulations (orange diamonds) and simulations without hadronic rescatterings (yellow squares) are shown.  Results only from core (red triangles) and corona (blue circles) components as a breakdown of results of simulations without hadronic rescatterings are shown simultaneously.

Second order of anisotropic flow coefficients obtained from two-particle correlation with eta gap ($\Delta \eta > 1.0$) as a function of transverse momentum from $Pb$+$Pb$ collisions in $0-5\%$ (orange squares) and $30-40\%$ (red triangles) in centrality. Comparisons with the ALICE experimental data <cit.> (open black symbols) for each centrality is made.

Figure <ref> shows the corresponding results with above in $Pb$+$Pb$ collisions.
Upper and lower figures are results without and with eta gap, respectively.
Results in 0-20% and 40-60% in centrality are shown in left and right panels, respectively.
One sees that there is a clear different behavior in the corona components while statistics is not sufficient for results with eta gap:
for results without eta gap, $\vtwtw$ from corona component is smaller in low $p_T$ (around $p_T<2$) while one shows rapid increase towards larger $p_T$ regime compared to core component.
We can infer some points from the above results:
first, the corona components in very low $p_T$, that behave as a correction to core components, has smaller anisotropy in momentum space.
Second, the back-to-back 2-particle anisotropy from non-flow is manifested especially above $p_T \approx 2$ GeV in results without eta gap.

Finally I show data to model comparison of $\vtwtw$ with eta gap ($\Delta \eta> 1.0$) as a function of $p_T$ in 0-5% and 30-40% in centrality of $Pb$+$Pb$ collisions in Fig. <ref>.
The results from DCCI2 clearly overestimate the experimental data in both centrality.
The overestimation is expected to be reconciled by introducing shear viscosity in hydrodynamic calculations because shear stress induced by gradient of flow in a perpendicular direction of the direction of flow suppress the growth of anisotropy.
I admit that the introduction of such dissipative currents in hydrodynamic simulation is necessary to reproduce the experimental data and it is indispensable to explore the QGP transport.
I leave this update for future work.

§.§ Ridge structure

In anisotropic flow coefficients,
we can quantify amplitude of each order of anisotropy.
In this section, I discuss details of two-particle correlation by plotting as a functions of $\Delta eta$ and $\Delta \phi$.
Details of the technique used in this calculation is explained in Appendix. <ref>.

As I discussed in Sec. <ref>, the ridge structure,
two-particle correlation at $\Delta \approx 0$ along wide $\Delta \eta$ range,
is a proxy of hydrodynamic response due to QGP formation.
The ridge structure has been observed not only in heavy-ion but also in $p$+$p$ collisions.
Thus, the general question is, “what is the origin of the ridge structure observed in the experimental data of $p$+$p$ collisions at high multiplicity?”.
Giving an answer within DCCI2 is one of the big goals in this study.
However, as we see from Sec. <ref> to
DCCI2 has not reached to the level in which I can compare the anisotropic flow results from the model with experimental data.
Admitting that fact,
I discuss whether we can see the ridge structure in high-multiplicity events in $p$+$p$ collisions.
In addition, I show if the simple expectation:
core components as hydrodynamic response should give a long range correlation
and that corona component should give a short range correlation,
is realized within DCCI2
in the followings.

§.§.§ $p$+$p$ collision results

Two-particle correlations as functions of $\Delta \eta $ and $\Delta \phi$ obtained from charged particles produced within $|\eta|<1.5$ and with $1<p_T<2$ GeV for both triggers and associates in 0-0.95% of V0M multiplicity class of $p$+$p$ collisions at [proton]=13 TeV. Results from full simulation (left), simulations without hadronic rescatterings (right) are shown at the top, while breakdowns of produced particles from simulations without hadronic rescatterings into core (left) and corona (right) components are shown at the bottom. The centrality and the $p_T$ range are determined based on the method used in the ALICE paper <cit.>.

Two-particle correlations as functions of $\Delta \eta $ and $\Delta \phi$ obtained from charged particles produced within $|\eta|<1.5$ and with $1<p_T<3$ GeV for both triggers and associates in 0-100% of V0M multiplicity class of $p$+$p$ collisions at [proton]=13 TeV. Results from full simulation (left), simulations without hadronic rescatterings (right) are shown at the top, while breakdowns of produced particles from simulations without hadronic rescatterings into core (left) and corona (right) components are shown at the bottom. The centrality and the $p_T$ range are determined based on the method used in the ALICE paper <cit.>.

I discuss multiplicity dependence of two-particle correlation in $p$+$p$ collisions at [proton]=13 TeV
with 2 different centrality/multiplicity cuts: the events are classified following the ALICE <cit.> and CMS collaborations <cit.>, respectively.
First I start with the ALICE kinematic settings.
Figure <ref> shows
two-particle correlation functions of $\Delta \eta$ and $\Delta \phi$
from high-multiplicity events (0-0.95% of V0M multiplicity class) of $p$+$p$ collisions at [proton]=13 TeV. Results from full simulation and simulations without hadronic rescatterings are shown in left and right at the top, while breakdowns of produced particles from simulations without hadronic rescatterings into core and corona components are shown in left and right at the bottom.
On the other hand, the corresponding results in low-multiplicity events (0-100% of V0M centrality)
are shown in Fig. <ref>.
It should be noted that, to gain statistics, two-particle correlations are obtained from charged particles produced within $|\eta|<1.5$ instead of $|\eta|<0.8$, the eta range covered by the detector at mid-rapidity introduced in the ALICE collaboration.
Hence the results shown here cannot be compared directly to the results from the ALICE experiment.
From the comparisons of full simulations between the two centrality classes,
one sees that the ridge structure, correlation around $\Delta \phi \approx 0$ along $\Delta \eta $ axis, is slightly enhanced in high-multiplicity class:
for instance the small peak around $\Delta \phi \approx 0$ and $\Delta \eta \approx -1$ is seen more clearly in high-multiplicity compared to low-multiplicity events.
One also sees that there is not so significant difference of results between full simulations and simulation without hadronic rescatterings, which is consistent with the observation for the $p_T$ integrated $\vtwtw$.
The difference between the two-particle correlation of core and corona components is
relatively clear in both multiplicity classes.
The near side ridge-like structure appears in core components in both multiplicity classes
while such structure is not that manifested in corona components.
With these results, I conclude that the ridge-like structure in near-side can be seen up to $|\Delta \eta| \approx 1$, which originates from the core components, within DCCI2 with the current parameter set.
However, it should be noted that the away-side peak at $\Delta \approx \pi$ stretched along $\Delta \eta$ is so small that the structure is almost absent
in high-multiplicity events of full simulation results
while that is observed in experimental data.
This means that back-to-back correlation of mini-jets that originally should be produced due to momentum conservation dismiss due to the dynamical core–corona initialization.
As we see in the comparison of $p_T$ spectra between DCCI2 and experimental data in Sec. <ref>,
results from DCCI2 underestimate experimental data above $p_T\approx 1$ GeV at the highest multiplicity class.
Thus, these results in total suggest that there is too much energy loss of partons
at this $p_T$ regime,
which should be reconciled in future updates of DCCI.

Two-particle correlations as functions of $\Delta \eta $ and $\Delta \phi$ obtained from charged particles produced within $|\eta|<2.4$ and with $1<p_T<3$ GeV for both triggers and associates in high-multiplicity $p$+$p$ collisions ($60<\Ntrk<85$) at [proton]=13 TeV. Results from full simulation (left), simulations without hadronic rescatterings (right) are shown at the top, while breakdowns of produced particles from simulations without hadronic rescatterings into core (left) and corona (right) components are shown at the bottom. The same kinematic cuts and event classification used in the CMS paper <cit.> are adopted.

Two-particle correlations as functions of $\Delta \eta $ and $\Delta \phi$ obtained from charged particles produced within $|\eta|<2.4$ and with $1<p_T<3$ GeV for both triggers and associates in almost minimum-bias $p$+$p$ collisions ($0<\Ntrk<\infty$) at [proton]=13 TeV. Results from full simulation (left), simulations without hadronic rescatterings (right) are shown at the top, while breakdowns of produced particles from simulations without hadronic rescatterings into core (left) and corona (right) components are shown at the bottom. The same kinematic cuts and event classification used in the CMS paper <cit.> are adopted.

Second, I investigate the corresponding results with the CMS kinematic cuts.
Figure <ref> and <ref> show
two-particle correlation functions as functions of $\Delta \eta$ and $\Delta \phi$
obtained from charged particles produced within $|\eta|<2.4$ and with $1<p_T<3$ GeV
in high-multiplicity ($60<\Ntrk<85$) and minimum-bias ($0<\Ntrk<\infty$) events, respectively.
Note that the range of $\Delta \eta$ is different from the one in Figs. <ref> and <ref>
because the two-particle correlation function is now obtained within wider eta range.
Only in high-multiplicity events in Fig. <ref>, one sees that there is the double peak structure appearing at $\Delta \phi \approx \pi$ in the results of full simulations and simulations without hadronic rescatterings,
which is not observed in experimental results.
Because the double peak structure is seen in results of core components but corona components,
it can be expected that the structure seen in full simulations mostly comes from the correlation in core components.
On the other hand, the results of minimum-bias events in Fig. <ref>
show that the strong double peak structure that we see in Fig. <ref> is diluted.
The result of core components shows much smaller correlation only in near-side by 2-order compared to those from corona components.
It should also be mentioned that the hadronic rescatterings enhances correlations in both multiplicity events,
which indicates that hadronic interactions pronounce flows
even in $p$+$p$ collisions.

Kinematic cut dependence of two-particle correlations as functions of $\Delta \eta $ and $\Delta \phi$ obtained from charged particles produced within $|\eta|<2.4$ in minimum-bias (events with at least one charged particle produced within $|\eta|<1.0$) $p$+$p$ collisions at [proton]=13 TeV.
The two-particle correlation is obtained for associates limited to core components by triggering corona components. (Left) Transverse momentum ranges of both triggers and associates are $1.0<p_T<3.0$ GeV. (Middle, right) Trigger-particle kinematic cuts dependence of distribution of core components with regarding to corona components.

Figure <ref> shows
kinematic cut dependence of two-particle correlation function as functions of $\Delta \eta $ and $\Delta \phi$ obtained from charged particles produced within $|\eta|<2.4$ in minimum-bias (events with at least one charged particle produced within $|\eta|<1.0$) $p$+$p$ collisions at [proton]=13 TeV.
The two-particle correlation is obtained for associates limited to core components by triggering corona components. Transverse momentum ranges of both triggers and associates are $1.0<p_T<3.0$ GeV in left figure.
Comparing the results in left figure and Fig. <ref>,
the double peak structure and the near-side peak are more obvious in the former result.
Thus, this means that within this $p_T$ range there is a tendency of that fluidization takes place along corona mini-jets when a mini-jet is produced
with a certain rapidity distance.
In the fluidization mechanism in Fig. <ref>,
partons with high $p_T$ tend to keep their energy and momentum,
hence the probability that
mini-jets produced with almost no rapidity gap are fluidized
is highly suppressed because the larger $p_T$ a particle obtains, the more mid-rapidity it is produced.
The middle and right figures show comparisons of trigger-particle kinematic cuts dependence of distribution of core components with regarding to corona components.
As one sees, the peaks seen in the left figure remain in the middle but disappear in the right.
Thus, the fluidization explained just above produces low $p_T$ core components.

§.§.§ $Pb$+$Pb$ collision results

Two-particle correlations as functions of $\Delta \eta $ and $\Delta \phi$ obtained from charged particles produced within $|\eta|<2.4$ and with $3.0<p_T<3.5$ for triggers and $1.0<p_T<1.5$ for associates in 5-10% centrality of $Pb$+$Pb$ collisions at =2.76 TeV. Results from full simulation (left), simulations without hadronic rescatterings (right) are shown at the top, while breakdowns of produced particles from simulations without hadronic rescatterings into core (left) and corona (right) components are shown at the bottom. The same kinematic cuts used in the CMS paper <cit.> are adopted.

Two-particle correlations as functions of $\Delta \eta $ and $\Delta \phi$ obtained from charged particles produced within $|\eta|<2.4$ and with $3.0<p_T<3.5$ for triggers and $1.0<p_T<1.5$ for associates in 20-40% centrality of $Pb$+$Pb$ collisions at =2.76 TeV. Results from full simulation (left), simulations without hadronic rescatterings (right) are shown at the top, while breakdowns of produced particles from simulations without hadronic rescatterings into core (left) and corona (right) components are shown at the bottom. The same kinematic cuts used in the CMS paper <cit.> are adopted.

Contrast to $p$+$p$ collisions, it has been understood and concluded that ridge structures observed in $Pb$+$Pb$ collisions originate from hydrodynamic responses against geometric anisotropy of initial profile of matter produced just after a collision.
Here, I discuss whether I can conclude so even in DCCI2
where both core and corona components are incorporated and initial profiles are dynamically generated, which are different from other hydro-based dynamical models.
Figures <ref> and <ref> show two-particle correlations as functions of $\Delta \eta $ and $\Delta \phi$ obtained from charged particles produced within $|\eta|<2.4$ and with $3.0<p_T<3.5$ for triggers and $1.0<p_T<1.5$ for associates in 5-10% and 20-40% of centrality from $Pb$+$Pb$ collisions at =2.76 TeV. The same kinematic cuts used in the CMS paper are adopted.
In contrast to $p$+$p$ collision results, ridge structure spanning long range in $\Delta \eta$ is seen in $Pb$+$Pb$ collisions, and that is more clearly seen in peripheral (20-40%) compared to central (5-10%) collisions.
This is consistent with the picture that almond-shaped matter generated early time in peripheral collisions produces more second order of anisotropic flow in momentum as a hydrodynamic response compared to central collisions.
It should be also noted that the ridge structure is seen only in core but in corona components.

§ EVOLUTION OF TRANSVERSE ENERGY

As shown in Fig. <ref>, I reproduced centrality dependence of charged particle multiplicity in $p$+$p$ and $Pb$+$Pb$ collisions within DCCI2.
Although the default 8 (or Angantyr model in heavy-ion modes) works reasonably well, reproduction of multiplicity within DCCI2 can be attained only after the considerable change of a parameter $p_{\mathrm{T0Ref}}$ from its default value as mentioned in Sec. <ref>.
A nontriviality in DCCI2 stems from different competing mechanisms of how the transverse energy changes during the evolution of the system.
In this subsection, I discuss the effects of string formation/fragmentation and longitudinal $pdV$ work on the transverse energy and explain why I needed to change this parameter in DCCI2.

The transverse energy per unit rapidity $dE_{T}/d\eta$ is a basic observable in high-energy nuclear collisions and contains rich information on the dynamics of an entire stage of the reactions.
The transverse energy changes mainly in the initial and the expansion stages of the reactions.
In the initial stage, the two energetic hadrons and/or nuclei form color flux tubes between them as they pass through each other.
The chromo-electric and magnetic fields in the color flux tubes possess the energy originating from the kinetic energy of colliding hadrons or nuclei.
The decays of color flux tubes into partons and subsequent rescatterings among them are supposed to lead to the QGP formation <cit.>.
Thus, how much energy is deposited in the reaction region is a fundamental problem of the QGP formation and depends on the initial dynamics of high-energy nuclear collisions.
On the other hand, in the expansion stage, the $pdV$ work associated with the longitudinal expansion after the QGP formation reduces the energy produced in the initial reaction region <cit.>.
The amount of reduced energy is sensitive to viscosity and other transport properties of the QGP <cit.>.
In addition, it can be possible to access the information on the initial state <cit.>.
Therefore $dE_{T}/d\eta$ can be a good measure to scrutinize modeling in the initial and the expansion stages of the reactions.

In 8, partons are first generated through hard scatterings and then, together with partons from initial and final state radiations, form hadron strings which eventually fragment into hadrons.
The transverse energy per unit rapidity of final hadrons is always larger than that of initially generated partons around midrapidity.
To understand this enhancement around midrapidity, suppose a hadronic string formed from a di-quark in the forward beam rapidity region and a quark in the backward beam rapidity region as an extreme case.
Although the partons lie only around beam rapidity regions and the transverse energy of them vanishes around midrapidity, that string fragments into hadrons almost uniformly in rapidity space.
Thus, the emergence of the transverse energy at midrapidity is a consequence of the formation of a color string between such partons around beam rapidity.
Since parameters in 8 are so tuned to reproduce the final hadron spectra, the initial parameters are highly correlated with parameters in the fragmentation as a whole.
Therefore the default parameter set should not be used if the subsequent hydrodynamic evolution, which reduces the transverse energy from its initial value of generated partons, is incorporated in DCCI2.

(Color Online)
Pseudo-rapidity distribution of transverse energy in (a) minimum-bias $p$+$p$ collisions at $\sqrt{s} = 7$ TeV
and (b) minimum-bias $Pb$+$Pb$ collisions at = $2.76$ TeV from DCCI2 and 8.
A comparison of transverse energy distribution between parton and hadron levels is made for results from DCCI2 and 8.
Results from the parton level in DCCI2 (green crosses) and 8 (blue pluses) and the ones from the hadron level (red crosses) in DCCI2 and 8 (orange pluses) are shown for comparison. For $Pb$+$Pb$ collisions, 8 Angantyr is used to obtain the results.

The infrared cut-off of the MPI, $p_{T0}$, as a function of center-of-mass energy of nucleon-nucleon collisions. Comparisons of $p_{T0}$ in default value used in 8 (dashed blue line), $p$+$p$ (solid orange line), and $Pb$+$Pb$ collisions (solid green line) in DCCI2 are shown. The red circle is each $p_{T0}$ at [proton] = 7 TeV, which corresponds to $p_{\mathrm{T0Ref}}$.

Figure <ref> shows $dE_{T}/d\eta$
of the initial partons before the string hadronization or the dynamical initialization, and that of the hadrons in the final state from both 8 and DCCI2 in minimum-bias (a) $p$+$p$ collisions at $\sqrt{s} = 7 $ TeV and (b) $Pb$+$Pb$ collisions at = $2.76$ TeV.
For the results from 8 in Fig. <ref> (a),
the transverse energy per unit pseudorapidity, $dE_{T}/d\eta$, of the final hadrons is always larger than that of the initial partons in the whole rapidity region except around the beam rapidity.
Since the final hadron yield is dominated by the corona components in minimum-bias $p$+$p$ collisions in DCCI2, the above behavior is qualitatively the same as the ones from 8.
The difference of the results between DCCI2 and 8 appears in the absolute value of $E_T/d\eta$ in parton level,
which is a consequence of partial formation of QGP fluids.
To obtain the same amount of the transverse energy in the final state in DCCI2, the transverse energy must be deposited initially $\approx 1.5$ times as large as that in 8 at midrapidity to reconcile the reduction of transverse energy due to $pdV$ work.
This is exactly possible by considerably decreasing the parameter $p_{\mathrm{T0Ref}}$.

For the results from 8 Angantyr in Fig. <ref> (b),
the behavior of the transverse energy in $Pb$+$Pb$ collisions is again the same as in $p$+$p$ collisions.
In contrast, the transverse energy in initial state  from DCCI2 shows $\approx 3$ times as large as that in 8 Angantyr.
The reason is the same as mentioned on the results of DCCI2 in $p$+$p$ collisions.
The decrease of the transverse energy due to $pdV$ work is clearly seen in $Pb$+$Pb$ collisions.

The parameter $p_{\mathrm{T0Ref}}$ is used to regulate infrared divergence of the QCD cross section.
The actual infrared cut-off is  $p_{T0}$ which is given as a function of center-of-mass energy of nucleon+nucleon collisions, [proton].
The $p_{\mathrm{T0Ref}}$ gives the reference $p_{T0}$ value at a reference [proton], which is 7 TeV by default.
Figure <ref> shows the infrared cut-off $p_{T0}$ as a function of [proton].
Comparisons of $p_{T0}$ in default value used in 8 (dashed blue line), $p$+$p$ (solid orange line), and $Pb$+$Pb$ collisions (solid green line) in DCCI2 are shown. The red circle is each $p_{T0}$ at [proton] = 7 TeV, which corresponds to $p_{\mathrm{T0Ref}}$.
The collision energy dependence of $p_{T0}$ is calm within the range of the LHC energy. Thus, one can see that $p_{\mathrm{T0Ref}}$ determines the overall factor of $p_{T0}$.

One can also interpret $p_{\mathrm{T0Ref}}$ as a parameter $p_{\perp\mathrm{min}}$ to separate soft from hard scales [In actual simulations in , a parameter $p_{\mathrm{T0Ref}}$ provides a scale to make a smooth turnoff of hard scattering rather than the sharp separation <cit.>.], and controls the number of multiparton interactions in  <cit.>.
The smaller the separation scale $p_{\perp\mathrm{min}}$ is, the larger the number of multiparton interaction $\langle n_{\mathrm{MPI}}(p_{\perp\mathrm{min}}) \rangle = \sigma_{2\rightarrow 2}(p_{\perp\mathrm{min}})/\sigma_{\mathrm{nd}} $ is.
Here $\sigma_{2 \rightarrow 2}$ and $\sigma_{\mathrm{nd}}$ are the perturbative QCD $2\rightarrow 2$ cross section and the inelastic non-diffractive cross section, respectively.
By increasing $\langle n_{\mathrm{MPI}}(p_{\perp\mathrm{min}}) \rangle$ as decreasing $p_{\mathrm{T0Ref}}$, initial partons are generated more and bring the sufficient amount of transverse energy in the final hadron state as shown in Fig. <ref> (b).
In this work, I use $p_{\mathrm{T0Ref}}=$1.8 and 0.9 for $p$+$p$ and $Pb$+$Pb$ collisions as mentioned in Sec. <ref>, which are smaller than the default values,
2.28 and 2.0 for
in .

So far, we have found that
hydrodynamics and string fragmentation have different evolution of transverse energy. Thus, the multiplicity of final hadrons from such a two-component model is sensitive to a fraction of each component in a system.
To, at least, reproduce multiplicity in DCCI2, I need to change the parameter $p_{\mathrm{T0Ref}}$ from its default value.
However, as I mentioned in Sec. <ref> in chapter 3, the other parameter $\sigma_0$ in Eq. eq:four-momentum-deposition
has a non-trivial correlation with $p_{\mathrm{T0Ref}}$,
which means that I need to tune both the parameters at the same time.
Suppose that one firstly tries to reproduce multiplicity by tuning $p_{\mathrm{T0Ref}}$.
Since a small $p_{\mathrm{T0Ref}}$ gives a rise to the number of multiparton interactions and
deposited transverse energy is enhanced, it affects to enhance final hadron multiplicity.
On the other hand,
since the number of initial partons produced in midrapidity increases,
this causes more fluidization in dynamical core–corona initialization.
Once a fraction of the core is enhanced, the multiplicity of final hadrons can also decrease since the initial transverse energy deposited in midrapidity region is used for $pdV$ work.
As a result of competition between two effects, multiplicity cannot linearly enhance or decrease by decreasing or increasing $p_{\mathrm{T0Ref}}$.
Secondly, suppose that one tries to reproduce particle yield ratios as functions of multiplicity by tuning $\sigma_0$.
Since changing $\sigma_0$ means changing a fraction of core and corona, final multiplicity is easily altered, too.
This is the reason why I need to fix both parameters by taking into account multiplicity and particle yield ratios at the same time.

Note that, if I made viscous corrections in the hydrodynamic evolution, the resultant change of  $p_{\mathrm{T0Ref}}$ from its default value could have been modest due to the less reduction of transverse energy <cit.>, which is beyond the scope of the present paper, but which should be investigated in the future work.

A string melting version of A Multi Phase Transport (AMPT) model  <cit.> and the hydrodynamic models using it for generating initial conditions <cit.> avoid this issue of the transverse energy in an “ad hoc" way: The hadrons decaying from a  string are re-decomposed into their constitutive quarks and antiquarks, and then form high-energy density partonic matter.
Although it is possible to count the energy stored along the string contrary to considering the generated partons directly, this prescription lacks gluons from melting strings.
Therefore, I do not pursue this idea in the present paper.

CHAPTER: CONCLUSION

I studied the interplay between core and corona components establishing the DCCI2, which describes the dynamical aspects of core–corona picture under the dynamical initialization scheme from $p$+$p$ to $A$+$A$ collisions.
To develop the DCCI2, I put an emphasis on the reconciliation of open issues of dynamical models, mainly relativistic hydrodynamic models, toward a comprehensive description of a whole reaction of high-energy nuclear collisions.
One of the important achievements is to generate the initial profiles of hydrodynamics by preserving the initial total energy and momentum of the collision systems.
This is achieved by adopting hydrodynamic equations with source terms on initial partons obtained from 8, one of the widely accepted general-purpose Monte-Carlo event generators.
Consequently, in addition to the equilibrated matter (core) described by relativistic hydrodynamics, I also consider the existence of non-equilibrated matter (corona) through dynamical initialization with the core–corona picture.

I have updated DCCI that I established in the Master's thesis.
The updates include sophistication of four-momentum deposition of initial partons in dynamical core–corona initialization, samplings of hadrons from hypersurface of fluids with , hadronic afterburner for final hadrons from both fluids and non-equilibrated partons with a hadronic transport model , and modification on color string structures of non-equilibrated partons due to the co-existence with fluids in coordinate space.

Discussion on the interplay between core and corona components is made once
fixing major parameters so that our model reasonably describes multiplicity as a function of multiplicity or centrality class and omega baryon yield ratios to charged pions as functions of multiplicity.
First I extracted the fractions of core and corona components to the final hadron yields as functions of multiplicity and centrality classes. I found that, as increasing multiplicity, the core components become dominant at $\langle dN_{\mathrm{ch}}/d\eta \rangle_{|\eta|<0.5} \approx 20$, which roughly corresponds to the highest multiplicity classes in $p$+$p$ collisions at $\sqrt{s} = 7$ and $13$ TeV and $\approx80\%$ centrality class in $Pb$+$Pb$ collisions at $\sqrt{s_{NN}} = 2.76$ TeV.

Next, I showed the fractions of core and corona components in $p_T$ spectra with and without particle identification.
In the charged particle spectra from minimum-bias $Pb$+$Pb$ collisions, the fraction of core components is dominant below $p_T\approx 5.5$ GeV, while that of corona components is dominant above that.
Interestingly, I found that there was an enhancement in the fraction of corona contribution with $R_{\mathrm{corona}} \approx 0.2$ at most in $p_T \lesssim 1$ GeV even in minimum-bias $Pb$+$Pb$ collisions. From this, the fraction of the corona contribution is anticipated to increase in peripheral collisions.
This brings up a problem in all conventional hydrodynamic calculations in which low $p_T$ soft hadrons are regarded purely as core components.
Since the fraction of each component would exist finite for a wide range of multiplicity and, as a result, there should be interplay between them, I suggest that
both small colliding systems and heavy-ion collisions
should be investigated in a unified theoretical framework by incorporating both core and corona components.
The PID and centrality classification revealed the dynamics in more detail.
In central collisions, the radial flow is also enhanced
while the $p_T$ integrated core components increase.
As a result, the core component dominance at intermediate $p_T$ (around $1<p_T<3$ GeV) is clearly pronounced at central collisions.
Because heavier particles are pushed towards higher $p_T$ due to the radial flow,
the core dominance at the intermediate $p_T$ and the corona enhancement at the very low $p_T$ become obvious in heavier particle spectra.
the fraction of core components is almost $100\%$ from $0$-$60\%$ in intermediate $p_T$ for proton spectra, and the corona component reaches around $50\%$ at the very low $p_T$ from central to peripheral collisions.

It is known that transverse momentum spectra solely from hydrodynamics or hybrid (hydrodynamics followed by hadronic cascade) models do not perfectly reproduce the experimental data below $p_T \approx 0.5$ GeV <cit.>,
although hydrodynamics is believed to provide a better description in the low $p_T$ region in general.
To investigate if the deviation between pure hydrodynamic results and the data in the low $p_T$ region can be filled with the corona components,
I made comparisons of $p_T$ spectra at the very low $p_T$ (smaller than $p_T \approx 1$ GeV) between experimental data and DCCI2 results from full simulations and those only from core components with hadronic rescatterings.
The slopes of $p_T$ spectra at the very low $p_T$ observed in the experiment cannot be reproduced only with core components with hadronic rescatterings, and the full simulation results show a better description of the slopes by including corona component.
Therefore I conclude that there is possibly  the “soft-from-corona” components in heavy-ion collisions,
and which would be indispensable to be considered for an adequate interpretation of the particle productions obtained from experimental data.

To investigate the effects of co-existence of core and corona components on observables, I showed $\langle p_T \rangle$ and $v_{2}\{2\}$ as functions of $N_{\mathrm{ch}}$.
In particular, in $Pb$+$Pb$ collisions,
I found that the finite contribution of corona components at midrapidity gives a certain correction to the results obtained purely from core components, which is described by hydrodynamics.
The correction is $\approx 5$-$11$% for $\langle p_T \rangle$ below $N_{\mathrm{ch}}\approx 200$,
while it is
$\approx 15$-$30$% for $\vtwtw$ below $N_{\mathrm{ch}}\approx 400$.
The former correction leads to the reasonable agreement of $\langle p_T \rangle$ with the experimental data.
This suggests that one needs to incorporate corona components in hydrodynamic frameworks to extract transport coefficients from comparisons with experimental data.

In the results of $\ctwofour$ as a function of $\Nch$ obtained from $Pb$+$Pb$ collisions,
I saw that the $\ctwofour$ obtained from core components is diluted due to the existence of the corona components, while the corona components show the zero-consistent $\ctwofour$.
This result indicates that anisotropic flows obtained from experimental data cannot be directly compared with pure hydrodynamic calculations if there is a certain portion of non-equilibrium components.

I also explored the effects of radial flow based on violation of $m_T$ scaling with hadron rest mass by classifying events into high- and low-multiplicity ones.
Noteworthy, I found that it is difficult to discriminate the radial flow originating from hydrodynamics from collectivity arising from color reconnection in 8.

The investigation of ridge structure within DCCI2 was performed to reveal its origin.
The analysis of $p$+$p$ collision results is performed with the event classification used in the ALICE and CMS experiment, respectively.
With the method used in the ALICE experiment,
I saw that the ridge-like structure appears in high-multiplicity events, and it originates from the core components.
On the other hand, such a ridge-like structure was not seen with the method used in the CMS experiments, which is inconsistent with the CMS experimental results.
I also checked that there is a strong momentum correlation between core and corona components,
which would be a consequence of the dynamical initialization.

I finally discussed the evolution of transverse energy in the DCCI2.
In string fragmentation, final transverse energy is larger than initial transverse energy as producing hadrons around midrapidity. While in hydrodynamics, transverse energy just decreases from its initial value during the evolution due to the longitudinal $pdV$ work.
To obtain the same amount of transverse energy in the final state in DCCI2 with default 8 Angantyr in minimum-bias $Pb$+$Pb$ collisions, it is necessary to have three times larger initial transverse energy than the one of default 8 Angantyr.

For a more quantitative discussion on transport properties of the QGP fluids, I admit an absence of viscous corrections to fully equilibrium distribution in our analysis.
I leave this as one of our future works.
Nevertheless, I emphasize that the corrections from corona components mean the ones from “far from” equilibrium components which should exist nonetheless and would more significantly affect the final hadron distributions than the viscous corrections.

With this model, I anticipate that it would be interesting to explore
planned $O$+$O$ collisions at LHC <cit.> since the collision system can provide data around the “sweet spot”
in which
the core components are to be dominant, and yet corrections from corona components cannot be ignored at all
In addition, investigation of strangeness enhancement in forward or backward regions might provide some insights into ultra-high-energy cosmic ray measurements <cit.>.
Incorporation of a dynamical description of kinematic and chemical pre-equilibrium stage <cit.> and investigation of medium modification of jets <cit.> are in our interests as well.
I leave the discussion on those topics as future work.

§ QUANTUM CHROMO DYNAMICS (QCD)

Quantum Chromo Dynamics (QCD) is a fundamental theory to describe strong force.

There are two notable characteristics which are "confinement of color" and "asymptotic freedom".

§.§ Color confinement

It is said to be that quarks cannot be extracted as one unit, since they are confined to form color singlet.

§.§ Asymptotic freedom

Asymptotic freedom of quarks was proposed by Wilchek and Gross.

§.§ Sigma model

Then one can rewrite the QCD Lagrangian density as follows,

By expressing the mass term with a combined field of isoscalar $\sigma$ and isovector $\vec{\pi}$, one can get a form of Lagrangian in a chiral-symmetry form.

The potential term, for the case that $\nu^2>0$, becomes a wine-bottle shape sitting on $\sigma$-$\vec{\pi}$ plane. Thus, the origin is no longer a lowest energy state, but the $\sigma^2+\vec{\pi}^2=\nu^2$ is. Now one can choose a particular vacuum from this edge of the wine bottle. By choosing a vacuum at $\sigma=\nu$ and $\vec{\nu}=0$ and shifting the Lagrangian with  $\sigma\rightarrow\sigma+\nu$ and $\vec{\pi}\rightarrow\vec{\pi}$, one can get the following form.
This Lagrangian is now not invariant under axial transformation, which means, chiral symmetry is spontaneously broken.
The important thing to notice here is that a mass term of $\vec{\pi}$ does not appear while that of $\sigma$ does. From this consequence, it is said that the isovector $\vec{\pi}$ is a Nambu-Goldstone boson, mass-less boson which appears due to the spontaneous chiral symmetry breaking.

Milne coordinate is adopted in hydrodynamic simulations. Since Milne coordinate is so called a curved space, it is convenient to have some knowledge on general relativity.

§ BASICS

§.§ Metric and transformation matrix

\begin{align}
ds^2 = \eta_{\mu\nu} dx^{\mu} dx^{\nu},
\end{align}

\begin{align}
\eta_{\mu\nu} =
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & -1 & 0 & 0 \\
0 & 0 & -1 & 0 \\
0 & 0 & 0 & -1
\end{pmatrix}
\end{align}

\begin{align}
ds^2 = g_{\alpha\beta} d\tilde{x}^{\alpha} d\tilde{x}^{\beta},
\end{align}

\begin{align}
g_{\alpha\beta} =  \frac{\partial x^{\mu}}{\partial\tilde{x}^{\alpha} }\frac{\partial x^\nu}{\partial\tilde{x}^{\beta}} \eta_{\mu\nu},
\end{align}

\begin{align}
d\bm{x} = dx^\mu \bm{e}_\mu = d\tilde{x}^\alpha \tilde{\bm{e}}_\alpha
\end{align}

\begin{align}
\bm{e}_\mu = \epsilon^\alpha_\mu \tilde{\bm{e}}_\alpha
\end{align}

\begin{align}
d\tilde{x}^\alpha = \epsilon^\alpha_\mu dx^\mu
\end{align}

\begin{align}
\epsilon^\alpha_{\mu} = \frac{\partial\tilde{x}^\alpha}{\partial x^\mu}
\end{align}

\begin{align}
ds^2 = (dx^\mu \bm{e}_\mu) \cdot (dx^\nu \bm{e}_\nu)  = dx^\mu dx^\nu g_{\mu\nu}
\end{align}

\begin{align}
g_{\mu\nu} = \bm{e}_{\mu} \cdot \bm{e}_\nu
\end{align}

§.§ Christoffel symbol

\begin{align}
\Gamma^\alpha_{\mu\beta} = \frac{\partial \tilde{x}^\alpha}{\partial x^\nu} \frac{\partial^2 x^\nu}{\partial x^\mu \partial \tilde{x}^\beta}
\end{align}

\begin{align}
\Gamma = \Gamma_{\beta \gamma \alpha } g^{\delta \alpha} = \frac{1}{2} g^{\gamma \alpha} (g_{\gamma \beta, \delta} - g_{\delta \gamma , \beta}+g_{\beta \gamma , \delta})
\end{align}

\begin{align}
A^\mu_{;\mu} &= A^{\mu}_{,\mu} + \Gamma^{\mu}_{\mu\xi} A^\xi \nonumber \\
& = \frac{1}{\sqrt{-g}} \frac{\partial}{\partial x^\mu} (\sqrt{-g}A^\mu)
\end{align}

Note that

\begin{align}
\det |\eta_{\mu\nu}| &= \det \left|\frac{\partial x^\mu}{\partial \tilde{x}^\alpha} \frac{\partial x^\nu}{\partial \tilde{x}^\beta} g_{\alpha \beta} \right| \nonumber \\
\Leftrightarrow   -1 &= \det \left| \frac{\partial x}{\partial \tilde{x}}\right|^2 \det |g| \nonumber \\\
\Leftrightarrow &  \det |g| \equiv g <0
\end{align}

\begin{align}
B^{\mu\nu}_{;\mu} = \frac{1}{\sqrt{-g}}\partial_\mu (\sqrt{-g}B^{\mu\nu}) + \Gamma_{\xi\mu}^\nu B^{\mu\xi}
\end{align}

§ MILNE AND CARTESIAN

CHAPTER: ANALYSIS METHODS

Details of analysis, which I skipped explanations in the body of this thesis, are explained in this Appendix.

§ TWO-PARTICLE CORRELATION

Two-particle correlation function is a distribution of distance of two-particle pairs in the momentum space consist of one trigger particle and another associate particle.
Definition of the two-particle correlation function introduced in Sec. <ref> and shown in Sec. <ref>
is given by

\begin{align}
C(\Delta \eta, \Delta \phi) = B(0,0) \frac{S(\Delta \eta, \Delta \phi)}{B(\Delta \eta, \Delta \phi)},
%\frac{1}{N_{\mathrm{trig}}} \frac{d^2N_{\mathrm{pair}}}{d\Delta \eta d\Delta \phi} = B(0,0) \frac{S(\Delta \eta, \Delta \phi)}{B(\Delta \eta, \Delta \phi)},
\label{eq:TwoParticleCorrelationRidge}
\end{align}

where $S(\Delta \eta, \Delta \phi)$ and $B(\Delta \eta, \Delta \phi)$ are raw two-particle correlations obtained within the same events and constructed from
randomly chosen two events, respectively.
It should be noted that there are several definitions besides Eq. eq:TwoParticleCorrelation used in the CMS collaboration.
However, the idea for those are the same:
extraction of two-particle correlation within each event contrast to backgrounds.

The raw two-particle correlations within the same events are obtained as

\begin{align}
S(\Delta \eta, \Delta \phi) =  \frac{1}{N_{\mathrm{trig}}} \frac{d^2N_{\mathrm{pair}}^{\mathrm{same}}}{d\Delta \eta d\Delta \phi},
\label{eq:Same_TwoParticle}
\end{align}

where $N_{\mathrm{trig}}$ is the number of trigger particles, and
$d^2N_{\mathrm{pair}}^{\mathrm{same}}/d\Delta \eta d\Delta \phi$ is the distribution of two-particle pair obtained within a single event as a function of $\Delta \eta$ and $\Delta \phi$.
For example, if kinematic cut is different for trigger and associate particles, $N_{\mathrm{pair}}^{\mathrm{same}}$ is $N_{\mathrm{trig}}\times N_{\mathrm{assoc}}$ where $N_{\mathrm{assoc}}$ is the number of associate particles.
The two-particle correlations constructed randomly chosen from two events are obtained as

\begin{align}
B(\Delta \eta, \Delta \phi) =  \frac{1}{N_{\mathrm{trig}}} \frac{d^2N_{\mathrm{pair}}^{\mathrm{mix}}}{d\Delta \eta d\Delta \phi}.
\label{eq:Mix_TwoParticle}
\end{align}

Here, $N_{\mathrm{pair}}^{\mathrm{mix}}$ is the number of two-particle pairs
where associate particles are randomly chosen pair-by-pair within the same centrality/multiplicity class.
Note that $N_{\mathrm{pair}}^{\mathrm{mix}}$ can differ depending on how many the dummy associate is chosen.
However, one can choose a certain proper number of dummies here
because $B(0,0)$ in Eq. eq:TwoParticleCorrelation can re-scale obtained $B(\Delta \eta, \Delta \phi)$.
In the results shown in Sec. <ref>,
20 dummy associates are picked up.
It should also be noted that $N_{\mathrm{trig}}$ is the same as in Eq. eq:Same_TwoParticle.

Top: examples of $S(\Delta \eta, \Delta \phi)$ (left) and $B(\Delta \eta, \Delta \phi)$ (right) obtained from Pb+Pb collisions. Bottom: Two-particle correlation function, $C(\Delta \eta, \Delta \phi)$, obtained from above $S(\Delta \eta, \Delta \phi)$ and $B(\Delta \eta, \Delta \phi)$.

Figure <ref> shows examples of $S(\Delta \eta, \Delta \phi)$, $B(\Delta \eta, \Delta \phi)$, and two-particle correlation function, $C(\Delta \eta, \Delta \phi)$, obtained with the $S(\Delta \eta, \Delta \phi)$ and $B(\Delta \eta, \Delta \phi)$ in Pb+Pb collisions.
One would notice the meaning of subtracting background from this.
In high-energy nuclear collisions, because particles are produced around midrapidity, there should be an obvious strong correlation in midrapidity.
In this method, by dividing the raw two-particle correlation with background,
such an obvious correlation in longitudinal direction is suppressed
to clearly see the azimuthal anisotropy.

§ SUBEVENT METHOD

The subevent method is proposed to subtract the back-to-back dijet correlation
It can be expected that the dijets tend to have short range correlation in longitudinal direction.
Under the subevent method,
one first prepares some pseudorapidity ranges which have a certain longitudinal distance for each.
Then, in one single event, one can obtain subsets of particles produced in
each of those pseudorapidity range, and the subsets of particles are called subevents.
Multi-particle correlations are calculated
by picking up particles from those subevents.

In this section,
I explain some details of the calibration of the 2-subevent method for 2-particle correlation and the 3-subevent method for 4-particle correlation
which are calculated in this thesis.
For more details, see i.e., Ref. <cit.>.

§.§.§ 2-subevent method for 2-particle correlation

In the 2-subevent method, one usually impose eta gap, $\Delta \eta$,
between the 2 subevents.
For instance, within a certain acceptance $|\eta|<\eta_{\mathrm{accept.}}$, the 2 subevents are obtained in $-\eta_{\mathrm{accept.}}<\eta<- \eta_{\mathrm{gap}}/2 $ and $ \eta_{\mathrm{gap}}/2 <\eta<\eta_{\mathrm{accept.}}$ with an eta gap $\Delta \eta>\eta_{\mathrm{gap}}$.
For the later explanation, I name these two subsets as $A$ and $B$.

Equation eq:TwoParticleCorrelation in Sec. <ref> of chapter 1 shows the definition of 2-particle correlation.
On the other hand,
the 2-particle correlation with 2-subevent method is obtained as,

\begin{align}
\label{eq:TwoParticleCorrelation_2sub}
\langle 2 \rangle = \frac{1}{N_A N_B}\sum_{i_A j_B}^{N_A,N_B} e^{in(\phi_{i_A} -\phi_{j_B} )} = \frac{Q_{n,A}Q_{n,B}^*}{N_AN_B},
\end{align}

where $N_A$ and $N_B$ are the number of particles produced in each subevent of $A$ and $B$.
On the other hands, $Q_{n,A}$ and $Q_{n,B}$ are the Q-vector defined in Eq. eq:Qvector and obtained only within a single subevent $A$ and $B$, respectively.
The change of the definition from Eq. eq:TwoParticleCorrelation to Eq. eq:TwoParticleCorrelation_2sub can be clearly understood once one thinks about its meaning:
in the 2-subevent, the number of particle pairs becomes $N_A N_B$ and the subtraction of self correlation is no longer needed because the 2-particle correlation is calculated by picking up 2-particle from each different subevent.

§.§.§ 3-subevent method for 4-particle correlation

Under the 3-subevent for 4-particle correlation, the four particles are picked up from 3 rapidity ranges dividing the acceptance range into 3, for instance, $-\eta_{\mathrm{accept.}}<\eta<\eta_{B}$, $\eta_B<\eta<\eta_C$, and $\eta_C <\eta<\eta_{\mathrm{accept.}}$,
where $\eta_B<0<\eta_C$.
From the backward (negative) to forward (positive) rapidity, I name those subevents as $B$, $A$, and $C$.

The 4-particles are picked up as follows. Two of four particles are picked up from the subevent $A$ located at the center and one particle is picked up from each $B$
and $C$ subevent, respectively.
Here, although there are two ways to calculate 4-particle correlation from 3-sub event
considering the symmetry,
I adopt the one which can reduce the effects of dijet correlations <cit.>.
The explicit form is expressed as

\begin{align}
\langle 4 \rangle &= \frac{1}{N_A(N_A-1)N_BN_C} \sum_{i_A, j_A, k_B, l_C}^{N_A, N_A, N_B, N_C} e^{in(\phi_{i_A}+\phi_{j_A}-\phi_{k_B}-\phi_{l_C})} \nonumber \\
\end{align}

For more discussions on multi-particle correlation, see Ref. <cit.>.

CHAPTER: STRING MODIFICATION ALGORITHM

In Sec. <ref> in Chapter 2, I explained the color string modification due to the spacial overlap between strings and QGP fluids.
Here, I would like to describe the algorithm used in the color string modification, which is the string tracing and string cutting.

§ STRING TRACING

As the name tells, this algorithm is used to trace color strings. Under the string tracing, I need to check if there exist any crossing points between color strings and hypersurfaces.
In Fig. <ref>, I show an example of the situation to deal with.
A color string consisting of a quark and an anti quark on a 2-dimensional plane of discretized space is shown.
The discretized space corresponds to the
hydrodynamic-simulation space.
While the actual simulation is performed
in the Milne coordinate, for space component $(x, y, \eta_s)$,
I will explain the algorithm with a general 3-D coordinate of cell number, i.e., direction 1, 2, and 3,
because the algorithm is constructed in terms of cell number rather than coordinate with physical scale.
The two leading quark and anti-quark are located at
$\vec{p}^A$ and $\vec{p}^B$, where these are vectors expressed with cell numbers.
The color string is linearly stretched in-between the leading quark and anti-quark.
It does not matter which of $\vec{p}^A$ and $\vec{p}^B$ corresponds to the quark and anti-quark,
so, hereafter, I call them just particle A and B.
Here, cells that have spacial overlap with the string are highlighted with yellow.
Thus, to look for crossing points between the color string and QGP fluids, one is required to check temperature of each yellow cell.

A color string consisting of a quark and an anti quark (red and blue) on a 2-dimensional discretized space.
Cells that have spacial overlap with the string are highlighted with yellow.

The problem we need to solve is
how we trace these cells with a reasonable computational time.
The computational cost on this algorithm matters because the number of and the length of strings generated in one heavy-ion collision event are enormous.
The answer is as follows:
one needs to trace a color string
one by one in the smallest-slope direction
of the vector connecting the particle A and B.
I will describe details in the followings.

Under this algorithm, there is one assumption which makes the calculation easy:
all partons are assumed to be located at the center of each cell and a color string stretches in between partons, which is not the case, for example in Fig. <ref>,
because all partons can be at any point within the simulation space of hydrodynamics.
However, since we take small cell widths as much as possible, the difference can be small enough to be neglected.

§.§.§ Determination of the direction to step

In the hydrodynamic simulation,
it would be a usual case that
one numbers each cell as $0, 1, 2, 3 \cdots, N_i$ where $N_i-1$ is the number of discretized cells in $i$-direction.
Under this algorithm,
I treat these cell number as the coordinate.
Here I define that the center of a cell
located at the $k_1$-th in direction 1,
$l_2$-th in direction 2,
and $m_3$-th in direction 3 corresponds to
a coordinate $(k_1, l_2, m_3)$.
In Fig. <ref>,
for instance, $\vec{p}^A$ and $\vec{p}^B$ are expressed as

\begin{align}
p^A_1 = 0.0, \enskip p^A_2 = 0.0, \enskip p^A_3 = 0.0,
\end{align}

\begin{align}
p^B_1 = 4.0, \enskip p^B_2 = 2.0, \enskip p^B_3 = 0.0.
\end{align}

First one needs to find a cell that each parton is located.
The set of cell number is obtained in the actual simulation as

\begin{align}
k_1 &= \mathrm{floor} (x/\Delta x) + N_x, \\
l_2 &= \mathrm{floor} (y/\Delta y) + N_y, \\
m_3 &= \mathrm{floor} (\eta_s/\Delta \eta_s) + N_{\eta_s}, \\
\end{align}

where $\mathrm{floor}(x)$ returns the largest integer that is not greater than $x$,
$\Delta x$, $\Delta y$, and $\Delta \eta_s$ are cell widths determined in Tab. <ref>.
Note that the origin of the coordinate $(x,y,\eta_s)$ is set at the center of the simulation space.

Here, the unit vector, $\vec{r}$, starting from the particle A and pointing the particle B is obtained as

\begin{align}
\vec{r} = \frac{\vec{p}^B-\vec{p}^A}{|\vec{p}^B-\vec{p}^A|}.
\label{eq:rAB}
\end{align}

Here, I name the direction
in which a string has the smallest slope as
$d_{\mathrm{max}}$ for the later explanation.
One can determine $d_{\mathrm{max}}$ as

\begin{align}
|r_{\mathrm{d_{\mathrm{max}}}}| = \max (|r_1|, |r_2|, |r_3|).
\label{eq:dmax}
\end{align}

In the case of Fig. <ref> where
$|r_3|=0$ is assumed, one can get

So far, there is no conditions for
$\vec{p}^A$ and $\vec{p}^B$,
but it is convenient to take the tracing direction always from small to large cell numbers.
Therefore, I define the starting point
$\vec{S}$ as

\begin{align}
\vec{S} = \left\{
\begin{aligned}
\  & \vec{p}^A \ (p_{d_{\mathrm{max}}}^A <p_{d_{\mathrm{max}}}^B)  \\
& \vec{p}^B \ (p_{d_{\mathrm{max}}}^B <p_{d_{\mathrm{max}}}^A).
\end{aligned}
\right. .
\label{eq:Svec}
\end{align}

With these variable determined above, we consider one-by-one steps in $d_{\mathrm{max}}$ direction shown in the left of Fig. <ref>.
The green arrow is the step with a width of $1$ in $d=1$ which corresponds to $d_{\mathrm{max}}$ in this case.
The right figure of Fig. <ref> shows the one-by-one step in the $d=2$ direction which is not $d_{\mathrm{max}}$.
The comparison between the left and right figures would clearly explain the reason
I take the $d_{\mathrm{max}}$ direction to see information of each cell.
For the later explanation,
I define the step with the width of $1$ in the $d_{\mathrm{max}}$ as $i_{\mathrm{cell}} = 0, 1, 2, 3 \cdots$,
then the position of the $d_{\mathrm{max}}$ direction becomes $p^A_1 + i_{\mathrm{cell}}$ after $i_{\mathrm{cell}}$ steps.

(Left) One-by-one step in the $d_{\mathrm{max}}$ direction (green arrows).
Cells that can be checked with the algorithm are highlighted with yellow.
(Right) One-by-one step in the $d\neq d_{\mathrm{max}}$ direction.
The dotted lines are grid lines in each direction.

§.§.§ Coordinates at each step

Now I need to obtain values of the other dimension on a color string.
In general expression, we can obtain the position after $i_{\mathrm{cell}}$ steps $\vec{p}^{\mathrm{now}}$  as

\begin{align}
\label{eq:pnow}
\vec{p}^{\mathrm{now}}  = \vec{S} + \frac{\vec{r}}{r_{d_{\mathrm{max}}}} (i_{\mathrm{cell}}-S_{d_{\mathrm{max}}} ).
\end{align}

This gives one the full coordinates of $\vec{p}^{\mathrm{now}}$ that is accessible with the one-by-one step illustrated in Fig. <ref>.
Basically, I trace a color string and find crossing points between hypersurfaces and color strings
by checking information at each $\vec{p}^{\mathrm{now}}$.
However, as one sees in Fig, <ref>,
there is another issue to consider:
some positions that is accessed with the one-by-one step are located at boundaries.
In order not to miss the crossing points,
I check both cells that is sharing the boundary.

Each point of $\vec{p}^{\mathrm{now}}$ that is accessible with the one-by-one step (blue crosses).
The crosses located at boundaries are assumed to span two cells sharing the boundary (filled with orange diagonal stripe pattern).

To realize that,
I need to judge where coordinates are located in a cell and which cell one should check in string tracing at every step with the following procedure.
There are three cases for steps:
1) a step is at the center of a cell,
2) a step is on the boundary of a cell, and
3) otherwise.
Coordinates are defined to take $(k_1, l_2, m_3)$ at
the center of a cell
located at the $k_1$-th in direction 1,
$l_2$-th in direction 2,
and $m_3$-th in direction 3.
The case 1) is satisfied when a value of $\vec{p}^{\mathrm{now}}$ in a direction have zero decimal part.
In this case, the coordinate exactly corresponds to the cell numbering.
The case 2) is satisfied when
a twice value of $\vec{p}^{\mathrm{now}}$ in a direction have zero decimal part.
The cell numbering to check is obtained by rounding and flooring the value.
For the case 3),
one can obtain the cell number to check by rounding the value.

More explicitly, I explain the case of Fig. <ref>.
In the direction of $d=1$, there is no ambiguity to determine cells to check since the one-by-one step is set to $d=1$.
Thus, all the coordinate of $p^{\mathrm{now}}_1$ obtained in Eq. eq:pnow
has a one-to-one correspondence with cell numbering.
For the direction of $d=2$, the $\vec{p}^{\mathrm{now}}$ at $i_{\mathrm{cell}}=0, 2, 4$ correspond to the case 1),
while that of $i_{\mathrm{cell}}=1, 3$ correspond to the case 2).
This is because in the latter case, $p^{\mathrm{now}}_2(i_{\mathrm{cell}}=1) = 1.5$ and $p^{\mathrm{now}}_2(i_{\mathrm{cell}}=3) = 2.5$.
Here, one needs to check the cells located at $p^{\mathrm{now}}_2=1$ and $p^{\mathrm{now}}_2=2$ at $i_{\mathrm{cell}}=1$,
and $p^{\mathrm{now}}_2=2$ and $p^{\mathrm{now}}_2=3$ at $i_{\mathrm{cell}}=3$ by rounding and flooring the $p^{\mathrm{now}}_2$ on a boundary.

§.§.§ Program flow

I show the program flow of the string tracing that I have explained so far in List. <ref>.
Note that this just illustrates the flow of program, not the program itself.
For a simplicity and generality, I omit accurate notation of each variables and formats of programming language
Since the actual program is written in C++, some notations are exactly expressed in C++.
I would like readers of this thesis to extract the main idea of the algorithm if one needs.

Line 8 and 9 correspond to Eqs. eq:rAB and eq:dmax.
From line 11 to 14, Eq. eq:Svec is obtained.
From line 18, the loop of one-by-one step begins.
At line 22, Eq. eq:pnow is obtained.

[language=C++, caption=Program flow of the string tracing]

//...several definitions of each function are above.

//1. Determination of the direction to step
r_vec = (p_Avec - p_Bvec)/absolute(p_Avec - p_Bvec)
r_dmax = max(r_vec[0], r_vec[1], r_vec[2])

S_vec =  p_Avec
S_vec =  p_Bvec

//One-by-one step

//2. Coordinates at each step
p_now = S_vec + (r_vec/r_dmax)*(i_cell - S_vec[dmax])

//Check location of rho_now for each direction,
//and find cell numbering to check

if (case 1)
else if (case 2)
else //case 3

//Check if there is crossing points between hypersurface and color strings.
if (crossing point)
//Archive the information!

§ CUTTING AND PAIRING

With the information obtained at crossing points between hypersurface and color strings with the string tracing method explained above,
I modify configuration of color strings as explained in Sec. <ref> in chapter 2.

In order to explain the algorithm to perform the color string modification, it should be instructive to first explain
how  handles colors (and anti-colors).
In , which partons are connected and confined as a color-singlet string is determined by assigning color and anti-color tags.
These tags are just random positive integers, not colors of ($R$, $G$, $B$) with $N_\mathrm{c}=3$.
Thus, a color and an anti-color tag having the same value can form a color singlet.
For instance, a quark having color and anti-color ($101$, $0$) and
an anti-quark with ($0$, $101$) form a color singlet string.
This treatment is valid under the leading color limit ($N_\mathrm{c} \rightarrow \infty$) <cit.>.
In , all perturbative processes are handled under this limit <cit.>, and this enables us to deal with color string modification with relative ease.

§.§.§ String cutting algorithm

Let me consider a simple case: one string that consists of a quark and anti-quark.

An example of the string cutting algorithm. A quark (red) and anti-quark (blue) is forming a string. A matter profile within hypersurface is illustrated with pale orange oval. Two crossing points (#1 and #2) between the string and the surface are expressed with yellow crosses.The process of the algorithm is shown from (a), (b), (c), to (d). The open red and blue circles are quarks and anti-quarks picked up from the crossing points.

Figure <ref> shows the flow of the string cutting algorithm performed at $\tau = \tau_s$ with a simple example.
The process starts from (a) in Fig. <ref>.
With the string tracing algorithm explained in Sec. <ref> in this Appendix B,
one can find crossing points between the hypersurface and the string as shown with the yellow crosses
by tracing the string in the direction of the green arrow.
One can also obtain the information of the fluid cell of the hypersurface at the crossing point such as thermodynamic variables and flow velocity.

Figure <ref> (b) shows how to handle the string cutting at the crossing point #1.
Since which partons are confined and form which color strings are determined with color and anti-color tags,
one can “cut” a string by modifying the tags.
To cut the string at the crossing point #1,
I add a quark and an anti-quark pair sampled from the hypersurface at the crossing point.
Here, I assign tags ($10002$,$0$) and ($0$,$10001$) for the sampled quark and anti-quark, respectively
For the color and anti-color tags of sampled partons, in order to
avoid assigning the same value as the originally exist one,
large integers ($\leq10^8$) are used in the actual simulations.
To let the program know that the leading quark ($101$, $0$) and the sampled anti-quark ($0$,$10001$) form a string,
the original color tag assigned to the leading quark is changed to $10001$.

Figure <ref> (c) shows the string cutting at #2.
Here, color and anti-color tags ($10003$,$0$) and ($0$, $10002$) are assigned to the newly sampled parton pair.
Note that the same value for color and anti-color cannot be used more than once.
The sampled anti-quark has the same anti-color value as the color tag assigned to the quark at #1.
Since this is the last crossing point appearing on this string,
the quark sampled at #2 and the original leading anti-quark with ($0$, $101$) should form a string.
This is done by replacing the anti-color tag of the original leading anti-quark with the color tag of the quark sampled at #2.

Finally Fig. <ref> (d) shows the color flow obtained just after the modification.
The middle piece of the string is discarded by identifying a string that consists only of a quark and an anti-quark assuming that such strings should be inside of the hypersurface.

Figure <ref> shows
an example of numerical treatment of the string cutting just explained above. Left figure shows the case that two strings are penetrating thermal medium having $T>T_{\mathrm{sw}}$. On the other hand, the right figure shows the modified strings due to string cutting adopted in the left hand case.
Four crossing points between hypersurface and color strings are found in this case, and
four thermal partons are picked up from the crossing points,
which cuts the original two strings into four strings.

Example of numerical treatment of string cutting. Left: Two strings (blue blocks) are penetrating thermal medium having $T>T_{\mathrm{sw}}$ (pink blocks). Right: Modified strings (green lines) due to string cutting adopted in the left hand case. Thermal partons (stars) are picked up from crossing points between hypersurface and color strings (orange blocks). Note that these are not obtained from an actual simulation. The axes $x, y,$, and $z$ are numbering of each hydro cell in $x, y$, and $\eta_s$ direction.

§.§.§ Parton pairing algorithm

The parton pairing algorithm is much simpler compared to the string cutting explained above
because one does not need to trace strings any more.
It is assumed that this kind of partons coming out from hypersurface with finite energy would form a string by picking up a single thermal parton in vacuum and are hadronized via string fragmentation.
This is realized as follows.
In Fig. <ref> (a),
a quark that comes out from a hypersurface and has a color and anti-color tag ($101$, $0$) is shown.
In order to make the string color singlet,
an anti-quark is picked up from the hypersurface in Fig. <ref> (b).
To avoid using the same value for color and anti-color tags,
new value is assigned.
The sampled anti-quark has ($0$, $10101$),
and the color tag of the quark that is about to come out from the hypersurface is changed from $101$ to $10101$.

For the case that an anti-quark comes out from a hypersurface,
the anti-quark picked up a quark from the hypersurface likewise the above case.
The case that a gluon comes out from medium is shown in Fig. <ref> (a).
The color and anti-color tags of the gluon is ($101$, $102$).
For the case of gluons, it is assumed that
a gluon picks up another gluon from a hypersurface and form a color singlet string.
In Fig. <ref> (b),
the gluon picked up from the hypersurface has color and anti-color tags, ($10102$, $10101$).
Here again, the color and anti-color tags of the gluon that comes out from the hypersurface are changed from ($101$, $102$) to ($10101$,$10102$), which form a color singlet.

The parton pairing algorithm when a quark comes out from a hypersurface.
(a) A quark with color and anti-color tag ($101$,$0$) is about to come out from the hypersurface. (b) An anti-quark with ($0$,$10101$) is sampled at the hypersurface, and the color tag of the quark is changed to $10101$.

The parton pairing algorithm when a gluon comes out from a hypersurface.
(a) A gluon with color and anti-color tag ($101$,$102$) is about to come out from the hypersurface. (b) A gluon with ($10101$,$10102$) is sampled at the hypersurface, and the color and anti-color tag of the gluon that was about to come out from the hypersurface is changed to ($10102$, $10101$).

Note that the examples shown so far are very simple cases.
In the actual simulations, for instance,
junctions or gluon loops are produced.
The former is a Y-shaped color singlet object that end points of three strings are converged, which can look like a Y shape.
The latter is a color singlet object which only consists of gluons.
To make it color singlet, gluons form a loop.
In these cases, string modifications are performed in different ways from ones explained above because of the different configurations of the color flows.

As I explained in Sec. <ref> of chapter 2, relativistic hydrodynamics is used to describe the space-time evolution of the QGP.
In this Appendix, I would like to access the details of the numerical treatment of hydrodynamic equations.

The code of (3+1)-D hydrodynamics used in this thesis is built by Tachibana in Ref. <cit.>.

In this code, the (3+1)-D hydrodynamics is solved with the Peacewise Parabolic Method (PPM).
The PPM is a robust algorithm to solve the Riemann problem with shock waves.
The first implementation to numerical simulations of relativistic heavy-ion collisions was done by Hirano in Ref. <cit.>.

Now we need to know how to solve them numerically in a discretized-cell space.
The problem one has to handle is the space-time evolution of a distribution of the density $U$ as shown in Fig.
Since we discretize space into cells in numerical simulation of hydrodynamics, $U$ is assumed to be constant in each cell.
Thus, there is a discontinuity of $U$ between two adjacent cells.
The initial value problem with such a discontinuity in numerical simulation of hydrodynamics is called, Riemann problem.

§ PPM

The numerical method for solving hydrodynamics is a huge topic as itself. Especially, one has to be careful in treatment of compressible fluids because shock wave can arise in a simulation, and which can make a simulation numerically unstable.

The Peacewise Parabolic Method (PPM) is a robust algorithm to solve the Riemann problem with shock waves. The PPM was employed to numerical simulations of relativistic heavy-ion collisions at first time in Ref. <cit.>.
This work uses the PPM in (3+1)-D hydrodynamic simulations.

§.§ Interpolation function, $U(x)$

As a first step of the PPM, one starts with interpolating density of each cell so that one can analytically obtain solutions of Riemann problem.
Suppose that the $j$th cell has density $U_j$, where $j$ is a cell index.
The $j$th cell is defined to be located within $x_{j-\frac{1}{2}}< x < x_{j+\frac{1}{2}}$ with a cell width $\Delta x = x_{j+\frac{1}{2}} - x_{j-\frac{1}{2}}$.
Here, we replace the $U_j$ with a smooth interpolation function $U_j(x)$ defined as,

\begin{align}
\label{eq:interpolated_func}
\dis U_{j} = \int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}} U_j(x) dx.
\end{align}

In the PPM, interpolation function is given as a hyperbolic form.
If we replace variable as $X=\frac{x-x_{j-\frac{1}{2}}}{\Delta x}$,  the interpolation function $U_j(x)$ can be expressed as,

\begin{align}
\label{eq:interpolation_func}
U_j(X) = aX^2 + bX + c,
\end{align}

where $a, b,$ and $c$ are unknown constants.
The boundary conditions of $U_j(X)$ given as $U_j(X=0) = U_{j,L}$ and $U_j(X=1) = U_{j,R}$ and integration of Eq. eq:interpolation_func from $X=0$ to $X=1$ gives simultaneous equations. Solving the equations, one can obtain

\begin{align}
a &= -6\left[  U_j - \frac{1}{2} \left( U_{j,R} + U_{j,L} \right)\right] \equiv U_{j,6}, \\
b &= U_{j,6} + (U_{j,R}-U_{j,L}) \equiv U_{j,6} + \Delta U_j, \\
c&= U_{j,L}.
\end{align}

Then finally the interpolation function $U_{j}(x)$ becomes,

\begin{align}
\label{eq:INTERPOLATION_Ux}
U_j(x) = - U_{j,6} \left(\frac{x-x_{j-\frac{1}{2}}}{\Delta x}\right)^2  + (U_{j,6} + \Delta U_{j} ) \frac{x-x_{j-\frac{1}{2}}}{\Delta x} + U_{j, L}.
\end{align}

§.§ Density at edge of cell, $U_{j,R}$ and $U_{j,L}$

The next problem is how we obtain the $U_{j, R}$ and $U_{j, L}$ as a function of density $U$ which is known value.
Here we analytically obtain each value of the $U_{j, R}$ and $U_{j, L}$ from 4-adjacent cells.
Suppose that there is the distribution of density $U$ shown in Fig. The cells are labeled as $j-1, j, j+1$ and $j+2$ respectively. First we assume that the integral of $U(x)$ can be parametrized as a quartic function,

\begin{align}
\label{eq:interpolation_fourcells}
\mathcal{U}(x) &= \int_{x_{j-\frac{3}{2}}}^x d\xi U(\xi)   \\
\label{eq:interpolation_fourcells_ana}
& = ax^4 + bx^3 + cx^2 + dx + e,
\end{align}

where $a,b,c,d$ and $e$ are unknown constants. Note that Eq. eq:interpolation_fourcells is defined within $x_{j-\frac{3}{2}} < x < x_{j+\frac{5}{2}}$ so the integrated interpolation function $U(\xi)$ is assumed to be defined globally within this region. Then values at cell boundaries can be calculated as,

\begin{align}
\label{eq:boundary_value}
U_{j,R} = U(x_{j+\frac{1}{2}})  = \left. \frac{d\mathcal{U}}{dx}\right|_{x=x_{j+\frac{1}{2}}}.
\end{align}

The quartic function is employed because we want to obtain $U_{j, R}$ and $U_{j, L}$ from $U$ of 4-adjacent cells.
One can solve above equations by assuming $x_{j+\frac{1}{2}}=0$ and $\mathrm{U}(x_{j-\frac{3}{2}})=0$ without loss of generality.

Since we know $U$ at each cell, integral of Eq. eq:interpolation_fourcells is easily obtained by taking sum of products of each $U$ and $\Delta x$. One would obtain the following equations:

\begin{align}
\label{eq:coupledeqA}
\mathcal{U}(x_{j-\frac{3}{2}}) &= \int_{x_{j-\frac{3}{2}}}^{x_{j-\frac{1}{2}}} d\xi U(\xi) = 0, \\
\label{eq:coupledeqB}
\mathcal{U}(x_{j-\frac{1}{2}}) &=
\int_{x_{j-\frac{3}{2}}}^{x_{j-\frac{1}{2}}} d\xi U(\xi) = U_{j-1}\Delta x, \\
\label{eq:coupledeqC}
\mathcal{U}(x_{j+\frac{1}{2}}) &=
\int_{x_{j-\frac{3}{2}}}^{x_{j+\frac{1}{2}}} d\xi U(\xi) = (U_{j-1}+U_{j})\Delta x, \\
\label{eq:coupledeqD}
\mathcal{U}(x_{j+\frac{3}{2}}) &=
\int_{x_{j-\frac{3}{2}}}^{x_{j+\frac{3}{2}}} d\xi U(\xi) = (U_{j-1}+U_{j}+U_{j+1})\Delta x,  \\
\label{eq:coupledeqE}
\mathcal{U}(x_{j+\frac{5}{2}}) &=
\int_{x_{j-\frac{3}{2}}}^{x_{j+\frac{5}{2}}} d\xi U(\xi) = (U_{j-1}+U_{j}+U_{j+1}+U_{j+2})\Delta x.
\end{align}

On the other hand, by inserting $x_{j- \frac{3}{2}}=-2\Delta x$, $x_{j-\frac{1}{2}}=-\Delta x$, $x_{j+\frac{1}{2}}=0$,  $x_{j+ \frac{3}{2}}=\Delta x$, , and $x_{j+\frac{5}{2}}=2\Delta x$ into Eq. eq:interpolation_fourcells_ana, one obtain,

\begin{align}
\label{eq:coupled_eqA}
\mathcal{U}(-2\Delta x) &= 16a \Delta x ^4 - 8b\Delta x^3 + 4c\Delta x^2 -2 d\Delta x + e, \\
\label{eq:coupled_eqB}
\mathcal{U}(-\Delta x) &=
a\Delta x^4 - b \Delta x^3 + c\Dx^2 - d \Dx + e,  \\
\label{eq:coupled_eqC}
\mathcal{U}(0) &= e, \\
\label{eq:coupled_eqD}
\mathcal{U}(\Dx) &=
a \Dx^4 + b\Dx^3 + c\Dx^2 + e\Dx,  \\
\label{eq:coupled_eqE}
\mathcal{U}(2\Delta x) &=
16 a \Dx^4 + 8b \Dx^3 + 4c \Dx^2 + 2d \Dx + e.
\end{align}

Then, Eq. eq:coupled_eqA to eq:coupled_eqE and Eq. eq:coupledeqA to eq:coupledeqE lead to 5 simultaneous equations which determine 5 unknown constants in Eq. eq:interpolation_fourcells_ana.

Finally, $a, b, c, d$, and $e$ are obtained as,

\begin{align}
a &= (-U_{j-1} + 3U_j -3U_{j+1} +U_{j+2})/24\Dx^3, \\
b &= -(-U_{j-1} + U_j +U_{j+1} -U_{j+2})/12\Dx^2, \\
c &= -(-U_{j-1} + 15U_j -15U_{j+1} +U_{j+2})/40\Dx, \\
d &= (U_j + U_{j+1})/2 + (-U_{j-1} + U_j +U_{j+1} -U_{j+2})/12, \\
e&= (U_{j-1} + U_{j}) \Dx.
\end{align}

Finally, $U_{j, R}$ becomes,

\begin{align}
\label{eq:UR}
U_{j,R} &= U(x_{j+\frac{1}{2}}) \nonumber \\
&= d \nonumber \\
&=\frac{U_j+U_{j+1}}{2} + \frac{(U_{j+1} - U_{j-1}) - (U_{j+2} - U_{j})}{12}  \nonumber \\
& = \frac{7}{12} (U_{j}+ U_{j+1}) -\frac{1}{12} (U_{j+2} + U_{j-1}).
\end{align}

On the other hand, $U_{j,L}$ becomes,

\begin{align}
\label{eq:UL}
U_{j,L} &= U(x_{j-\frac{1}{2}}) \nonumber \\
& = \frac{7}{12} (U_{j-1}+ U_{j}) -\frac{1}{12} (U_{j+1} + U_{j-2}).
\end{align}

§.§ Monotonicity of $U_R$ and $U_L$

There are some cases where monotonicity of $U_R$ and $U_L$ is lost in $x_{j-\frac{1}{2}}<x<x_{j+\frac{1}{2}}$, and which we need to fix.
(a) First case is when $U_j$, density of the $i$th cell, does not take the value between $U_{j,R}$ and $U_{j,L}$.
(b) Second case is when $U_j$ is in between $U_{j, R}$ and $U_{j,L}$ but the interpolation function $U_j(x)$ has minimum or maximum value in $x_{j-\frac{1}{2}}<x<x_{j+\frac{1}{2}}$.

To impose monotonicity on the interpolation function, we replace the values of $U_{j,L}$ and $U_{j,R}$ as follows.

\begin{align}
\label{eq:Mono_a}
U_{R,j} \rightarrow U_{j}, \ U_{L, j} \rightarrow U_{j} \nonumber \\
\mathrm{if} \  (U_{R,j}-U_{j})(U_{L,j}-U_{j})\leq0,
\end{align}

\begin{align}
\label{eq:Mono_bm}
U_{j,L} \rightarrow 3U_{j}-2U_{R,j} \nonumber \\
\mathrm{if} \ \Delta U_{j}U_{j,6} > \Delta U_{j}^2,
\end{align}

\begin{align}
\label{eq:Mono_bp}
U_{j,R} \rightarrow 3U_{j}-2U_{L,j} \nonumber \\
\mathrm{if} \ \Delta U_{j}U_{j,6} < -\Delta U_{j}^2.
\end{align}

The condition in Eq. eq:Mono_a corresponds to the case of (a) while ones in Eq. eq:Mono_bm and eq:Mono_bp correspond to the case of (b). Each situation is shown in Fig.().
For the case of (a), $U_{R,j}$ and $U_{L,j}$ are replaced so that the interpolation function becomes constant within $x_{j-\frac{1}{2}}<x<x_{j+\frac{1}{2}}$.
On the other hand, for the case of (b), $U_{R,j}$ or $U_{L,j}$ is replaced so that the interpolation function takes its maximum or minimum at the edge of the $j$th cell keeping the value of $U_j$. If $U_{R,j}$ is closer to the maximum or minimum point $U_{L,j}$ is moved while $U_{R,j}$ and $U_{j}$ is fixed, and vice versa.

Just in case, note that the conditions in Eq. eq:Mono_bm and eq:Mono_bp are obtained as follows.
Remember that the interpolation function is expressed as Eq. eq:INTERPOLATION_Ux. By replacing $\frac{x-x_{j-\frac{1}{2}}}{\Dx}$ with $X$ and taking derivative with $X$, one would obtain

\begin{align}
\frac{dU_j(X)}{dX} = -2U_{6,j}X+U_{6,j}+\Delta U_j, \ (0\leq X\leq 1).
\end{align}

This leads to the following conditions of non-monotonicity,

\begin{align}
& 0\leq\frac{U_{6,j}+\Delta U_j}{2U_{6,j}}\leq 1 \nonumber \\
\Leftrightarrow \  & -U_{6,j}\leq\Delta U_j \leq U_{6,j} \nonumber \\
\Leftrightarrow \  & |\Delta U_j |>| U_{6,j}|.
\end{align}

The conditions in Eq. eq:Mono_bm and Eq. eq:Mono_bp corresponds to the cases where $\Delta U_j<0$ and $0<\Delta U_j$, respectively.

§.§ Contribution of propagating $U_R$ and $U_L$ in $\Delta t$

Now we have the interpolation function of Eq. eq:INTERPOLATION_Ux as a function of density $U$ of each cell at a fixed time.
However, in the form of Eq. eq:UR and eq:UL, there
How fluids propagate can be evaluated by using a signal velocity. During $\Delta t$, discontinuity flows from the right and left into the constant region, as shown in Fig.().
Average values of interpolation functions $U(x)$ that flow into the boundary is obtained as follows.

\begin{align}
\label{eq:ave_UL}
\bar{U}_{j+1,L}& = \frac{1}{|b_{j+\frac{1}{2},r}|\Delta t} \int_{x_{j+\frac{1}{2}}}^{x_{j+\frac{1}{2}} + |b_{j+\frac{1}{2}, r}|\Delta t} U_{j+1} (x) dx, \\
\label{eq:ave_UR}
\bar{U}_{j,R} &= \frac{1}{|b_{j,l}|\Delta t} \int^{x_{j+\frac{1}{2}}}_{x_{j+\frac{1}{2}} - |b_{j+\frac{1}{2}, r}|\Delta t} U_{j} (x) dx
\end{align}

\begin{align}
b = \frac{v+c}{1+vc}
\end{align}

The Eq. eq:ave_UL becomes,

\begin{align}
\bar{U}_{L, j+1} &= \frac{1}{|b|\Delta t } \int_{x_{j+\frac{1}{2}}}^{x_{j+\frac{1}{2}} + |b|\Delta t} dx U_{j+1} (x) \nonumber \\
& = \frac{1}{|b|\Delta t} \int_{x_{j+\frac{1}{2}}}^{x_{j+\frac{1}{2}} + |b|\Delta t} dx U_{L, j+1} + \frac{x-x_{j+\frac{1}{2}}}{\Dx}\left[ \Delta U_{j+1} + U_{6, j+1} \left(1-\frac{x-x_{j+\frac{1}{2}} }{\Dx}\right) \right] \nonumber \\
& = U_{L, j+1} + \frac{1}{|b|\Delta t} \left[
\int_{x_{j+\frac{1}{2}}}^{x_{j+\frac{1}{2}} + |b|\Delta t} dx \frac{x-x_{j+\frac{1}{2}}}{\Dx} \Delta U_{j+1} \right. \nonumber \\
& \hspace{30mm} \left. + \int_{x_{j+\frac{1}{2}}}^{x_{j+\frac{1}{2}} + |b|\Delta t} dx \left[ \frac{x-x_{j+\frac{1}{2}}}{\Delta x} U_{6, j+1} - \left( \frac{x-x_{j+\frac{1}{2}}}{\Dx}\right)^2 U_{6, j+1} \right] \right] \nonumber.
\end{align}

The first term in the integral can be computed as,

\begin{align}
\int_{x_{j+\frac{1}{2}}}^{x_{j+\frac{1}{2}} + |b|\Delta t} dx  \frac{x-x_{j+\frac{1}{2}}}{\Dx} \Delta U_{j+1}
&= \int_{x_{j+\frac{1}{2}}}^{x_{j+\frac{1}{2}} + |b|\Delta t} dx \frac{\Delta U}{\Delta x} - \int_{x_{j+\frac{1}{2}}}^{x_{j+\frac{1}{2}} + |b|\Delta t} dx \frac{\Delta U}{\Dx} x_{j+\frac{1}{2}} \nonumber\\
& = \left[ \frac{1}{2} \left[(x_{j+\frac{1}{2}} + |b|\Delta t)^2 - x_{j+\frac{1}{2}}^2 \right] - (|b|\Delta t x_{j+\frac{1}{2}}) \right] \frac{\Delta U}{\Dx} \nonumber \\
& = \frac{1}{2} |b|^2 \Delta t^2 \frac{\Delta U}{\Delta x}
\end{align}

The second term becomes,

\begin{align}
\dis & \int_{x_{j+\frac{1}{2}}}^{x_{j+\frac{1}{2}} + |b|\Delta t} dx \left[ \frac{x-x_{j+\frac{1}{2}}}{\Delta x} U_{6, j+1} - \left( \frac{x-x_{j+\frac{1}{2}}}{\Dx}\right)^2 U_{6, j+1} \right] \nonumber \\
&= \frac{U_6}{|b|\Delta t} \left[\frac{1}{2\Dx} |b|^2 \Delta t^2  \right] - \frac{U_6}{|b|\Delta t \Dx^2} \left[\frac{1}{3} |b|^3 \Delta t^3\right] \nonumber \\
& = \frac{U_6}{2\Dx}|b|\Delta t \left[ 1- \frac{2|b|\Delta t}{3\Dx}\right].
\end{align}

Finally, one would obtain

\begin{align}
\bar{U}_{L,j+1}=U_{L,j+1} + \frac{b\Delta t }{2\Dx} \left[\Delta U + \left(1-\frac{2b\Delta t}{3\Dx}U_{6} \right)\right].
\end{align}

In the same way, eq:UR is expressed as

\begin{align}
\bar{U}_R = U_R - \frac{b\Delta t }{2\Dx} \left[\Delta U - \left( 1 - \frac{2b\Delta t}{3\Delta x}U_6 \right)\right].
\end{align}

§.§.§ (3+1)-D relativistic hydrodynamics

* Reference:

* Language:

The code of (3+1)-D hydrodynamics used in this thesis is
built by Tachibana in Ref. <cit.>,
therefore I leave the detailed explanation for the reference. The (3+1)-D hydrodynamics is solved with the Piecewise Parabolic Method (PPM).
The PPM is a robust algorithm to solve the Riemann problem with shock waves. The first implementation to numerical simulations of relativistic heavy-ion collisions was done by Hirano in Ref. <cit.>.

§ STATISTICS

§.§ Standard errors

When it comes to deal with event-by-event fluctuating data, it must be inevitable to get familiar with statistical analysis.

Although it would be quite straight-forward to get average and error-bar if you have simple event-by-event data set, that is not always the case.
For instance, if one uses Monte-Carlo event generator it is sometimes necessary to generate weighted events to reduce statistical errors. Suppose that one needs to analyze centrality dependence of certain observable. However, because a cross-section of colliding nuclei is proportional to an impact parameter, i.e., $d\sigma(b)/db \propto b$, it would be hard to gain statistics of collisions with small impact-parameter. One feasible way to solve this issue is to generate weighted events by sampling impact parameter according to some probability distributions, for instance, Gaussian function centered at zero. An inverse of the factor of Gaussian function corresponds to a weight of the event. This enables us to obtain enough statistics of collisions with small impact parameters.
On the other hand, in order to get event-average, one needs to take into account of the weights in statistical analysis.

In this section, I summarize how to calculate standard errors in weighted events. Definitions of standard errors for two different weighting methods are explained below. First, error bars in weighted events where each event is able to be counted as one complete event. Second, those where each event is $\it{not}$ able to be counted as one complete event.

Before going to a discussion of weighted event, I define average $\bar{x}$, unbiased sample variance $\sigma^2$ (here after I call this variance for simplicity), and corresponding standard error $\mathrm{SE}$ in un-weighted events.

\begin{equation}
\begin{aligned}
\label{eq:statval_unweighted}
\bar{x} &=& \frac{\dis\sum_i^\Nev x_i}{\dis\sum_j^\Nev 1}, \hspace{5mm}
\sigma^2 &=& \frac{\dis\sum_i^\Nev (x_i - \bar{x})^2}{\dis\sum_j^\Nev 1}, \hspace{5mm}
\mathrm{SE} &=& \sqrt{\frac{\sigma^2}{\Nev}},
\end{aligned}
\end{equation}

where $x_i$ is a observable from a $i$th event and $\Nev$ is the number of event that one wants to analyze.
It is often the case with that there is not so much problem if one uses unbiased sample variance instead of biased one as an error bar in plots, since we often perform a number of event such as $\Nev >> 1$.

If the events are weighted and each produced event is able to be counted as one complete event, each definition in Eq. eq:statval_unweighted becomes as follows.

\begin{equation}
\begin{aligned}
\label{eq:statval_weighted_normal}
\bar{x} &=& \frac{\dis\sum_i^\Nev x_i w_i}{\dis\sum_j^\Nev w_j}, \hspace{5mm}
\sigma^2 &=& \frac{\dis\sum_i^\Nev w_i (x_i - \bar{x})^2}{\dis\sum_j^\Nev w_j}, \hspace{5mm}
\mathrm{SE} &=& \sqrt{\frac{\sigma^2}{\Nev}},
\end{aligned}
\end{equation}

where $w_i$ is a weight of the $i$th event. For the variance, it is convenient to write down the value only with event-by-event variables to avoid adding an extra for-loop obtaining average in an analysis code.
Inserting the definition of average $\bar{x}$, one can readily obtain

\begin{align}
\sigma^2 = \frac{\dis\sum_i^\Nev w_i x_i^2 \sum_k^\Nev w_k
-  \LSq \sum_i^\Nev x_iw_i \RSq^2}
{ \dis \LSq \sum_j^\Nev w_j  \RSq^2}.
\end{align}

Next I discuss the case if one generated event is not able to be counted as one complete event.
In this case, we have to re-define each value since definition of “event-average” is different from the above one.

To make a discussion clear, let us think about one specific example.
When we try to get event-average of $\pT$ spectra using event-generator, it is often the case that it is hard to get enough statistics of events with high $\pT$ particle productions since a cross-section of such a high $\pT$ event is small.
To resolve this difficulty, one can specify transverse momentum produced in a hard process such as 2$\rightarrow$2, which is often denoted as $\hat{p}_T$, to generate high $\pT$ events. However, in order to get event-average of $\pT$ spectra, since a cross-section of such an event should be smaller than a total cross-section, a weight must be assigned to each data set generated with a certain range of $\hat{p}_T$. A weight is a fraction of an average of cross-sections of a certain range of $\hat{p}_T$ event to a total cross-section, $w_i = \sigma_{i}/\sigma_{\mathrm{tot}}$ [Ideally speaking, standard errors for the average of cross-section should be considered. In this discussion I just neglect them for simplicity.].
Suppose that one generates $N_\pThat$ event for $n$ each range of $\pThat$, and the number of total generated event is $\Nev$, i.e., $\Nev = \sum_i^n N_\pThat = n N_\pThat$. Note that the weight is normalized as $\sum_i^n w_i = \sum_i^n \sigma_i/\sigma_{\mathrm{tot}} = 1$.
In this case, to get straight to the point, each value in Eq. eq:statval_unweighted becomes

\begin{equation}
\begin{aligned}
\label{eq:statval_weighted}
\bar{x} &=& \frac{\dis\sum_i^\Nev x_i w_i}{\dis\sum_j^\Nev w_j}, \hspace{5mm}
\sigma^2 &=& \frac{\dis \sum_i^\Nev  \LSq x_i w_i \RSq^2}{\dis\sum_j^\Nev w_j} - \LSq \frac{\dis \sum_i^\Nev x_i w_i}{\dis \sum_j^\Nev w_j} \RSq^2, \hspace{5mm}
\mathrm{SE} &=& \sqrt{\frac{\sigma^2}{\dis \sum_i^\Nev w_i}}.
\end{aligned}
\end{equation}

These can be understood when one back to the definition of event-average.
In the above case, one event is composed of each range of $\pThat$ and each data set of $\pThat$ has a corresponding weight.
This means that an observable value of one complete event corresponds to $x_i^\prime \defeq \sum_j^n x_{ij} w_j$ where $x_{ij}$ is a observable value of the $i$th event generated with the $j$th range of $\pThat$. This $x_i^\prime$ is related to $x_i$ in Eq. eq:statval_weighted as

\begin{align}
\dis \sum_i^\Nev x_i w_i = \sum_i^{N_\pThat}\sum_j^n x_{ij} w_j =  \sum_i^{N_\pThat} x_i^\prime.
\end{align}

It would get clear if one rewrite Eq. eq:statval_weighted with these definition.
Remembering that $\sum_i^\Nev w_i = N_\pThat$, average, variance and standard error become

\begin{align}
\label{eq:average_weighted_clear}
\bar{x} = \frac{\dis \sum_i^{N_{\pThat}} x_i^\prime}{N_\pThat},
\end{align}

\begin{equation}
\begin{aligned}
\label{eq:variance_weighted_clear}
\sigma^2 &=& \frac{\dis \sum_i^{N_{\pThat}} \sum_j^{n}  \LSq x_{ij} w_i \RSq^2}{N_\pThat}
- \LSq \frac{ \sum_i^{N_\pThat} x_i^\prime }{N_\pThat} \RSq^2 \\
&\approx& \frac{\dis \sum_i^{N_{\pThat}} {x_i^\prime}^2 }{N_\pThat}
- \LSq \frac{ \sum_i^{N_\pThat} x_i^\prime }{N_\pThat} \RSq^2,
\end{aligned}
\end{equation}

\begin{align}
\label{eq:SE_weighted_clear}
\mathrm{SE} = \sqrt{\frac{\sigma^2}{N_\pThat}}.
\end{align}

From first to second line in Eq. eq:variance_weighted_clear, an approximation, $\sum_j^{n}  \LSq x_{ij} w_i \RSq^2 \approx {x_i^\prime}^2$, is used since it is a usual case that values of weights satisfy $w_0 \sim 1>> w_1 >> w_2 >> \cdots >> w_n $
when a corresponding average of $\pThat$ is $ {\hat{p}_{T0}} < {\hat{p}_{T1}} < {\hat{p}_{T2}} < \cdots < {\hat{p}_{Tn}} $.
The meaning of Eq. eq:average_weighted_clear, eq:variance_weighted_clear, and eq:SE_weighted_clear would be clear now. These are having forms of average, variance, and standard error for $N_\pThat$ events.

Finally I mention some applicability of each definition of standard error.
As I already mentioned at the beginning of this section, one has to use these two different definitions depending on whether one generated event is considered as one complete event or not.
Besides that, notice that weights should be normalized if one uses Eq. eq:statval_weighted while that is not necessary to be done for Eq. eq:statval_weighted_normal.