Yang Wang

0000-0002-1512-5653]Yang Wang (汪洋)

0000-0002-5823-0349]Chenxing Dong (董辰兴)

Hengxin Ruan (阮恒心)

Qiufan Lin (林秋帆)

0000-0002-9300-2632]Yucheng Zhang (张宇澄)

Shupei Chen (陈树沛)

The star formation history (SFH) is a key issue in the evolution of galaxies.
In this work,
we developed a model based on a Gaussian and gamma function mixture to fit SFHs with varying numbers of components.
Our primary objective was to use this model to reveal the shape of SFHs and the corresponding physical driving factors.
Specifically, we applied this model to fit SFHs from the TNG100-1 simulation.
Our study led to the following findings: 1) Our model fits with TNG star formation histories well, especially for high-mass and red galaxies; 2) A clear relationship exists between the number and shape of fitted components and the mass and color of galaxies, with notable differences observed between central/isolated and satellite galaxies. 3) Our model allowed us to extract different episodes of star formation within star formation histories with ease and analyze the duration and timing of each star formation episode. Our findings indicated a strong relationship between the timing of each star formation episode and galaxy mass and color.

§ INTRODUCTION

Star formation is one of the most important processes driving the evolution of galaxies.
Well measured star formation history (SFH) can help us understand many aspects in galaxy evolution,
such as the underlying physics of boosting and quenching of star formation,
and the timing and time scale of these mechanisms
[Thomas et al., 2005, Gallazzi et al., 2005, Panter et al., 2007, Choi et al., 2014, Conroy et al., 2014, Pacifici et al., 2016, Carnall et al., 2018, Schreiber et al., 2018].
On the other hand, well modeled SFH is important when inferring galaxy properties, such as stellar mass, SFR, metallicity, and dust contents, from the observed spectra or spectral energy distribution (SED)
[Conroy, 2013, Kauffmann, 2014, Janowiecki et al., 2017, Leja et al., 2017, Telles & Melnick, 2018, Zhou et al., 2020].

Currently, two commonly used methods exist for modeling the SFH: parametric and nonparametric.
Parametric models employ simple analytical formulas with a small number of parameters to describe the SFH.
Despite the inherent complexity of SFHs in real galaxies, these simple functions demonstrate their ability to derive the SFHs and constrain galaxy properties [Carnall et al., 2019].
Researchers explore various forms of parameterized formulas to enhance predictive accuracy, including exponential decline ($\tau$ model) [Reddy et al., 2012, Wuyts et al., 2011], delayed exponential decline [Gavazzi et al., 2002, Behroozi et al., 2010, Lee et al., 2010], rising form [Buat et al., 2008, Maraston et al., 2010, Papovich et al., 2011], log-normal [Gladders et al., 2013, Abramson et al., 2015, Diemer et al., 2017], double power laws [Carnall et al., 2018, Pacifici et al., 2016], modified exponentially declining [Simha et al., 2014, Ciesla et al., 2016], $\Gamma$ function [Lu et al., 2015, Zhou et al., 2020, Lu et al., 2016], Gaussian function [Bellstedt et al., 2020], and Bessel function [Iyer & Gawiser, 2017].
Parametric methods offer the advantages of computational efficiency, conceptual simplicity, and strong physical motivation.

Parametric models, however, cannot capture all SFHs and may introduce bias if inappropriate functions are utilized [Simha et al., 2014, Carnall et al., 2019].
Consequently, researchers have introduced nonparametric models to accommodate the increased complexity and variability of SFH shapes.
The simplest nonparametric model employs a piecewise constant function to fit the $SFR(t)$ [Ocvirk et al., 2006, Kelson, 2014, Leja et al., 2017, Chauke et al., 2018].
The piecewise function method can be enhanced by incorporating adaptive time binning [Tojeiro et al., 2007].
Other methods include direct fitting of SFH patterns from theoretical models [Finlator et al., 2007, Pacifici et al., 2012], polynomial expansion of the SFH [Jiménez-López et al., 2022], and Gaussian process regression [Iyer et al., 2019].
Although nonparametric methods offer improved flexibility and accuracy,
they are computationally demanding and are susceptible to various degeneracies due to the large number of involved free variables [Walcher et al., 2011, Conroy, 2013, Lower et al., 2020].

Regardless of the methods employed, the shape of SFHs is always a significant consideration.
There hasn't been universal agreement on the best form(s) of SFH models.
Furthermore, [Carnall et al., 2019] and [Leja et al., 2019] have shown that the form of prior significantly impacts the posterior distributions of SFHs.
Thus, selecting an appropriate prior, which is based on our comprehension of SFH shapes, is crucial in accurately reconstructing SFHs.

With the emergency of hydrodynamical cosmological simulations (e.g., Illustris[Vogelsberger et al., 2014], IllustirsTNG[Nelson et al., 2019, Pillepich et al., 2018], EAGLE[Schaye et al., 2015, Crain et al., 2015], FIRE[Wetzel et al., 2016, Hopkins et al., 2018]), it is possible to get numerous realistic and physics motivated SFHs.
These simulated SFHs provide  references and training sets for studies aimed at recovering SFHs <cit.>.
Analyzing simulated SFHs will assist in selecting appropriate priors for observational tasks.
Moreover, analyzing the SFHs is necessary to infer the underlying physics of galaxy formation.
The simulations encompass numerous processes that play significant roles in galaxy SFHs. These processes include environment-dependent accretion, stochastic variations, minor and major mergers, gas-to-star conversion dictated by physical conditions in the interstellar medium, ejection of gas through galactic winds, and subsequent recycling of ejected material through accretion.
Explicit functions can provide insights into the roles and quantities of these physical processes in galaxy formation.
For these reasons, the parameterization of SFHs derived from simulations is a topic worthy of investigation <cit.>.

In this study, our approach involves fitting the SFHs from using a combination of two basis forms: Gaussian distribution and $\Gamma$ distribution.
Our idea aligns closely with the dense basis approach advocated by [Iyer & Gawiser, 2017], with the distinction that we utilize a reduced number of basis forms.
Our goal is to uncover the underlying patterns of SFHs and offer valuable insights for studies focused on recovering SFHs.

The paper is organized as following.
In Section <ref>, we present the dataset employed in this study.
In Section <ref>, we analyze the performance of fitting SFHs.
In Section <ref>, we classify the SFHs based on our fitting results and investigate the relationship between SFH types and galaxy properties.
In Section <ref>, we examine the specific features of SFHs, including peak positions, peak widths, and peak separations. Additionally, we analyze the relationship between SFH characteristics and galaxy properties.
Finally, we present our conclusions and engage in additional discussions in Section <ref>.

§ DATA AND METHODOLOGY

This section provides a brief introduction to the data used in our study.
We explore the star formation history of simulated galaxies using simulation.

§.§ The Simulation and Samples

To conduct our analysis, we select galaxies from the simulation and reconstruct their SFHs.
The simulation is a large-scale, cosmological simulation which uses the AREPO moving mesh code developed by [Springel, 2010] and includes the effects of gravity and magnetohydrodynamics. Its cosmological parameters, in agreement with Planck2015 [Collaboration et al., 2016], are $\Omega_{\Lambda} = 0.6911$, $\Omega_m = 0.3089$, $\Omega_b = 0.0486$, $\sigma_8 = 0.8159$, $n_s = 0.9667$, and $h = 0.6774$. The simulation covers a box size of $110.7^3 \Mpc^3$ with $1820^3$ dark matter particles and $1820^3$ initial hydrodynamic cells. The dark matter particle mass resolution is $7.5 \times 10^6 M_{\odot}$ and the initial mass resolution of baryons is $1.4 \times 10^6 M_{\odot}$. The simulation has 100 output snapshots, covering the redshift range 127 to 0. The TNG100-1 simulation features an updated physical model that encompasses the revised recipes for star formation and evolution, chemical enrichment, cooling, and feedbacks [Weinberger et al., 2017, Pillepich et al., 2018, Pillepich et al., 2018, Nelson et al., 2018]. It also features a modified AGN feedback model to regulate massive galaxies [Weinberger et al., 2017] and a galactic winds model to shape low-mass galaxies [Pillepich et al., 2018]. More details about the simulation and its data release can be found in the introductory paper series of TNG100-1 [Pillepich et al., 2018, Springel et al., 2018, Nelson et al., 2018, Naiman et al., 2018, Marinacci et al., 2018] and in the data release itself [Nelson et al., 2019].

A total of $43,440$ galaxies were identified using the SUBFIND algorithm [Springel et al., 2001].
These galaxies satisfy two criteria, with $M_*(z=0)>10^8M_{\odot}/h$ at redshift $0$, and a merger tree length constraint such that it is longer than 10 snapshots.
We obtained the merger trees of these galaxies from the TNG project supplementary file, which are generated using the SubLink algorithm [Rodriguez-Gomez et al., 2015].

The SFH of a galaxy is determined based on its stellar population.
Each stellar particle records its age, measured in terms of lookback time $t_L$ in units of giga-years (Gyr), as well as its initial stellar mass.
We divide their initial stellar masses into different time bins based on their age, using a time bin width of $\Delta t = 0.1\ {\rm Gyr}$.
The SFH of each galaxy is then calculated as $\SFR(t_L)=(M_{*,initial}(t_L-\Delta t) - M_{*,initial}(t_L))/\Delta t$.
This definition of SFH closely resembles those reconstructed through observations, with numerous studies deriving SFR estimates using identical methods [Donnari et al., 2019, Hahn et al., 2019, Matthee & Schaye, 2019, Conroy, 2013, Johnson et al., 2013].
The time bin width of $0.1\ {\rm Gyr}$ is widely employed in previous works.
The chosen of time bin may affect the reconstruction of SFH.
We represent the results of tests on different time bin widths in Appendix <ref>.
We found that the time bin width of $0.1\ {\rm Gyr}$ is appropriate for the SFH reconstruction based on the TNG100-1 data.
Different from previous studies [Benítez-Llambay et al., 2015, Digby et al., 2019, Joshi et al., 2021], we employ the initial stellar mass rather than the stellar mass at $z=0$ to reconstruct the SFH.
In the simulation, the stellar particles return approximately $40\%$ of in their mass to inter stellar medium during their evolution [Pillepich et al., 2018] .
Hence, utilizing the initial stellar mass allows us to account for stellar mass loss and obtain a more accurate estimation of the SFH.
Appendix <ref> shows the difference between the SFH derived from the initial stellar mass and the SFH derived from the stellar mass at $z=0$.

