The Thickness of Galaxy Disks from z = 5 to 0 Probed by JWST

Although a thick disk is a structure prevalent in local disk galaxies and also present in our home Galaxy, its formation and evolution are still unclear. Whether the thick disk is born thick and/or gradually heated to be thick after formation is under debate. To disentangle these two scenarios, one effective approach is to inspect the thickness of young disk galaxies in the high-redshift Universe. In this work we study the vertical structure of 191 edge-on galaxies spanning redshift from 0.2 to 5 using JWST NIRCAM imaging observations. For each galaxy, we retrieve the vertical surface brightness profile at 1 R e and fit a sech2 function that has been convolved with the line spread function. The obtained scale height of galaxies at z > 1.5 shows no clear dependence on redshift, with a median value in remarkable agreement with that of the Milky Way’s thick disk. This suggests that local thick disks are already thick when they were formed in early times and secular heating is unlikely to be the main driver of thick disk formation. For galaxies at z < 1.5, however, the disk scale height decreases systematically toward lower redshift, with low-redshift galaxies having comparable scale height with that of the Milky Way’s thin disk. This cosmic evolution of disk thickness favors an upside-down formation scenario of galaxy disks.


INTRODUCTION
The thick disk was first discovered in external galaxies Burstein (1979); Tsikoudi (1979), and then also identified in the Milky Way by counting stars in the vertical direction at the solar radius (Gilmore & Reid 1983;Jurić et al. 2008).Such geometric thick component was later on found to be prevalent in local disk galaxies (Yoachim & Dalcanton 2006;Comerón et al. 2012Comerón et al. , 2018)), albeit suspicion that it may be artificial feature caused by the scatter light (de Jong 2008; Sandin 2014).With a careful treatment of the extended wings of the point spread function (PSF), Comerón et al. (2018) confirmed the prevalence of thick disk in local edge-on galaxies.
While the thick disk was first discovered in external galaxies, it is more extensively studied in the Milky Way.Early studies based on solar neighborhood observations have shown that the geometric thick disk of the Milky Way hosts stellar populations that are generally old, chemically metal-poor and α-enhanced, and kinematically hot (e.g., Fuhrmann 1998;Bensby et al. 2005;Reddy et al. 2006;Lee et al. 2011;Adibekyan et al. 2012;Haywood et al. 2013;Bensby et al. 2014).More recent massive stellar spectroscopic surveys, however, reveal more complexity of the geometric thick disk which may be composed of different stellar populations at different radii (Martig et al. 2016;Nidever et al. 2014;Hayden et al. 2015;Lian et al. 2020;Queiroz et al. 2020;Lian et al. 2022).In particular, the thick disk at larger radii tend to contain stars that are relative younger and less α-enhanced.For the thick disks in local galaxies, their stellar population properties have not been extensively studied yet.Based on integral field unit observations of a few nearby edge-on galaxies, the local extra-galactic thick disks are likely similar to that of the inner Milky Way in both stellar population and kinematic properties (Comerón et al. 2016(Comerón et al. , 2019;;Pinna et al. 2019).
One popular explanation of the thick disk formation assumes that the thick disk stars observed today were initially formed within a thin disk and then vertically heated to form a thick disk via interactions with disk symmetries (e.g., spiral arms and Giant Molecular Clouds, Villumsen 1985) and/or galaxy companions (Quinn et al. 1993;Wyse et al. 2006;Villalobos & Helmi 2008).In this scenario, stars born earlier have more chance to be heated and therefore are located in a thicker component, predicting a positive correlation between the age and the vertical extension of stellar populations.On the contrary, another scenario that was proposed more recently assumes the local thick disks were initially formed to be thick (Brook et al. 2004;Bournaud et al. 2009).This scenario is supported by the fact that galaxies in the early Universe, when the local thick disk stars are expected to form, are much more turbulent that the local galaxies given their large velocity dispersion (e.g., Weiner et al. 2006;Shapiro et al. 2008;Simons et al. 2016;Übler et al. 2019) and disturbed and clumpy morphology (e.g., Elmegreen et al. 2005Elmegreen et al. , 2009)).
In galaxy simulations, a mixture of results has yet been reported.In some simulations, the stars in a thick disk today are initially born in a thin component which is then gradually thickened by vertical heating after birth (Meng & Gnedin 2021;Yi et al. 2023).In other simulations (Clarke et al. 2019;Bird et al. 2021), however, the stars already form a thick disk in the early times, but are also then vertically heated afterwards to form a thicker disk.Whether the local thick disks are born already thick or heated to be thick after formation is still under debate.
While both scenarios can broadly explain the old nature of the thick disks in the local galaxies and the Milky Way, the composite stellar populations of the Milky Way's thick disk revealed by recent massive stellar spectroscopic surveys impose challenges to the heating scenario.In particular, Lian et al. (2022) reported the thickness of mono-abundance populations of the Milky Way and found signature of two disconnected phases of thick disk formation separated by a few Gyrs.These results suggest a non-monotonic thickening history of the Milky Way and thus the vertical heating is not likely the major mechanism for the Milky Way's thick disk formation.
To understand the thick disk formation of galaxies in general, studying the disk thickness of high-redshift edge-on galaxies provides an alternative and intriguing perspective.A few early attempts have been made in this direction using HST observations (Elmegreen & Elmegreen 2006;Elmegreen et al. 2017;Hamilton-Campos et al. 2023).However, even with HST spatial resolution, high-redshift edge-on galaxies are marginally resolved.
In this work, we study the thickness of edge-on galaxies beyong the local Universe using the latest JWST observations, which have higher spatial resolution and are deeper than the HST observations that allow us to study the vertical structure of more distant galaxies.Throughout this paper, we assume the cosmological parameters (Ω M , Ω Λ , h 0 ) = (0.27, 0.73, 0.71) and adopt AB magnitude system.