Four examples of the fittings to the SFHs.
The blue solid line represents the actual SFH data, and the orange dashed line represents the corresponding fitted curve.
The dotted lines show the decomposition of components in the fitted curve.
The goodness of fit for the overall fitting is displayed in the legend.
We illustrates the situations for different level of fitting goodness: good fit containing multiple components (top left), good fit containing one component (top right), median fit (bottom left), and bad fit (bottom right).

§.§ Fitting Method

We employed a range of combinations of Gaussian and Gamma distributions to establish the best-fit to the SFH curves.
The reason we use Gaussian and Gamma distributions as basic functions is because these two patterns have been widely applied in previous works on SFH reconstruction and spectral energy distribution (SED) fitting [Iyer & Gawiser, 2017, Zhou et al., 2020].
Moreover, the shapes of these two distributions can respectively fit symmetric and asymmetric distributions, thereby enabling the overall fitting function to have better universality while limiting the number of basis.
The general fitting formula was expressed as follows:

\begin{equation}
SFR(t_L) = \sum_{i=1}^{N_c} \varphi_i(t_L),
\end{equation}

where $\varphi_i(t_L)$ could be either a Gaussian distribution,

\begin{equation}
\varphi_i(t_L) = C_i\frac{e^{-(t_L-\mu_i)^2/2\sigma_i^2}}{\sqrt{2\pi}\sigma_i},
\end{equation}

or a Gamma distribution,

\begin{equation}
\varphi_i(t_L) = C_i\frac{\beta_i^{\alpha_i} (t_L-\nu_i)^{\alpha_i-1} e^{-\beta_i(t_L-\nu_i)}}{\Gamma(\alpha_i)},
\end{equation}

Here, $C_i$, $\mu_i$, $\sigma_i$, $\alpha_i$, $\beta_i$, and $\nu_i$ are fitting parameters.

Also, we utilized the $R^2$ to evaluate the goodness of our fitting through all steps.
For a sequence $y_i$ and its corresponding fitting values $y_{fit,i}$, the $R^2$ is calculated as:

\begin{equation}
R^2=1-\frac{\sum_i(y_{fit,i}-y_i )^2}{n-p-1} / \frac{\sum_i(y_i - \overline{y})^2}{n-1},
\end{equation}

where $\overline{y}$ represents the mean of the sequence $y_i$,
$n$ is the number of data points involved in the fitting,
and $p$ is the number of variables, which is set to $1$ in this case.
The $R^2$ ranges between $-\infty$ and $1$, with the results closer to $1$ indicating better agreement between data and model.
While generally ranging from $0$ to $1$, $R^2$ may be negative if the fitting is severely poor.

The fitting process is given in pseudo-code Algorithm <ref>.
The main idea is to obtain the best fit of each component independently first and then combine them together.
Each component is adjusted to match either the local or global curve of the SFH, and both Gaussian and gamma distribution shape of fittings are tested. The best component fit is determined by the maximum $R^2$ value.
The overall fit is obtained by summing all the components together.
Overall fits with different numbers of maximum components $N_c$, ranging from $1$ to $6$, are created.
Then, the best fit is found based on the $R^2$ values of overall fit candidates.
We did not test the SFH fittings with $N_c$ greater than $6$ due to the very low fraction of these samples.
Our fitting results indicate a significant decrease in the number of SFHs with an increasing number of components.
In current set, fraction of SFH fits with $6$ components is $4\%$.
The fraction of SFH fits with more than $6$ components should make up less than this value.

$x=t_L$, $y=SFR(t_L)$
parameters of fiting $\hat{y}$
$N_c \in [1, 6]$
Try different maximum number of components $N_c$
$i \in [1, N_c]$
Fitting the $i$th component
find the maximum point of $y(x)$ as $(x_i, y_i)$
select local data $y_{local} = y( x_i - 1 Gyr < x \le x_i + 1 Gyr) $
fit $y_{local}$ with Gaussian distribution, get $\hat{y}_{local,gau}$
fit $y_{local}$ with Gamma distribution, get $\hat{y}_{local,gam}$
select global data $y_{global} = y(all range) $
fit $y_{global}$ with Gaussian distribution, get $\hat{y}_{global,gau}$
fit $y_{global}$ with Gamma distribution, get $\hat{y}_{global,gam}$
the priors for Gaussian fitting are $\mu_i = x_i$, $\sigma_i = 1$
the ranges of parameters are $\mu_i \in [0, 14]$, $\sigma_i \in (0,10)$
the priors for Gamma fitting are $\nu_i = x_i$, $\beta_i = 1$, $\alpha_i=1$
the ranges of parameters are $\nu_i \in [0, 14]$, $\beta_i \in [0.1,+\infty], \alpha_i \in [0,100]$
calculate $R^2$ of $y_{local,gau}$, $y_{local,gam}$,$y_{global,gau}$ and $y_{global,gam}$
find best component fit $\hat{y}_i = \arg\max\limits_{R^2}(y_{local,gau}, y_{local,gam},y_{global,gau}, y_{global,gam}$)

remove fitted component $y = y - \hat{y}_i$
change negative point to zero $y[y<0] = 0$

combine all components  $\hat{y}'_{N_c} = \sum_{i}^{N_c}{\hat{y}_i}$
use $\hat{y}'_{N_c}$ as prior to redo the fitting to $SFR(t_L)$ to get tuned $\hat{y}''_{N_c}$
calculate $R^2$ of $\hat{y}''_{N_c}$

select the best fitting $\hat{y} = \arg \max \limits_{R^2}(\hat{y}''_{1}, \hat{y}''_{2}, ... , \hat{y}''_{6})$
Fitting algorithm for one SFH

FigFitAll illustrates the fittings to four galaxies SFHs.
We select these four SFHs to cover a range of scenarios with regard to different fitting components and goodness.
As it shows, our method is capable of accurately capturing the trends of SFHs, even when the $R^2$ value is low.
Furthermore, the bottom right sub-plot reveals that the extremely small $R^2$ value is attributed to the significant fluctuations of the origin SFH.
By applying the F-test to determine significance, we find that only $281$ SFH fittings fail to reach the $5\%$ significance level (i.e., $p$ value is greater than $0.05$),
indicating that only these $281$ fittings should be considered unreliable.
The $R^2$ values of thess unreliable fittings are from $-0.067$ to $0.021$.
However, in order to improve the credibility of following analysis, we further define a group of well-matched samples by selecting SFH fittings that satisfy an arbitrary criterion of $R^2 \ge 0.5$.
There are $23,229$ SFHs remained in this group.
In the following sections, we will specify whether we use all samples or well-matched samples for our investigations.

The SFHs of galaxies are classified based on the properties of the fittings.
Six categories labeled as “C1” to “C6” are assembled based on the number of fitting components.
The SFHs are also divided into three additional categories based on the major shapes of the fitting components.
These categories are “G” type SFH with more Gaussian components, “$\Gamma$” type SFH with more Gamma components, and “G$\Gamma$” type SFH with an equal number of both types of components.
TabFitNumR presents the number of SFHs of each type.

As TabFitNumR shows,
fittings with fewer components are more prevalent than those with multiple components.
This is because the number of galaxies increases exponentially as the mass decreases, and smaller mass galaxies dominate in terms of quantity.
The star formation process in these smaller mass galaxies is usually shorter, therefore it can generally be described with fewer components.
On the other hand, there are also numerical reasons for this, such as our method tending to achieve fitting with the fewest possible components.
TabFitNumR also indicates that SFHs that match our fitting model well contain somewhat more Gaussian components than Gamma components, with the exception of single-component galaxies which exhibit a tendency towards a Gamma-shaped distribution.

The numbers of SFHs with different fitting components. The three middle columns show the statistics for all SFHs, while the rightmost columns show the results for SFHs that are well fitted with the Gaussian+Gamma format.
Types “C1” to “C6” mean the SFH types classified by component number from $1$ to $6$.
Types “$\Gamma$”, “G$\Gamma$” and “G” mean the SFH types classified by the major component shape, corresponding to SFHs with more Gamma components, with an equal number of two components, or with more Gaussian components.

4c|all samples     4c$R^2 \ge 0.5$ samples

types                                                 $\Gamma$                                      G$\Gamma$     G                   $\Gamma$     G$\Gamma$     G

C1               19462                                10290                                         0             9172      8162      4733         0             3429

C2               10897                                2933                                          5140          2824      5786      1259         2626          1901

C3               5719                                 2618                                          0             3101      3680      1358         0             2322

C4               3450                                 888                                           1181          1381      2446      462          856           1128

C5               2121                                 822                                           0             1299      1671      596          0             1075

C6               1791                                 449                                           533           809       1484      346          443           695

total number     43440                                18000                                         6854          18586     23229     8754         3925          10550

§ GOODNESS OF FIT

Histograms of $R^2$ of SFH fittings.
The overall distribution is shown as a gray histogram.
The histograms of $R^2$ values for SFHs of different types of fittings are shown in different colors.
The SFH fittings are classified by the number of components (left plot) or by the shape (Gamma or Gaussian distribution) of  dominant component(right plot).

The two-dimensional distribution of $R^2$ of SFH fittings and stellar mass of corresponding galaxies at $z=0$.
The gray line displays the median value of $R^2$ for galaxies of different masses, while the error bars represent the range of $R^2$ values from $20\%$ to $80\%$ in each mass bin.

The two-dimensional distribution of $R^2$ of SFH fittings and color of corresponding galaxies at $z=0$.
The gray line displays the median value of $R^2$ for galaxies of different colors, while the error bars represent the range of $R^2$ values from $20\%$ to $80\%$ in each color bin.

The two-dimensional distribution of $R^2$ of SFH fittings and the aga of them.
The gray line displays the median value of $R^2$ for galaxies of different SFH lengths, while the error bars represent the range of $R^2$ values from $20\%$ to $80\%$ in each age bin.

The cumulative density function (CDF) of $R^2$ of SFH fittings for both central and satellite galaxies.
The blue lines indicate statistics of central and isolated galaxies.
The oranges lines indicate statistics of satellite galaxies.
Solid lines denote the statistics of all samples.
The dashed lines show the CDF of sub-samples with stellar masses between $10^9M_{\odot}h^{-1}$ and $2\times10^9M_{\odot}h^{-1}$.
The dotted lines show the CDF of sub-samples with stellar masses between $10^{10}M_{\odot}h^{-1}$ and $2\times10^{10}M_{\odot}h^{-1}$.

Before conducting analysis,
it is important to evaluate the ability of our fitting method to accurately describe the SFHs of galaxies.
We examined the distribution of the goodness of fit for all samples, which is shown in FigR2.

For the total trends (grey line in FigR2), the median of $R^2$ value is  $0.53$, the peak of $R^2$ distribution is at $R^2=0.70$.
The median value here is consistent with [Iyer & Gawiser, 2017].
To quantitatively evaluate the performance at reconstructing the SFHs from mock catalogs, e.g., SAMs [Somerville & Davé, 2015],  hydrodynamic simulation MUFASA [Davé et al., 2017] and stochastic SFHs [Kelson, 2014], [Iyer & Gawiser, 2017] used $R^2$, which is in the same definition as this work, to quantify the accuracy of reconstruction.
The $R^2$ distribution plot indicates that the ability of reconstruction of our model reaches acceptable accuracy, therefore it can be used to fit the majority of SFHs of galaxies.
Nonetheless, the long tail in the $R^2<0$ region suggests that there is a small number of SFHs that cannot be fitted with our model.
Generally, fittings with multiple components exhibit good accuracy in fitting, where the fitting quality of Gaussian components is slightly superior to that of Gamma components, as shown consistently in the results presented in TabFitNumR.

The shape of the SFH curve of a galaxy is tightly related to its underlying physical processes.
In our study, the "goodness of fit" value can indicate to some extent whether the shape of the SFH curve of a galaxy approximates a certain type of function/distribution, specifically a combination of Gaussian and Gamma distributions.
Hence, exploring how the $R^2$ value relates to the properties of galaxies at $z=0$ will be useful in our investigation of star formation history.

FigR2Mass depicts the relationship between the $R^2$ value of SFH fittings and the stellar mass of their corresponding galaxies at $z=0$.
It is evident from the figure that higher mass galaxies exhibit better fitting quality.
This implies that the star formation history of high-mass galaxies conforms more to Gaussian or Gamma distribution shapes.
This $R^2 - M_*$ relation is partially ascribed to the resolution effect present in the simulation.
Low-mass galaxies contain a lower number of particles, leading to reduced accuracy and greater noise in reconstructing SFHs.
Consequently, the overall precision of the fitting is lowered.
Although, we have already set a minimum mass threshold to ensure sufficient galaxy particles to ease the resolution effect.
A more comprehensive assessment of this point would necessitate additional simulations with different resolutions, which are beyond the scope of this work.

FigR2Color depicts the relationship between the $R^2$ value of SFH fittings and the color of their corresponding galaxies at $z=0$.
It indicates that the correlation between SFH fitting goodness and color is weak for blue galaxies, where fitting goodness is relatively low.
Conversely, at the red end, our model yields better fitting results, as redder colors correspond better with SFH peak functions that resemble Gaussian or Gamma distributions.
We attribute this phenomenon to the quenching process that red galaxies underwent during stellar formation history, leading to distinctive SFH with peak-shaped functions that resemble Gaussian or Gamma functions more closely.

FigR2Len illustrates the relationship between the $R^2$ values and the age of galaxies.
Following [Nelson et al., 2018], the age of a galaxy is defined as the mass weighted average age of all stellar particles within it.
It represents the overall distribution of star formation time, as well as a effective length of the SFH
FigR2Len exhibits that young and old galaxies have better fitting goodness to our model.
This result is straightforward to comprehend.
Young galaxies have shorter and simpler shapes, which results from fewer physical mechanisms in the star formation process.
Old galaxies corresponds to those massive and red galaxies, which can achieve better fitting results to our model, as have been discussed in previous paragraphs.
For both old or young galaxies, their SFHs should have relatively significant peaks at early or late stages, which makes the fittings easier. On the other hand, the mass weighted age reaches half of the SFH length, approximately $6\ Gyr$ here, when the SFR distributes evenly accross the whole SFH. In this case, the SFH is hard to be fitted by Gaussian or Gamma distributions, thus results in the extremely low $R^2$ value at $Age\simeq 6 Gyr$. In subsequent content, the FigTypeSFHlen illustrates that the age distribution of SFHs with more components tends to concentrate around $6\ Gyr$, which also verifies the argument here.

We evaluated the fitting goodness of different galaxy types, as depicted
in FigR2Sat.
The galaxies were categorized into two groups: satellite galaxies and central/isolated galaxies.
In the simulation, we define the biggest galaxy in a dark matter halo as the central or isolated galaxy, and other galaxies as satellites.
The results show that satellite galaxies have a relatively better fitting.
We also plot the cumulative density function (CDF) for sub-samples with a stellar mass around $10^{9} M_{\odot}/h$ (dashed lines) and $10^{10} M_{\odot}/h$ (dotted lines).
The discrepancy between centrals and satellites still exists, negating the possibility of mass effect.
A plausible explanation is that satellites are more likely to be quenched in their history, resulting in SFH curves that are similar to Gaussian or Gamma-like peaks.
This explanation aligns with the color dependence of $R^2$.

Generally, only half of the galaxies possess SFH shapes that are consistent with the Gaussian-Gamma combined distribution.
High-mass galaxies, red galaxies, and satellite galaxies generally demonstrate a better fitting for these shapes.
We hypothesize that these results are mainly due to mechanisms that can shape distinctive peaks in SFHs, such as rapid increases and quenching of star formation.

§ DEPENDENCY OF SFH TYPES

The cumulative stellar mass function of galaxies at $z=0$.
Different colors represent galaxies with different SFH types.
The gray line shows the cumulative stellar mass function of all samples.
The SFHs types are classified by the number of components (left plot) or by the shapes of dominant components (right plot).
In the top sub-panel of each plot, the residuals to the cumulative mass function of all samples for the samples with different SFH types are given.

The cumulative density function of galaxies color ($g-r$) at $z=0$. Different colors represent different SFH types. The gray line shows the color distribution of all samples.
The SFHs types are classified by the number of components (left plot) or by the shapes of dominant components (right plot).
The top sub-panel of each plot shows the residuals to the color distribution of all samples for samples with different SFH types.

The cumulative density function of average galaxy age with different SFH types. Different colors of histograms represent different SFH types. The gray histogram shows the average galaxy age distribution of all samples.
The SFHs types are classified by the number of components (left plot) or by the shapes of dominant components (right plot).
The top sub-panel of each plot shows the residuals to the SFH length distribution of all samples for samples with different SFH types.

The SFHs that have been fitted are classified into different types based on the number of components and their major shapes, as explained in sec:fit.
This section aims to investigate how galaxy properties differ depending on the type of SFHs.
In order to ensure the validity of our study, we will only investigate well-matched samples of SFHs in this section.
This is because it is not convincing to attribute poorly fitted SFHs to specific types as defined by our classification.

FigTypeMass shows the stellar mass function of galaxies at $z=0$ corresponding to the different types of SFHs.
High-mass galaxies tend to have SFHs with more decomposition components, which aligns with our intuition.
Concerning the component shapes, high-mass galaxies tend to exhibit a prevalence of Gaussian components, whereas low-mass galaxies tend to be dominated by Gamma components.
The first two rows of FigExample1, FigExample2 and FigExample3 show SFH examples of low mass and high galaxies.

shows the distribution of galaxy color at $z=0$ corresponding to different types of SFHs.
Multiple-component SFHs are predominantly found in blue galaxies, while single-component SFHs are largely associated with red galaxies.
The color of galaxies shifts towards the blue area as the number of decomposed components in their history increases.
A simple understanding of this is that galaxies with fewer stellar formation epochs are less likely to experience star formation in their later stages, when compared to galaxies with multiple stellar formation epochs.
Consequently, the colors of the former tend to be red.
Galaxies with more Gamma components tend to have a slightly bluer color than others.
This is primarily because the typical pattern of star formation for many blue galaxies is a slow but continuous rise until recent time.
The Gamma distribution provides a good fit for this type of curve.
The third and fourth rows of FigExample1, FigExample2 and FigExample3 show SFH examples of red and blue galaxies.

FigTypeSFHlen shows the cumulative density distribution of the age of galaxies.
The age of a galaxy is defined as the mass weighted age of its all particles.
Galaxies with a single-component SFH are concentrated around two area, approximately age of 4 Gyr and 10 Gyr.
This corresponds to two types of galaxies.
The former represents young galaxies that have just formed and only experienced one period of star formation.
The latter refers to older quiescent galaxies that have remained in a quiescent state after an initial period of star formation,
hence having only one component in their early stage.
As the number of components increases, the distribution of galaxy ages gradually converges towards the central position of 6 Gyr.
The duration of star formation periods is determined by the underlying physical mechanisms thus will not extend infinitely.
Therefore, for SFHs with finite lengths (cannot exceed the age of the universe), there must be more effective star formation times to allow for additional star formation cycles and thus more SFH components.
When a galaxy maintains a constant rate of star formation from the beginning of the universe,
its average age of stellar particles tends to approach half of the age of the SFH length, approximately 6 Gyr in this case.
So, more components of SFHs indicates more effective star formation time across a galaxy's life.
The galaxy age does not show a strong correlation with the shape (Gaussian or Gamma) of the SFH components.
It implies that different physical drivers of star formation have relatively weak correlations with the galaxy age, indicating that they do not have a particular preference for the timing of occurrence.
The last two rows of FigExample1, FigExample2 and FigExample3 show SFH examples of  galaxies with large and small ages.

FigTypeCS shows the ratio of the number of central and isolated galaxies to satellite galaxies in galaxies with different SFH types.
The horizontal dashed line represents the overall sample ratio.
Data points above this line indicate fewer satellite galaxies in this subsample,
and vice versa.
The proportion of central and isolated galaxies increases with increasing decomposed SFH components, indicating more star formation epochs within central and isolated galaxies.
Because most of the satellite galaxies are quenched in late stages, they naturally have fewer star formation epochs.
The proportion of satellite galaxies with Gamma SFHs is slightly higher than those with Gaussian SFHs.
We believe that it is mainly due to the influence of the cluster stripping,  the star formation process in satellite galaxies is quenched more quickly, resulting in the asymmetric gamma distribution of the star formation history curve.