Sample selection
In this work, we select edge-on galaxies using public JWST/NIRCAM imaging data from Cosmic Evolution Early Release Science (CEERS) survey.CEERS survey (PI: S. Finkelstein) is designed to cover 100 arcmin 2 within Extended Groth Strip (EGS) field with JWST infrared imaging and spectroscopy using NIRCam, MIRI, and NIRSpec (Finkelstein et al. 2023).We use the JWST imaging data that are reduced with grizli1 and msaexp (Valentino et al. 2023) and released in the DAWN JWST Archive (DJA), which is a repository of public JWST galaxy data.
We focus on the JWST observations at one single band, F115W, which among all JWST filters has currently the longest exposure time in the CEERS survey and is a filter at short wavelength with relatively sharp PSF.In the released F115W images in the CEERS survey, a total number of 98149 sources are detected.For each source, we retrieve a 50×50 pixel (∼ 2 ′′ × 2 ′′ ) image cut and obtain basic shape and photometry parameters (e.g., minor-to-major axis ratio, position angle) using Sextractor.We then measure the effective radii of these galaxies by extracting light profile using a series of ellipses with the same shape parameters derived above.Since the JWST CEERS field overlaps with the HST EGS field, we adopt the physical properties of these galaxies, including photometric/spectroscopic redshifts, stellar mass, and star formation rate (SFR) from public EGS catalogs (i.e., mass and redshift from Stefanon et al. (2017) and SFR from Barro et al. (2019)).
We then apply the following criteria to select our sample of edge-on galaxies: • edge-on with axes ratio < 0.4; • relatively massive with stellar mass log(M * /M ⊙ ) > 8.5; • star forming with specific star formation rate log(sSFR/yr −1 ) > −11; • redshift between 0 and 5.
A number of 352 galaxies satisfy these criteria.By visual inspection, we further remove 161 galaxies which are either wings of bright stars, part of nearby spiral galaxies, or close to companion galaxies or bright stars.Finally, the remaining 191 galaxies consist of our edge-on galaxy sample.For three of them, spectroscopic redshift is available.Figure 1 shows the distribution of our edge-on galaxy sample in stellar mass and redshift.The majority of the sample have a stellar mass below 10 9.5 M ⊙ and redshift between 1 and 3.

Construction of the line spread function
To account for the effect of the instrument broadening on the light distribution in the vertical direction, we derive the line spread function (LSF) following the approach introduced in Elmegreen et al. ( 2017) and Comerón et al. (2018).We first generate NIRCAM 2D PSF model at F115W band oversampled in a 4 × 4 grid with pixel scale of 0.01 ′′ using WebbPSF2 (Perrin et al. 2014).We then convolve this PSF model with an infinitely thin edge-on galaxy that has a radially declining profile.From the mock galaxy image we measure the intensity profile perpendicular to the mid-plane at a certain radius which is assumed to be the line spread function.Figure 2 shows the obtained LSFs for input galaxy with various effective radius (R e ,left panel), measured at different radial position for an input galaxy with R e = 0.2 ′′ (middle panel), and convolved with the 2D PSF model rotated in different angles (right panel).
It can be seen that the LSFs show no clear dependence on the size of the input galaxy and weak dependence on the radial position where the LSF is measured and the orientation of the PSF.In this work we extract the vertical brightness profile of observed galaxies at 1 R e to measure the scale height.Therefore we adopt the LSF measured at 1 R e with input galaxy R e = 0.2 ′′ and non-rotated PSF.We have tested that the scale height measurements with different PSF orientations vary less than 10%.The FWHM of the adopted LSF is 0.039 ′′ , a factor of six smaller than that of the HST images at similar wavelength range in the near-infrared (Hamilton-Campos et al. 2023).This demonstrates the advanced ability of JWST to spatially resolve distant galaxies.