The ratio of central and isolated galaxies over satellite galaxies for different types of SFHs. The horizon line indicate the ratio for the whole sample.

Overall, high-mass, reddish, long-lived galaxies, and central and isolated galaxies tend to have more SFH components, suggesting more star formation periods.
Meanwhile, the shapes of the main components of SFHs in these galaxies tend to be slightly Gaussian like.
Nevertheless, there is considerable overlap in attribute distributions of galaxies with distinct SFH types, make it a challenging to establish a clear standard for strictly constraining the shape and number of SFH components to apply it in reconstructing the SFH of a given galaxy.

Examples of C1, C2 and C3 type SFHs with different galaxy mass, color and age. Each column show one kind of SFHs as labeled on the top. Each row represent one galaxy property range as labeled on the left. In each sub figure, the solid blue line represents the origin SFH, and the dashed orange line represents the best fit to it. The colored dashed lines represent decomposed components.

Examples of C4, C5 and C6 type SFHs with different galaxy mass, color and age. Each column show one kind of SFHs as labeled on the top. Each row represent one galaxy property range as labeled on the left.
The blue solid line represents the actual SFH data, and the orange dashed line represents the corresponding fitted curve.
The dotted lines show the decomposition of components in the fitted curve.

Examples of $\Gamma$, G$\Gamma$ and G type SFHs with different galaxy mass, color and age. Each column show one kind of SFHs as labeled on the top. Each row represent one galaxy property range as labeled on the left. In each sub figure, the solid blue line represents the origin SFH, and the dashed orange line represents the best fit to it. The colored dashed lines represent decomposed components.

§ EPISODES OF SFH

Both Gaussian and Gamma distributions have well-defined shapes,
hence, the fitting results can be used to characterize the episodes of each star formation history.
By exploring the parameters of the fitting function, the peak position and full width at half maximum (FWHM) of each episode can be easily identified and used to further understand the occurrence and time scale of star formation process in galaxies.
The position and width of SFH peaks are key parameters in this analysis, as they help us gain insight into the star formation history of galaxies.
We continue to constrain our samples with $R^2 \ge 0.5$ in this section.

§.§ Number of episodes

In this section, we recalculated the number of peaks in the SFH to determine the number of star formation episodes present in a given SFH.
This is because some fitted components may overlap with others.
We combine two fitted components into one peak when the distance between their peaks is smaller than the minimum FWHM of them.
The numbers of SFHs with varying numbers of peaks obtained after applying this procedure are presented in the third column of TabPeakNum, with the numbers of SFHs having the same number of components shown in the second column.
The table indicates that for each type of multi-component SFH, there are approximately $250$ SFHs whose components can be combined, which is quite a small fraction to the total number of SFHs.
This suggests that in most cases, the number of fitted components is consistent with the number of star formation episodes in a SFH.

Number of SFHs with different numbers of fitting components (left) and actual peaks (right), limited to the sample with fitting goodness of $R^2>0.5$.

$N$     lSFH with
$N$ components     lSFH with
$N$ peaks

1       8162                       8406

2       5786                       5868

3       3680                       3626

4       2446                       2388

5       1671                       1779

6       1484                       1162

TabPeakNum indicates that $64.9\%$ of the SFHs have a better fit with two or more peaks.
We examined the intrinsic peak numbers of SFHs in the simulation and confirmed that the reconstructed multi-peaks SFH fraction aligns with the intrinsic fraction.
However, these results appear to contradict the findings of [Iyer & Gawiser, 2017] and [Iyer et al., 2019], who reported a fraction of multi-episode SFHs ranging from $10\%$ to $20\%$.
They validated their method using hydro simulation and SAMs data, which derived a similar fraction from intrinsic peaks.

The main purpose of [Iyer & Gawiser, 2017] and [Iyer et al., 2019] was to constrain the SFH through SED fitting.
It is hard to extract information about all stages in a galaxy's life via SED.

The accuracy of fitting the episodes of SFH is highly sensitive to the timing, duration, and strength of the star formation rate, all of which are challenging to reconstruct solely from SED photometry/spectroscopy.
Due to the degeneracy of stellar population on spectrum and the noise, the reconstructed SFH from SED is always more smoothed-out compared to a direct fit to the original SFH.
Another possibility is that the discrepancy arises from the method used to count the intrinsic peaks.
However we speculate this possibility.
Considering the intrinsic peak number, it is unlikely that the simulation reproduces SFHs with significantly more peaks than those in the MUFASA simulation and SAMs.

We conducted tests using different sets for intrinsic peak counting, and the results are listed in Appendix <ref>.
In summary, [Iyer & Gawiser, 2017] and [Iyer et al., 2019] find less peaks, which is reasonable considering the nature of reconstructing SFHs from SED. Take advantages of direct fitting,
our approach attempts to distinguish more small peaks compared with the work in [Iyer & Gawiser, 2017, Iyer et al., 2019].
This approach only concerns the performance of direct fitting to SFH.
Whether the same performance can be achieved in observations and whether such capability can benefit studies are questions necessitate further exploration and analysis.

§.§ Timing and time scale of episodes

FigPeaks shows the distribution of peak positions and widths of SFHs.
Each data point in the figure corresponds to one peak.
In order to provide readers with a better understanding of the distribution of different intensities of peaks in the multiple-component SFH, we ranked the peaks within the same SFH according to their height and group peaks with same index as sub-samples.
Scatters and histograms with different colors in FigPeaks represent the results of peaks with different index.
However, we found that the position and width distributions of peaks across different orders are remarkably similar.
Therefore, in all subsequent discussions, unless otherwise stated, the subject of the discussion is all peaks.

The distribution of SFH peak position $t_{L,peak}$ reveals two conspicuous peaks. The first occurs at approximately $t_{L,peak} \approx 0$ Gyr, while the second arises at $t_{L,peak} \approx 10$ Gyr ($z\approx 1.8$).
The peak occurring at $t \approx 10\ Gyr$ is consistent with the general trend of cosmic star formation rate, which found that the star formation rate density peaks at $z \approx 2$ <cit.>.
The peak at $t \approx 0\ Gyr$ corresponds to SFH components keep rising.
[Iyer & Gawiser, 2017] likewise documented a significant fraction of recent peaks.

The FWHM of most peaks concentrate in the region of $0\ Gyr < FWHM < 2\ Gyr$, and reach the maximum of distribution at $FWHM \simeq 1.2\ Gyr$.
Many works emphasis the importance of SFR regulators in this timescale <cit.>, which correspond to various physical process like gas inflow and outflow, AGN feedbacks, galactic wind, giant molecular cloud life cycles, gas recycling, halo dynamical time-scale, etc. <cit.>.
The peak of FWHM distribution of $1.2\ Gyr$ coincides with the time scale of dark matter halo's dynamical timescale[Mo et al., 2010, Lilly et al., 2013].
This suggests that the dynamical process of dark matter halos significantly affects the star formation process of galaxies.

The FWHM distribution also includes a part of peaks featuring wider widths that correspond to long-term star formation processes, primarily associated with mergers[Robertson et al., 2006, Jiang et al., 2008, Boylan-Kolchin et al., 2008, Hani et al., 2020], metallicity evolution[Torrey et al., 2018], and galaxy quenching[Sales et al., 2015, Nelson et al., 2018, Rodríguez Montero et al., 2019, Wright et al., 2019].
However, the amount of the wide peaks is $1$ to $2$ orders of magnitude smaller than that of the narrow peaks.
It worth noting that, from the aspect of algorithm, our method is designed to get wider peaks to better represent the whole trend of the SFH of a galaxy.
But the FWHM distribution in FigPeaks still prefers narrow peaks.
Therefore, it demonstrates the importance of short time-scale processes in shaping the SFH of the galaxy.

The joint distribution of SFH width and position shows three major branches, as indicated by the three thick dashed lines in FigPeaks.
One of these branches is located vertically near $t_{L, peak}=0$, which corresponds to galaxies exhibiting a sustained increase in star formation rate until present time.
Another branch manifests in the form of a horizontal line with $2\ Gyr$ width near $FWHM \simeq 1.2\ Gyr$.
This branch suggests that a large number of star formation processes occurred within the timescale of $0$ to $2\ Gyr$.
These processes do not exhibit any clear preference for specific timing of occurring.
The third branch occurs like a negative correlation between $FWHM$ and $t_{L,peak}$.
For these components, the summation of looking back time of peaks position and $FWHM$ of peak width is close to the length of the SFH they belong to, namely $t_{L,peak} + FWHM \approx SFH\ Length$ .
These components are in line with the process of sustained SFR increase from the birth of the galaxy, succeeded by a decline at a specific point in time.
In FigPeaks, the line $FWHM=13.6-t_{L,peak}$ is shown as a reference, not meaning the regression for these group of components.

The position and width of peaks of SFHs.
$t_{L,peak}$ is the lookback time of the peaks.
$FWHM$ is the full width at half maximum of the peaks.
In the middle panel, each point represents data of one peak.
The upper panel shows the histogram of the peak positions,
and the right panel shows the histogram of the peak widths.
For each SFH, we sort the peaks by their height and distinguish the distribution of the peaks in different colors according to their order in the SFH.
'Peak 1' is the strongest peak, 'Peak 2' is the second strongest peak, and so on.
The three wide dashed lines in the figure correspond to $FWHM=13.6 - t_{L,peak}$, $t_{L, peak}=0$ and $FWHM=1$.

It is important to note that the bimodal distribution of peak positions does not imply the distance between adjacent peaks in a single SFH.
Within a single SFH, it is found that adjacent peaks still tend to occur at shorter time intervals.
This can be seen in Figure <ref>.
The median peak separation, accompanied by a $1 \sigma$ scatter, is $2.07^{+1.76}_{-0.92}\ Gyr$.
This value is significantly smaller than the one measured in the study by [Iyer et al., 2019].
In their study, the peak spacing was determined to be $0.42^{+0.15}_{-0.10}t_{univ}$, where $t_{univ}$ is the age of the universe at the time of observation.
Their estimate corresponds to $5.72^{+2.04}_{-1.36}\ Gyr$ in the context of this work.
However, the SFH reconstruction performed by [Iyer et al., 2019] used a Gaussian process, which is different from the methodology used in this work.
In their study, [Iyer et al., 2019] attempted to use the fewest possible number of time points, which resulted in the smoothing of the SFH and the omission of certain short-term fluctuations.
Additionally, as mentioned in the previous section, the SED fitting process naturally leads to fewer peaks being found when reconstructing the SFH.
As a consequence, larger intervals between peaks are expected.
On the other hand, although similar in methodology, the results from [Iyer & Gawiser, 2017] are not as convincing for comparison.
According to Figure 15 in [Iyer & Gawiser, 2017], the median peak separation is about $0.07 t_{univ}$, which is much smaller than the value found in this work and in [Iyer et al., 2019].
However, there are only a few tens of SFHs with two episodes in [Iyer & Gawiser, 2017], which reduces the statistical significance.
In summary, comparing peak separation between different studies is challenging due to differences in data and fitting methodologies.
These differences warrant further investigation.

The histogram of separations between neighboring peaks.

§.§ Mass dependency

The dependence of SFH peak position (top row), FWHM (middle row) and neighboring peak separation on stellar mass at $z=0$.
The colors correspond to logarithmic galaxy count in the corresponding SFH-properties-stellar-mass bins.
The median values of the peak position, FWHM, or peak separation in different stellar mass bins are displayed by the gray lines.
The error bars denote the range of the same value from $20\%$ to $80\%$ within the bin.
The statistics for all galaxies are presented in the left row, those for central and isolated galaxies are in the middle row, while satellite galaxies are on the right row.

The SFH of a galaxy is strongly associated with the galaxy properties.
Consequently, it is reasonable to anticipate the presence of correlations between the extracted SFH features(peak position, $FWHM$ and peaks separation) and the galaxy properties, including stellar mass and color.
Figure <ref> illustrates the relationship between the stellar mass at $z=0$ and the SFH features.
For each SFH peak, we record its position, $FWHM$ and the stellar mass of the galaxy it belongs to.
Also, we select all neighboring peak pairs from all SFHs and record their separation and the stellar mass of the galaxy they belong to.
Then we plot the 2D histograms of these SFH features and the stellar mass of corresponding galaxy.

The peak position of SFH in a galaxy exhibits a U-shaped relationship with the galaxy's mass.
The peak of SFH occurs earlier in galaxies with high and low masses, while it occurs later in galaxies with moderate masses.
Generally, the formation of high-mass galaxies requires larger dark matter halos with deeper potential wells to attract more gas.
The gas accretion into these strong potential wells triggers star formation earlier.
As the galaxy evolves, mechanisms such as AGN feedback and metal enrichment gradually strengthen to suppress the process of stellar formation.
This suppression effect is more pronounced in higher mass galaxies, preventing their star formation rate from surpassing the early levels in subsequent times.
The combination of these mechanisms leads to the earlier appearance of the SFH peaks in higher mass galaxies.
On the other hand, for galaxies at the lowest mass end, their limited gas reservoirs result in an early cessation of stellar formation, leading to smaller final masses.
Consequently, the peak position of their SFH remains in the early stages.
As for satellite galaxies, their stellar formation process stops earlier due to accretion and gas stripping from host galaxies, causing the peak position of SFH in moderate mass galaxies to shift towards earlier times

The correlation between the width of SFH peak and galaxy mass is relatively weak.
With the exception of galaxies at low and high mass ends, there is a slight decrease in the $FWHM$ as galaxy mass increases.
Analyzing the joint distribution of FWHM and mass, we observe a reduction in the number of broader peaks with the  mass range of approximately $9 \lesssim log_{10}(M_{*}/(h^{-1}M_{\odot})) \lesssim 11$.
This suggests a decreasing influence of the long time-scale star formation driving mechanisms on high-mass galaxies, implying that short time-scale process like feedbacks and gas recycling become more active in high-mass galaxies.
Both at the high and low mass ends, the median $FWHM$ experiences an increase, primarily due to the decline in the number of most narrow peaks in the SFHs.
The SFH of massive galaxies involves a large number of physical processes,
resulting in a more complex SFH shape and more decomposable components.
As a result, our method tends to exclude short time-scale components when there are too many pronounced components to be extracted.
On the other hand, the evolution process in low-mass galaxies is less complex with fewer fluctuations, leading to fewer instances of fitting with narrow peaks.

Excluding the highest and lowest mass galaxies, the peak separation of SFH decreases slightly with increasing galaxy mass.
Based on the context of galaxy SFR rejuvenation [Fang et al., 2012], it can be inferred that larger the galaxy mass results in a shorter time required to return to a state of star formation.
This is attributed to high-mass galaxies having stronger gas accretion rates and a higher chance of merging,
facilitating the acquisition of necessary fuel for the star formation process.
When distinguishing central/isolated galaxies from satellite galaxies, the correlation between peak separation and galaxy mass is weaker for satellite galaxies.
We propose two reasons for this.
On one hand, when satellite galaxies are influenced by central galaxies and experience gas stripping, the rise in their star formation rate is halted earlier, causing an earlier occurrence of the peak position in the star formation process.
This suppressing effect is more pronounced in low-mass galaxies.
Therefore, the larger peak separation of SFH at low mass is reduced more, resulting in a flatter relationship curve.
On the other hand, the frequency of disturbance and stripping events for satellite galaxies is primarily dependent on the external environment rather than their own properties.
This introduces more randomness and disturbance to their star formation process, thereby weakening the relationship between their SFH characteristics and their own stellar mass.

§.§ Color dependency

Figure <ref> illustrates the relationship between the galaxy color at $z=0$ and the SFH features.
For each SFH peak, we record its position, $FWHM$ and the $g-r$ color of the galaxy it belongs to.
Also, we select all neighboring peak pairs from all SFHs and record their separation and the $g-r$ color of the galaxy they belong to.
Then we plot the 2D histograms of these SFH features(peak position, $FWHM$ and peaks separation) and the color of corresponding galaxies.

The dependence of SFH peak position (top row), FWHM (middle row) and neighboring peak separation on galaxy color $g-r$ at $z=0$.
The colors on plots correspond to logarithmic galaxy count in the corresponding SFH-properties-color bins.
The median values of the peak position, FWHM, or peak separation in different color bins are displayed by the gray lines.
The error bars denote the range of the same value from $20\%$ to $80\%$ within the bin.
The statistics for all galaxies are presented in the left row, those for central and isolated galaxies are in the middle row, while satellite galaxies are on the right row.

According to Figure 14, there exists a clear linear correlation between the color of galaxies and the position of the SFH peak.
It indicates that galaxies with a redder color tend to have earlier peaks,  which is naturally attributed to the fact that galaxies formed in earlier tend to have older stellar ages, thus displaying a redder color.

The width of the SFH peaks exhibits a tri-modal distribution with the colors of galaxies.
At the blue end with $g-r<0.2$, the peaks have smaller FWHM.
This is because the most blue galaxies are still in the early stage of star formation, with their SFR rapidly increasing and their peak positions mainly located at a redshift of $0$.
Therefore, the fitted peaks have narrower widths.
When the $g-r$ value is within the range of $0.2$  to  $0.6$, the width of the peaks has no correlation with the colors of the galaxies.
When the galaxy color exceeds $0.6$, a large number of wide peaks appear in satellite galaxies.
These are SFH components in satellite galaxies whose star formation processes are regulated by self quenching.
This result is consistent with previous work which summarized the quenching time scale in other simulations.
[Sales et al., 2015] found that the quenching time scale of satellite galaxies is is about $2$ to $5$ Gyr in Illustris.
[Wright et al., 2019] found that the quenching time scale is about $2.5\sim 3.3$ Gyr, and extending out to $10$ Gyr, in EAGLE.

The peaks separation in the SFH also has a clear relationship with galaxy color.
In redder galaxies, the stellar formation process is concentrated in earlier periods, resulting in a more concentrated peak position of star formation.

§ CONCLUSION

To summarize, this study has undertaken the following objectives and obtained statistical findings concerning the SFHs of galaxies from the simulation:

* We employed a function combined by up to $6$ components, described by Gaussian or Gamma distributions, to construct fittings  for the SFHs in the simulation.
Our fitting method can well recover the trends of the SFH, though the goodness of fit can depend on galaxy mass, color, history length, and the type of galaxy.

* We analyzed the SFHs  that are well-fitted by the Gaussian and Gamma distribution combinations.
Our finding indicates that low-mass galaxies, red galaxies, and satellite galaxies tend to have fewer components and to be majored by Gamma distribution, while high-mass galaxies, blue galaxies and central galaxies tend to have more components and to be majored by Gaussian distribution.

* We calculated the distribution of the position, width, and separation of SFH peaks.
The peak position distribution is concentrated at $t_L = 0\ Gyr$ and $t_L = 10\ Gyr$.
The width distribution of peaks is mainly concentrated at $FWHM \sim 0\ Gyr - 2\ Gyr$.
The joint distribution of peak position and width is concentrated in three regions, corresponding to different types of star formation modes.
In addition, we found that the separation between adjacent peaks in the TNG simulation is concentrated around $2\ Gyr$, much smaller than previously reconstructed SFHs from observations.

* We analyzed the relationships between three temporal attributes of the SFH, namely peak position, peak width, and peak spacing, and two properties of galaxies, namely stellar mass and galaxy color.
Our findings indicate that the peak position and peaks separation are related with stellar mass and galaxy color.
The peak width shows limited relationship with stellar mass and galaxy color. This is reasonable because the peak width reflects the time scale of basic physical processes which drives star formation.

Our model can accurately reproduce the trend of galaxies' stellar formation evolution.
The method of multi-component decomposition can greatly help us understand the different physical driving factors of star formation in galaxies.
Although more than half of the SFHs have low goodness of fit, this is mainly due to the excessive fluctuations in the SFHs, which makes it extremely difficult to fit them for all kinds of functional templates.
Of course, the method we used is based on simulation data, the gaps between our method and SFH reconstruction methods suit for observational galaxies can not be neglected.
Since it is difficult to give an absolute resolution criterion for the shape types and number of fitted components, our model will require better algorithms and more advanced computing power to determine these factors when applied to observational data, which will be further developed and improved in future work.

§ ACKNOWLEDGEMENTS