Vertical surface brightness profile (VSBP)
To measure the disk thickness of our edge-on galaxies, we first retrieve the vertical surface brightness profile (VSBP) at 1 R e , which is represented by the median profile between 0.5 and 1.5 R e .The vertical density profiles of galactic disks are well described by sech 2 function (e.g., Xiang et al. 2018;Hamilton-Campos et al. 2023).We therefore fit the VSBP of our galaxies with a sech 2 function Σ = µ 0 × 4 (e ∆z/hz + e −∆z/hz ) 2 that has been convolved with the LSF obtained in §2.2.Here µ 0 and h z are the central surface brightness and the scale height, respectively.The sech 2 function approximates exponential form at large distance, but the scale height in the sech 2 function is one-half of the exponential scale height.
Figure 3 shows the JWST/F115W image and retrieved VSBPs of eight example edge-on galaxies at different redshifts.It can be seen that the observed VSBPs are significantly broader than the LSF, meaning that these galaxies are well resolved by JWST.The best-fitted sech 2 function are shown as orange line of which the scale height is shown at the bottom.
The uncertainties of the scale height measurement is estimated using Mento Carlo simulation.For each galaxy, we first compute a randomised VSBP based on the fiducial profile and associated uncertainties estimated from the weight images that have taken the sky background noise into account.We then redo the fitting for the new VSBP.After repeating this process 100 times, we take the standard deviation of the 100 measurements as the uncertainty of the scale height measurement.This uncertainty estimate, however, is not available for 65 out of 191 galaxies due to large uncertainties in the VSBP that stop the fitting from converging for the randomized VSBP.