This work is supported by “The Major Key Project of PCL”.
YW thanks the support from Kunlun server in SYSU.
The authors thanks the Illustris and Illustirs-TNG projects for providing the data.
The authors thanks the referee for the valuable comments.

[Abramson et al., 2015]
Abramson, L. E., Gladders, M. D., Dressler, A., et al. 2015, , 801, L12, 10.1088/2041-8205/801/1/L12

[Behroozi et al., 2010]
Behroozi, P. S., Conroy, C., & Wechsler, R. H. 2010, , 717, 379, 10.1088/0004-637X/717/1/379

[Behroozi et al., 2013]
Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, , 770, 57, 10.1088/0004-637X/770/1/57

[Bellstedt et al., 2020]
Bellstedt, S., Robotham, A. S. G., Driver, S. P., et al. 2020, , 498, 5581, 10.1093/mnras/staa2620

[Benítez-Llambay et al., 2015]
Benítez-Llambay, A., Navarro, J. F., Abadi, M. G., et al. 2015, , 450, 4207, 10.1093/mnras/stv925

[Boylan-Kolchin et al., 2008]
Boylan-Kolchin, M., Ma, C.-P., & Quataert, E. 2008, , 383, 93, 10.1111/j.1365-2966.2007.12530.x

[Buat et al., 2008]
Buat, V., Boissier, S., Burgarella, D., et al. 2008, , 483, 107, 10.1051/0004-6361:20078263

[Carnall et al., 2019]
Carnall, A. C., Leja, J., Johnson, B. D., et al. 2019, , 873, 44, 10.3847/1538-4357/ab04a2

[Carnall et al., 2018]
Carnall, A. C., McLure, R. J., Dunlop, J. S., & Davé, R. 2018, , 480, 4379, 10.1093/mnras/sty2169

[Chauke et al., 2018]
Chauke, P., van der Wel, A., Pacifici, C., et al. 2018, , 861, 13, 10.3847/1538-4357/aac324

[Choi et al., 2014]
Choi, E., Naab, T., Ostriker, J. P., Johansson, P. H., & Moster, B. P. 2014, , 442, 440, 10.1093/mnras/stu874

[Ciesla et al., 2016]
Ciesla, L., Boselli, A., Elbaz, D., et al. 2016, , 585, A43, 10.1051/0004-6361/201527107

[Collaboration et al., 2016]
Collaboration, P., Ade, P. A. R., Aghanim, N., et al. 2016, , 594, A13, 10.1051/0004-6361/201525830

[Conroy, 2013]
Conroy, C. 2013, , 51, 393, 10.1146/annurev-astro-082812-141017

[Conroy et al., 2014]
Conroy, C., Graves, G. J., & van Dokkum, P. G. 2014, , 780, 33, 10.1088/0004-637X/780/1/33

[Crain et al., 2015]
Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, , 450, 1937, 10.1093/mnras/stv725

[Davé et al., 2017]
Davé, R., Rafieferantsoa, M. H., Thompson, R. J., & Hopkins, P. F. 2017, , 467, 115, 10.1093/mnras/stx108

[Diemer et al., 2017]
Diemer, B., Sparre, M., Abramson, L. E., & Torrey, P. 2017, , 839, 26, 10.3847/1538-4357/aa68e5

[Digby et al., 2019]
Digby, R., Navarro, J. F., Fattahi, A., et al. 2019, , 485, 5423, 10.1093/mnras/stz745

[Donnari et al., 2019]
Donnari, M., Pillepich, A., Nelson, D., et al. 2019, , 485, 4817, 10.1093/mnras/stz712

[Fang et al., 2012]
Fang, J. J., Faber, S. M., Salim, S., Graves, G. J., & Rich, R. M. 2012, , 761, 23, 10.1088/0004-637X/761/1/23

[Finlator et al., 2007]
Finlator, K., Davé, R., & Oppenheimer, B. D. 2007, , 376, 1861, 10.1111/j.1365-2966.2007.11578.x

[Gallazzi et al., 2005]
Gallazzi, A., Charlot, S., Brinchmann, J., White, S. D. M., & Tremonti, C. A. 2005, , 362, 41, 10.1111/j.1365-2966.2005.09321.x

[Gavazzi et al., 2002]
Gavazzi, G., Bonfanti, C., Sanvito, G., Boselli, A., & Scodeggio, M. 2002, , 576, 135, 10.1086/341730

[Gladders et al., 2013]
Gladders, M. D., Oemler, A., Dressler, A., et al. 2013, , 770, 64, 10.1088/0004-637X/770/1/64

[Hahn et al., 2019]
Hahn, C., Starkenburg, T. K., Choi, E., et al. 2019, , 872, 160, 10.3847/1538-4357/aafedd

[Hani et al., 2020]
Hani, M. H., Gosain, H., Ellison, S. L., Patton, D. R., & Torrey, P. 2020, , 493, 3716, 10.1093/mnras/staa459

[Hopkins et al., 2018]
Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, , 480, 800, 10.1093/mnras/sty1690

[Iyer & Gawiser, 2017]
Iyer, K., & Gawiser, E. 2017, , 838, 127, 10.3847/1538-4357/aa63f0

[Iyer et al., 2019]
Iyer, K. G., Gawiser, E., Faber, S. M., et al. 2019, , 879, 116, 10.3847/1538-4357/ab2052

[Iyer et al., 2020]
Iyer, K. G., Tacchella, S., Genel, S., et al. 2020, , 498, 430, 10.1093/mnras/staa2150

[Janowiecki et al., 2017]
Janowiecki, S., Salzer, J. J., van Zee, L., Rosenberg, J. L., & Skillman, E. 2017, , 836, 128, 10.3847/1538-4357/836/1/128

[Jiang et al., 2008]
Jiang, C. Y., Jing, Y. P., Faltenbacher, A., Lin, W. P., & Li, C. 2008, , 675, 1095, 10.1086/526412

[Jiménez-López et al., 2022]
Jiménez-López, D., Corcho-Caballero, P., Zamora, S., & Ascasibar, Y. 2022, , 662, A1, 10.1051/0004-6361/202141338

[Johnson et al., 2013]
Johnson, B. D., Weisz, D. R., Dalcanton, J. J., et al. 2013, , 772, 8, 10.1088/0004-637X/772/1/8

[Joshi et al., 2021]
Joshi, G. D., Pillepich, A., Nelson, D., et al. 2021, , 508, arXiv:2101.12226, 10.1093/mnras/stab2573

[Katsianis et al., 2020]
Katsianis, A., Xu, H., Yang, X., et al. 2020, arXiv e-prints, arXiv:2010.08173.

[Kauffmann, 2014]
Kauffmann, G. 2014, , 441, 2717, 10.1093/mnras/stu752

[Kelson, 2014]
Kelson, D. D. 2014, arXiv e-prints, arXiv:1406.5191.

[Lee et al., 2010]
Lee, S.-K., Ferguson, H. C., Somerville, R. S., Wiklind, T., & Giavalisco, M. 2010, , 725, 1644, 10.1088/0004-637X/725/2/1644

[Leja et al., 2019]
Leja, J., Carnall, A. C., Johnson, B. D., Conroy, C., & Speagle, J. S. 2019, , 876, 3, 10.3847/1538-4357/ab133c

[Leja et al., 2017]
Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, , 837, 170, 10.3847/1538-4357/aa5ffe

[Lilly et al., 2013]
Lilly, S. J., Carollo, C. M., Pipino, A., Renzini, A., & Peng, Y. 2013, , 772, 119, 10.1088/0004-637X/772/2/119

[Lower et al., 2020]
Lower, S., Narayanan, D., Leja, J., et al. 2020, , 904, 33, 10.3847/1538-4357/abbfa7

[Lu et al., 2016]
Lu, Y., Benson, A., Mao, Y.-Y., et al. 2016, , 830, 59, 10.3847/0004-637X/830/2/59

[Lu et al., 2015]
Lu, Z., Mo, H. J., Lu, Y., et al. 2015, , 450, 1604, 10.1093/mnras/stv667

[Maraston et al., 2010]
Maraston, C., Pforr, J., Renzini, A., et al. 2010, , 407, 830, 10.1111/j.1365-2966.2010.16973.x

[Marinacci et al., 2018]
Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, , 480, 5113, 10.1093/mnras/sty2206

[Matthee & Schaye, 2019]
Matthee, J., & Schaye, J. 2019, , 484, 915, 10.1093/mnras/stz030

[Mo et al., 2010]
Mo, H., van den Bosch, F. C., & White, S. 2010, Galaxy Formation and Evolution (Cambridge University Press).

[Naiman et al., 2018]
Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, , 477, 1206, 10.1093/mnras/sty618

[Nelson et al., 2018]
Nelson, D., Pillepich, A., Springel, V., et al. 2018a, , 475, 624, 10.1093/mnras/stx3040

[Nelson et al., 2018]
Nelson, D., Kauffmann, G., Pillepich, A., et al. 2018b, , 477, 450, 10.1093/mnras/sty656

[Nelson et al., 2019]
Nelson, D., Springel, V., Pillepich, A., et al. 2019, Computational Astrophysics and Cosmology, 6, 2, 10.1186/s40668-019-0028-x

[Ocvirk et al., 2006]
Ocvirk, P., Pichon, C., Lançon, A., & Thiébaut, E. 2006, , 365, 46, 10.1111/j.1365-2966.2005.09182.x

[Pacifici et al., 2012]
Pacifici, C., Charlot, S., Blaizot, J., & Brinchmann, J. 2012, , 421, 2002, 10.1111/j.1365-2966.2012.20431.x

[Pacifici et al., 2016]
Pacifici, C., Kassin, S. A., Weiner, B. J., et al. 2016, , 832, 79, 10.3847/0004-637X/832/1/79

[Panter et al., 2007]
Panter, B., Jimenez, R., Heavens, A. F., & Charlot, S. 2007, , 378, 1550, 10.1111/j.1365-2966.2007.11909.x

[Papovich et al., 2011]
Papovich, C., Finkelstein, S. L., Ferguson, H. C., Lotz, J. M., & Giavalisco, M. 2011, , 412, 1123, 10.1111/j.1365-2966.2010.17965.x

[Pillepich et al., 2018]
Pillepich, A., Nelson, D., Hernquist, L., et al. 2018a, , 475, 648, 10.1093/mnras/stx3112