Evolution of disk thickness
The obtained scale height of our edge-on galaxies as a function of redshift is shown in Figure 4.The enlarged squares indicate the median trend with redshift bin width of 0.5 at z < 3 and bin width of 1 at z > 3.As shown in Hamilton-Campos et al. (2023), the scale height of galaxies selected with axis ratio < 0.4 would be on average overestimated by 22% because of projection effect in case of non-zero inclination to our line-of-sight.This correction factor has been applied to our scale height measurements as shown in Figure 4.For comparison, we include the scale height measurement of the Milky Way's chemical thick (dashed line) and thin discs (dash-dotted line) from Lian et al. (2022).Since they are derived with an exponential fit, we divide the measurements by a factor of 2 to make them comparable to the sech 2 scale height in this work.At z > 1.5, there is a large scatter of scale height among individual edge-on galaxies with no clear trend with redshift.The median sech 2 scale height is 0.38±0.13kpc, in remarkable consistence with that of the Milky Way's thick disk (0.43 kpc), suggesting that the Galactic thick disk is already thick in the early times.A similarly large scale height of high-redshift galaxies was also found in previous works using HST observations (Elmegreen & Elmegreen 2006;Hamilton-Campos et al. 2023).This result suggests that the secular heating is unlikely the main formation mechanism of thick disk.Interestingly, at z < 1.5, the galaxy disk thickness decreases towards lower redshift, reaching a comparable value of the Milky Way's thin disk at low redshift of z ∼ 0.4.This trend is significant with Pearson correlation coefficient of 0.36.However, no such trend was found before in previous relevant studies based on HST observations (Elmegreen & Elmegreen 2006;Hamilton-Campos et al. 2023).The lack of redshift evolution of disk thickness in these studies is possibly because their distant edge-on galaxies are marginally resolved given the much broader PSF of HST.However, it is worth noting that our edge-on galaxy sample at z < 1 is still small.We expect this could be greatly improved by future wider survey of JWST and forthcoming ground-and space-based wide-field imaging surveys in the near future (e.g., LSST).
The observed redshift evolution of the galaxy disk thickness since z ∼ 1.5 strongly favors an upside-down formation scenario of galactic disks.In this scenario, the early galaxy disk at high redshift is initially born thick because of high   turbulence in the disk that may be caused by, for example, frequent galaxy mergers and interactions, accretion of gas, and star formation feedback.Gas turbulence makes the gas disk thick and therefore the star formation and young stellar disk thick.As galaxies grow and evolve to the low redshift Universe, the gravitational potential well become deeper and merger events become less frequent.As a result, the disk turbulence decreases and star formation occurs in an increasingly thinner disk.This scenario has also been proposed to explain the thick disk formation and evolution of the Milky Way that is revealed in great detail with massive spectroscopic stellar surveys (Bovy et al. 2012(Bovy et al. , 2016;;Mackereth et al. 2017;Yu et al. 2021;Lian et al. 2022).
In contrast to upside-down formation scenario, the secular heating scenario suggests the galaxy disk to be born thin and heated to be thick gradually.These are the opposite of the results shown in Fig. 4 and therefore the heating scenario is unlikely the main driver behind the thick disk formation.To further investigate the role played by secular heating in thick disk formation and evolution, one need to decompose the thin and thick disk of low-redshift galaxies and connect the thick disk of galaxies at different redshifts.
To inspect the mass dependence of galaxy disk thickness, we plot the scale height of our edge-on galaxies against their stellar mass in Figure 5.The sample is split into eight redshift bins to minimize the variation of disk thickness due to redshift.For galaxies at z > 1, there is a tentative positive correlation between the disk thickness and stellar mass.This mass dependence is not clearly seen at low redshift.Different from our result, Elmegreen et al. (2017) found a strong positive correlation between the stellar mass and disk thickness for a sample of 107 edge-on galaxies at z < 1.5 using HST ACS observations.The measurements of disk scale height for all of our galaxies reported here are based on imaging data in a single band of F115W.Therefore for galaxies at different redshifts, the scale height are actually measured in different restframe wavelength range that may trace stellar populations of different ages.Elmegreen & Elmegreen (2006) found that, for the same galaxy, the sech 2 scale height measured based on observations at longer wavelength is systematically higher.To investigate whether and how our results are affected by the choice of the filter for the imaging data, we have measured the scale heights of our galaxies using JWST imaging data in longer wavelength filters, including F150W, F200W, F277W, F356W, and F444W following the same procedure as we adopt in F115W.We find that the obtained individual scale height measurements and redshift trend are remarkably consistent with that obtained in F115W band, especially for galaxies at z > 1.This is understandable because high-redshift galaxies are generally young and thus the stellar populations traced by observations in different wavelength range are not in large difference.We conclude that our results are not dependent on the choice of the filter used for the scale height measurement.

Figure 1 .
Figure 1.The distribution of stellar mass and redshift of our edge-on galaxy sample.
The 1D point spread function (PSF) that is used for modelling the vertical surface brightness profile of edge-on galaxies in our sample.The 1D PSF is constructed by convolving the 2D PSF of JWST F115W images with an infinitesimally thick disk model of a given size.The left panel shows the 1D PSF measured at one effective radii for disk models of different sizes.The right panel illustrates the 1D PSF measured at different position relative to the center for the same disk model with effective radius of 0.2 ′′ .

Figure 3 .
Figure3.JWST/F115W images and vertical surface brightness profiles of eight example galaxies from low redshift at top left to high redshift at bottom right.In the image, the two gray shaded regions cover the pixels within 0.5−−1.5 Re.The length of 1 kpc at each corresponding redshift is illustrated at the bottom right corner of the image.In the VSBP panel, the retrieved VSBP and best-fitted sech 2 model are shown as blue points and orange curve.The obtained sech 2 scale length is presented at the bottom.The line spread function obtained in §2.2 is also included as reference.

Figure 4 .
Figure4.Redshift evolution of galaxy disk thickness.Small blue circles indicate the sech 2 scale height measurements of invidividual edge-on galaxies with only photometric redshift estimate and orange circles for those with spectroscopic redshift.Filled circles are the galaxies with valid uncertainty estimate and empty circles without.The level of uncertainties is shown at the bottom in gray.Enlarged black squares and error bars show the median scale length and 1σ scatter at each redshift bin with bin width of 0.5 at z < 3 and bin width of 1 at z > 3. The converted sech 2 scale height of the Milky Way's thick and thin disks are denoted by the dashed and dash-dotted lines, respectively.

Figure 5 .
Figure5.Scale height as a function of stellar mass for galaxies in eight redshift bins.Each panel shows edge-on galaxies in one redshift bin.The symbols are the same as Figure4.