[Pillepich et al., 2018]
Pillepich, A., Springel, V., Nelson, D., et al. 2018b, , 473, 4077, 10.1093/mnras/stx2656

[Reddy et al., 2012]
Reddy, N. A., Pettini, M., Steidel, C. C., et al. 2012, , 754, 25, 10.1088/0004-637X/754/1/25

[Robertson et al., 2006]
Robertson, B., Bullock, J. S., Cox, T. J., et al. 2006, , 645, 986, 10.1086/504412

[Rodriguez-Gomez et al., 2015]
Rodriguez-Gomez, V., Genel, S., Vogelsberger, M., et al. 2015, , 449, 49, 10.1093/mnras/stv264

[Rodríguez Montero et al., 2019]
Rodríguez Montero, F., Davé, R., Wild, V., Anglés-Alcázar, D., & Narayanan, D. 2019, , 490, 2139, 10.1093/mnras/stz2580

[Sales et al., 2015]
Sales, L. V., Vogelsberger, M., Genel, S., et al. 2015, , 447, L6, 10.1093/mnrasl/slu173

[Schaye et al., 2015]
Schaye, J., et al. 2015, , 446, 521, 10.1093/mnras/stu2058

[Schreiber et al., 2018]
Schreiber, C., Glazebrook, K., Nanayakkara, T., et al. 2018, , 618, A85, 10.1051/0004-6361/201833070

[Simha et al., 2014]
Simha, V., Weinberg, D. H., Conroy, C., et al. 2014, arXiv e-prints, arXiv:1404.0402, 10.48550/arXiv.1404.0402

[Somerville & Davé, 2015]
Somerville, R. S., & Davé, R. 2015, , 53, 51, 10.1146/annurev-astro-082812-140951

[Springel, 2010]
Springel, V. 2010, , 401, 791, 10.1111/j.1365-2966.2009.15715.x

[Springel et al., 2001]
Springel, V., Yoshida, N., & White, S. 2001, , 6, 79, 10.1016/S1384-1076(01)00042-2

[Springel et al., 2018]
Springel, V., Pakmor, R., Pillepich, A., et al. 2018, , 475, 676, 10.1093/mnras/stx3304

[Tacchella et al., 2020]
Tacchella, S., Forbes, J. C., & Caplar, N. 2020, , 497, 698, 10.1093/mnras/staa1838

[Telles & Melnick, 2018]
Telles, E., & Melnick, J. 2018, , 615, A55, 10.1051/0004-6361/201732275

[Thomas et al., 2005]
Thomas, D., Maraston, C., Bender, R., & Mendes de Oliveira, C. 2005, , 621, 673, 10.1086/426932

[Tojeiro et al., 2007]
Tojeiro, R., Heavens, A. F., Jimenez, R., & Panter, B. 2007, , 381, 1252, 10.1111/j.1365-2966.2007.12323.x

[Torrey et al., 2018]
Torrey, P., Vogelsberger, M., Hernquist, L., et al. 2018, , 477, L16, 10.1093/mnrasl/sly031

[Vogelsberger et al., 2014]
Vogelsberger, M., Genel, S., Springel, V., et al. 2014, , 444, 1518, 10.1093/mnras/stu1536

[Walcher et al., 2011]
Walcher, J., Groves, B., Budavári, T., & Dale, D. 2011, , 331, 1, 10.1007/s10509-010-0458-z

[Weinberger et al., 2017]
Weinberger, R., Springel, V., Hernquist, L., et al. 2017, , 465, 3291, 10.1093/mnras/stw2944

[Wetzel et al., 2016]
Wetzel, A. R., Hopkins, P. F., Kim, J.-h., et al. 2016, , 827, L23, 10.3847/2041-8205/827/2/L23

[Wright et al., 2019]
Wright, R. J., Lagos, C. d. P., Davies, L. J. M., et al. 2019, , 487, 3740, 10.1093/mnras/stz1410

[Wuyts et al., 2011]
Wuyts, S., Förster Schreiber, N. M., van der Wel, A., et al. 2011, , 742, 96, 10.1088/0004-637X/742/2/96

[Zhou et al., 2020]
Zhou, S., Mo, H. J., Li, C., Boquien, M., & Rossi, G. 2020a, , 497, 4753, 10.1093/mnras/staa2337

[Zhou et al., 2020]
Zhou, Z.-B., Zhu, W., Wang, Y., & Feng, L.-L. 2020b, , 895, 92, 10.3847/1538-4357/ab8d32

§ TIME BIN RESOLUTION IN SFH RECONSTRUCTION

The SFHs of three galaxies obtained at different time bin resolutions (blue solid lines), along with the corresponding fitted curve (orange dashed lines) and the decomposition of components (colored dotted lines).
Each row contains the results for one time bin.
Each column contains the results for one galaxy.

The histograms of fitting goodness $R^2$.
The histograms of different colors represent the results of the SFHs obtained at different time bin resolutions.
The upper panel shows the histograms for high mass galaxies,
while the lower panel shows the histograms for low mass galaxies.

The histogram of the difference in the number of components in the SFH obtained at different time bin resolutions.
$C_{dt}$ represents the number of components in an SFH.
The x-axis shows the difference between the number of components at four other time bin resolutions and the number of components at a resolution of 100 Myr.
The y-axis represents the corresponding number of galaxies.
The upper panel shows the results for high-mass galaxies,
while the lower panel shows the results for low-mass galaxies.

The resolution of the time bin significantly influences the extraction of a galaxy's star formation history and the subsequent fitting work.
FigFitTimebin illustrates the SFHs (represented by the blue line) and the corresponding fitting curves (represented by the orange line) obtained using different time bin resolutions.
The figure reveals that a larger time bin results in a smoother SFH curve, albeit at the cost of losing finer details.

However, a smaller time bin is not always the optimal choice.
Given the finite number of particles in a galaxy during the simulation, an excessive number of time bins would lead to fewer particles in each bin, thereby causing significant numerical fluctuations.
This effect consequently reduces the goodness of fit for the SFHs.
As depicted in FigFitTimebin, when the width of the time bin is $10$ Myr, the noise in the SFH curve itself is considerably high.
Despite this, our fitting method can still identify the main components in the curve, although the goodness of fit is significantly reduced numerically.

Conversely, a larger time bin results in a fitted curve that aligns better with the SFH.
This is because using a larger time bin is akin to performing a smoothing process on the curve, which effectively eliminates noise.

In FigR2_timebin, we present the statistical goodness of fit of our method for SFHs with different time bins.
It can be observed that $R^2$ decreases as the time bin decreases.
This relationship is particularly pronounced in low-mass galaxies, primarily due to the numerical noise caused by the insufficient number of particles in these galaxies.

Another point of interest is whether our fitting results would be affected when different time bins are used.
According to our inspection results (partially shown in FigFitTimebin), although the number of components and the type of each component (Gaussian or Gamma) may vary, the major peaks can be well restored and maintain stable positions and widths.
The main differences lie in the handling of some small components, where different time bins can lead to different results.

We use the fitting results of a time bin of $100$ Myr as a benchmark and compare the number of components obtained by fitting the SFHs of other time bins with it.
As shown in FigCChange_timebin, most galaxies' SFH fitting does not change the number of components relative to when the time bin equals $100$ Myr ($C_{\Delta t}-C_{0.1}=0$), where $C_{\Delta t}$ is the number of fitting components of a single SFH when the time bin is $\Delta t$.

For high-mass galaxies, whether the time bin increases or decreases, the change in the number of components is mainly a decrease, i.e., the side where $C_{\Delta t}-C_{0.1}<0$ is higher than the other side.
As mentioned above, this is because the numerical noise brought about by the decrease in time bin can blur some small peaks, while the smoothing effect brought about by the increase in time bin may also erase these peaks.
This suggests that $100$ Myr is a good choice of time bin for high-mass galaxies in the TNG simulation.

For low-mass galaxies, we can see that the number of components tends to increase when the time bin increases. This suggests that the time bin of $100$ Myr may still be too small for these galaxies, and a slightly larger time bin matches their particle number better.

At last, readers should remember that the influence of time bin resolution mainly comes from the mass resolution limit in simulation.
In other theoretical galaxy models, i.e., SAMs, or observational data, the influence of mass resolution can be alleviated.

§ STELLAR MASS LOSS

The SFH obtained from the initial stellar mass (blue line) and the final stellar mass (orange line) of one galaxy.

Due to the loss of mass in stellar particles in simulation, the Star Formation Histories (SFHs) constructed based on the mass of particles at z=0 and those based on the initial mass of particles will differ.
FigSFHcompare illustrates the differences between these two types of SFHs.
As the age of each particle remains unchanged, with only the value of its stellar mass varying, the differences between the two SFHs are solely vertical shifts.
Consequently, the differences in the fitting results based on the two types of SFHs primarily lie in the amplitude of the components, with no changes in the position and width of each component.

§ FACTORS AFFECTING INTRINSIC SFH PEAK COUNTS

This table represents the factors that influence the intrinsic peak counts of SFH.
The last row provides the reference for the multi-episode fraction of reconstructed SFHs in this work.
Note that the multi-component method we proposed in this work does not require prominence parameter to define the peaks.
Note that the prominence parameter in this work is the difference of $SFR$ rather than $logSFR$, while the latter one  is used in [Iyer et al., 2019].

articles                      dataset                                z                        $M_*$                                    $N_{sample}$                 smooth method                                               peak finder prominence             $frac(N_{peak}>1$)

3*Iyer2017     MUFASA                                 1                        3*$\ge 10^9M_{\odot}$     1200                         3*not mentioned                              3*not mentioned     $12\% \sim 17 \%$

SAMs                                   1                                                                 1200                                                                                $15\% \sim 20 \%$

stochastic                             1                                                                 1200                                                                                $\sim 15 \%$

3*Iyer2019     3*MUFASA
& SAMs             $0.5$        3*$\ge 10^8M_{\odot}$     [vpos]3*$10000$  in
$0.5<z<3$     3*Gaussian
Processes     3*$log \frac{SFR_{peak}}{SFR_{min,local}} $
$> 1.5 + \frac{1.5}{4}log \frac{M_*}{10^8 M_{\odot}}$                   $\sim 20\%$

$1$                                                                                $10\% \sim 20\%$

$2$                                                                                $10\% \sim 15\%$

3*Test         3*TNG                   3*$0$     3*$\ge 10^8M_{\odot}$     3*$23229$     3*savgol_filter  a     Iyer2019 b                 $16.6\%$

$SFR_{max}/2$  c           $19.9\%$

$SFR_{max}/4 $                     $66.5\%$

This work                     TNG                                    0                        $\ge 10^8M_{\odot}$                      $23229$                      –                                                          –                                 $64.9\%$

a We use the savgol_filter function from scipy.signal module.

b This prominence is the same as that in the line of Iyer et al. 2019.

c $SFR_{max}$ is the maximum value of each smoothed SFH.

Two examples of the effects of prominence. In each figure, blue solid line represents a origin SFH from the TNG simulation.
The orange dashed line represents the smoothed curve using the scipy.signal.savgol_filter function with a window of $20$ points and polyorder $3$.
The peaks found with a prominence equivalent to that in [Iyer et al., 2019] are indicated by green circles.
Purple squares indicate the peaks found with a prominence of $SFR_{max}/2$,
while red triangles indicate the peaks found with a prominence of $SFR_{max}/4$.
To ensure the clarity of the markers, they have been slightly shifted upwards.
In the left figure, three different sets of markers can identify the same peaks in the smoothed SFH. However, in the right figure, only the smallest prominence is able to distinguish the smaller peak.

Section <ref> highlights a discrepancy between the ”true fraction" of multiple episodes SFHs in this study and those reported in [Iyer & Gawiser, 2017, Iyer et al., 2019].
To investigate this difference, we emulate the method of identifying intrinsic peaks as described in [Iyer & Gawiser, 2017] and [Iyer et al., 2019], and examine various sets of prominence for peak finding algorithm.
As demonstrated in TabPeakFrac, a higher prominence in peak finding algorithm results in a smaller fraction of multi-episode SFHs.
Based on the value of multi-episode fraction, the SFHs reconstructed by our multi-component method closely aligns with the intrinsic SFHs defined with with a prominence of $SFR_{peak}-SFR_{min,local}>SFR_{max}/4$, where $SFR_{max}$ is the max SFR value of one smoothed SFH.
Although the smoothing method and prominence in our test can not be aligned with [Iyer & Gawiser, 2017] and [Iyer et al., 2019] exactly, we can confirm that the prominence accepted in [Iyer et al., 2019] is larger than $SFR_{max}/4$ for the majority of SFHs in the TNG simulation.
We apply the prominence from [Iyer et al., 2019] to our data and find that the peaks defined by this prominence are more akin to the performance of the prominence of $SFR_{max}/2$ (refer to FigPromTest).


# Apples and Cars: a Comparison of Security

Zhendong Ma<EMAIL_ADDRESS>Bosch Engineering GmbHHolzkirchenGermany

###### Abstract.

Cybersecurity has gained importance for cars that increasingly rely on
software and networks. “Smartphone on wheels” is often used as an analogy to
highlight the need for security. As a high-value target of cyberattacks,
modern smartphones implement layers of protection. Automotive embedded systems
share many similarities with smartphones. We compare the security architecture
of an iPhone and a car to identify gaps and discuss the potentials for the
cars of the future.

security architecture, automotive, smartphone

††copyright: none††copyright: none

##  1\. Introduction

Modern cars are often referred to as “smartphone on wheels”, because the
embedded systems in a car are programmed and connected using similar
technologies that power a smartphone in a user’s pocket. Smartphones are a
high-value target of cyberattacks. Consequently, layers of security measures
are in place to protect the hardware, software, and data against a wide
spectrum of threats. The system and threat model of a connected car have many
similarities to a smartphone. Assuming that smartphones have more advanced
security, questions arise as what are the commonalities and differences in
securing a smartphone comparing to a connected car, and what can automotive
security learn from it? In this paper, we compare the security architecture of
an Apple iPhone(apple21, ) with that of typical connected cars. The comparison
shows gaps of the cars of today and potentials for the cars of tomorrow.

Figure 1 illustrates the system and threat model of both systems. The red
arrows indicate the attack surface. A car with a domain-based architecture
uses a gateway to divide Electronic Control Units (ECU) into domains according
to their functionalities. The gateway ECU routes and filters cross-domain
communications. The threat model includes attacks against the cloud backend,
network-based attacks via connectivity interface (i.e. the telematics ECU),
local attacks via on-board OBD port and the in-vehicle communication network,
and attacks on ECUs via physical access.

Figure 1. System and threat model of a car and an iPhone

An iPhone consists of hardware including a highly integrated System-on-Chip
(SoC) and peripherals, the operating system iOS, and Apps from Apple and
third-party. Its threat model includes attacks against the iCloud server,
networked attacks with or without user interaction (the so-called “zero-click
attack”), attacks via installed Apps, attacks on iOS, and hardware attacks via
physical access down to the silicon.

##  2\. Comparison of security

Hardware security iPhone’s threat model assumes that an attacker has physical
access to the internals of the SoC to intercept, extract and manipulate the
data on the bus and in the random-access memory (RAM) and non-volatile memory
(NVM). Inside the Apple SoC there is a coprocessor named _Secure Enclave_ with
its own microkernel operating system for cryptographic operation. It isolates
cryptographic secrets and operation from the application processor and the
main memory and protects the cryptographic module from side-channel attacks.
The Secure Enclave is equipped with tamper detection and power and clock
monitoring circuit to protect itself from tampering and dynamic power
analysis. Both the Secure Enclave and the application processor have its
hardware-based inline AES engine, which encrypts data written to memory. In
addition to encrypted memory, Secure Enclave also implements measures for
memory authentication and replay protection to ensure memory security within
the SoC. The build-in public key accelerator supports RSA and ECC encryption
and signing. An immutable Boot ROM provides a hardware root-of-trust for
system boot.

The ECUs inside the cars typically use microcontrollers ($\mu$C) for time-
critical functions such as steering or engine control, and microprocessors
($\mu$P) for computation intensive functions such as autonomous driving (AD)
or in-vehicle infotainment. The state-of-the-art $\mu$Cs are equipped with
Hardware Security Modules (HSM) for secure key storage and cryptographic
operation. A HSM has its own processor and RAM for a physically isolated
execution environment. It communicates with the main CPU via the system bus.
It has hardware accelerator for AES. Some models include hardware accelerator
for asymmetric algorithm such as RSA and ECC. Hardware security of $\mu$P is
typically achieved by trusted execution environment (TEE) technology, which is
a virtual execution environment created on the main processor and supported by
hardware-based access control.

In comparison, Secure Enclave is embedded in the SoC package. HSM is located
outside the main processors of the $\mu$C. TEE is a software solution
supported by hardware. The more a secure execution environment can be
isolated, the higher the security level it will achieve.

System and software security iPhone’s system security ensures that the device
runs only authentic software from Apple and the system maintains its integrity
during runtime. Secure boot ensures that the software in each stage verifies
the integrity of the next stage before handing over control. The iBoot
bootloader is a security-enhanced build that hardens the software against
typical C vulnerabilities such as buffer-overflow, heap exploitation, and type
confusion. iOS uses several integrity protection mechanisms. _Kernel integrity
protection_ creates a protected physical memory region for the kernel and its
extensions. The region is locked after boot. _Fast permission restrictions_
use a CPU register to restrict memory read/write/execute permissions on a per-
thread basis. During boot, iBoot loads coprocessor’s firmware into protected
memory region and the system prevents executable and writeable mapping on that
region. _Pointer authentication codes_ protect program control flows by
preventing modification of function pointers and return addresses. All
executable code are signed and verified before execution in the memory. _Page
protection layer_ ensures that only it can alter memory pages containing user
code such that an attacker cannot tamper the code after signature
verification. Third-party apps are sandboxed to restrict their access to other
files and system resources. Further apps memory protections include Address
Space Layout Randomization (ASLR) and Execute Never (XN) which prevent code
execution from certain memory pages. Apple over-the-air (OTA) update process
includes additional parameters to personalize the update package for each
device and to prevent replay and downgrade attack. To request an update, an
iPhone sends a nonce and its unique device ID to the server. The server signs
the update including the installation bundle, the nonce, and the device ID.

In automotive domain, secure boot gradually becomes a standard security
feature to ensure software authenticity and integrity inside ECUs. Although
some $\mu$Cs include memory protection unit that restricts read/write access
to certain memory ranges based on the initiator’s ID, it is rather a safety
feature. _Runtime Manipulation Detection_ takes a step further towards the
protection of the integrity of runtime memory. However, currently measures
mainly rely on reducing vulnerabilities in ECU software by emphasizing secure
coding practice (e.g. MISRA C/C++ coding standard) and extensive testing (e.g.
fuzzing). Many ECUs can be flashed with signed software. However, legacy ECUs
or those with limited capacities often rely on a tester or an OTA master for
signatures verification instead of doing it inside the target ECU, which
increases the risk of exploitation of race conditions.

Network security iPhones and cars use very similar network security
technologies for secure communication. TLS, VPN, or the latest Wi-Fi or
Bluetooth security standards are commonly implemented in both worlds. Apple
randomizes the addresses and identifiers in wireless communication to prevent
device tracking. It uses application firewall to filter communication on a
per-application basis. Network security measures in cars are often allocated
at entry point. For in-vehicle network, commonly a part of the Ethernet is
secured by TLS and a handful of bus signals are protected by Secure Onboard
Communication (SecOC) for integrity and freshness. A gateway ECU implements
firewall and packet/message filtering on cross-domain communications.

##  3\. Conclusion and discussion

The comparison shows that the “Apple” is several years ahead of the “cars”.
Despite being in different product categories, the security model of an iPhone
from today should be the goal of the cars of tomorrow. The security of iPhone
characterizes a zero trust model from the silicon level up, while cars still
follow a network-based perimeter model. As the car industry evolves towards
high performance vehicle computers and zone architecture, the security model
needs to be adapted to couple with increased system integration and function
mix. iPhones demonstrate that it is technologically feasible to reach a high
level of security without sacrificing performance and usability. Despite the
effort, iPhone security measures are still bypassed or defeated on a rare but
constant basis. Some of the attack methods will be applicable to automotive
embedded systems. Recognizing the evolving threat model will give the
automotive industry incentive for raising the bar on security for the cars of
tomorrow.

## References

  * (1) Apple Platform Security (May 2021), https://support.apple.com/guide/security/welcome/web