1 Introduction↩︎

The thick disk was first identified in 1979 as a diffuse stellar component in edge-on S0 galaxies, distinct from both the bulge and thin disk [1], and was subsequently recognized in the Milky Way in 1983 through star counts toward the South Galactic Pole [2]. These discoveries established the now-standard view that the Milky Way disk comprises two components: a thin and a thick disk [3]–[5].

The formation history of the Milky Way’s disk appears to have proceeded in two nearly disjoint phases [6]–[8]. The thin disk, with an exponential scale height of \(\sim\)​0.3–0.4 kpc [5], [9], is dominated by younger stars and was formed from the cold gaseous mid-plane during the phase of steady star formation over the past \(\sim\)​8 Gyr [7], [10], [11]. The thick disk, with an exponential scale height of \(\sim\)​0.7–1.2 kpc, consists primarily of old, [\(\alpha\)/Fe]-enhanced stars that originated during an earlier, bursty star formation phase within the first \(\sim\)​5 Gyr [4], [5], [9], [11], [12]. Thick disks are now known to be nearly ubiquitous among nearby disk galaxies [13]–[19].

Despite their prevalence, the physical origin of thick disks remains debated. They may have been born as inherently thick structures under turbulent, gas-rich conditions in the early Universe [11], [20]. Alternative, they may have grown from initially thin disks through heating by minor mergers or internal secular processes [21]–[24].

Studying galaxy disks across cosmic time provides crucial insight into the mechanisms that construct their vertical structure. In the HST era, the limited spatial resolution made it difficult to distinguish thin and thick disks, restricting studies to global scale-height measurements. Tracing rest-frame optical to UV wavelength, [25] reported a median global scale height of \(0.63\) kpc, with a scatter of \(0.24\,\)kpc among individual measurements, at \(z\sim2\) (see also [26]) and detected vertical color gradients that suggest the coexistence of thin and thick components at early epochs. A later analysis by [27] yielded a median scale height of \(0.74\) kpc, with a scatter of \(0.35\,\)kpc, after correcting for the overestimation caused by projection effect of deviations from perfect edge-on inclination; however, such corrections depend essentially on the ratio of scale height to scale length, a dependence explored in this study.

Progress in characterizing disk structure beyond \(z\sim1\) has recently accelerated thanks to the unprecedented sensitivity and resolution of JWST’s Near-Infrared Camera (NIRCam). JWST observations have revealed a substantial population of regular disk galaxies already in place within the first few Gyrs, suggesting that disk assembly and evolution began earlier than previously thought [28]–[38]. Using F115W imaging that traces rest-frame NIR at \(z\sim0\) but UV at \(z>3\) for a sample with a median stellar mass of \(10^9\,M_\odot\), [39] measured a median global scale height of \(0.38\) kpc, with a scatter of \(0.13\) kpc, at \(z=0.2\)–5, lower than the previous HST results. [40]¹, using F277W, F356W, and F444W filter to trace rest-frame 1–2 wavelength at \(z\leq 3\), obtained a median global scale height of \(0.66\pm0.07\) kpc (where all quoted values in this work follow the +/\(-\) notation for uncertainties) for stellar mass of \(10^{9.2}\,M_\odot\), consistent with earlier HST results. Moreover, they performed the first thick-thin disk decomposition at \(z>1\), finding median scale heights of \(0.33\pm0.04\) kpc and \(0.94\pm0.07\) kpc for the thin and thick components, respectively.

Still, as these pioneering JWST studies do not adopt a fixed rest-frame wavelength for the measurements, the derived scale heights may be biased by extinction from mid-plane dust lanes or by spatial variations in stellar populations. Furthermore, their analyses were based on relatively small samples (fewer than 200 galaxies). Expanding to larger, statistically robust samples is essential to establish the redshift evolution of disk scale heights. Equally important is correcting for the bias introduced by variation in rest-frame wavelength and developing a more accurate method to correct for the bias introduced by deviations from a perfectly edge-on inclination. Motivated by these needs, we present a comprehensive analysis of edge-on disk galaxies observed by the COSMOS-Web survey [41]. In this work, we focus on the global disk scale height, averaging thin and thick components when present, while future studies will address thin–thick disk decompositions.

This paper is organized as follows. Section 2 describes the sample selection and its properties. Section 3 details the image analysis and scale height measurements. Section 4 presents the main results. Section 5 discusses the implications for the formation and evolution of galaxy disks. Finally, Section 6 summarizes the main findings. Throughout this paper, we adopt AB magnitudes, [42] initial mass function, and assume a flat \(\Lambda\)CDM cosmology with \((\Omega_{\rm m}, \Omega_{\Lambda}) = (0.27, 0.73)\) and a Hubble constant of \(H_0 = 70\,{\rm km\,s^{-1}\,Mpc^{-1}}\).

2 Sample and Data↩︎

The COSMOS-Web survey is the largest galaxy survey conducted in JWST Cycle-1 [41]. Its NIRCam observations in four filters (F115W, F150W, F277W, and F444W) provide a deep and high-resolution dataset for our structural and photometric analysis. For this study, we use the first data release of COSMOS-Web NIRCam mosaics² [44]. The mosaics have a pixel scale of \(0\farcs03\) at all filters. [45] constructed the COSMOS-Web photometric catalog, by performing source extraction across 33 photometric bands spanning 0.3–8 \(\micron\) and measured photometry across available bands using a multi-band model-fitting approach with SourceX-Tractor++ [46], [47]. Using these measurements, [45] derived physical properties for each source through template spectral energy distribution fitting using LePhare [48], [49]. From the COSMOS-Web catalog (version 1.1), we obtain photometric redshifts (\(z\)), stellar masses (\(M_*\)), and absolute magnitudes in the Near-Ultraviolet (NUV), \(r\)-, and \(J\)-bands (\(M_{\rm NUV}\), \(M_r\), and \(M_J\)). Figure 1 illustrates the rest-frame wavelength probed by each NIRCam filter as a function of redshift. In particular, the F444W filter observes rest-frame 1 \(\micron\) at \(z=3.5\), above which the rest-frame pivot wavelength drops below 1 \(\micron\). To minimize the impact of dust extinction on our morphological measurements, we require that the reddest filter probes at least 1 \(\micron\). We therefore limit our analysis to galaxies within \(0 < z < 3.5\). The choice of the upper redshift limit also aligns with the epoch of the Milky Way’s thick disk formation [50].

We further apply the following selection criteria to construct our parent sample. To ensure a clean and reliable galaxy sample, we require the following flags provided in the COSMOS-Web catalog: type = \(0\) (select galaxies instead of stars), flag_star_hsc = \(0\) (exclude sources severely contaminated by nearby stars), and warn_flag = \(0\) (remove spurious detections such as hot pixels, sources detected in only a single filter, or other imaging artifacts). We restrict our analysis to galaxies with \(M_*\geq10^{10}\,M_{\odot}\), as recent studies suggest that lower-mass galaxies at \(z>1\) may be prolate rather than disky systems [51]–[53]. Applying these criteria yields a parent sample of 19199 galaxies.

We visually inspected each galaxy in the parent sample and excluded 1445 objects that were under interaction/merger or were too close to other galaxies or bright stars for reliable morphology measurements. This step results in a cleaned sample of 17754 galaxies. For each object, we use the Python package SEP [54], [55] on the F444W images to generate a target galaxy segmentation map and generate a mask file for removing contaminating sources. Mask generation was performed with both cold and hot detection modes, with parameters adjusted manually. These masks are used for subsequent model fitting.

Figure 2 shows \(M_{\rm NUV}-M_r\) vs. \(M_r-M_J\) color-color diagrams for four redshift bins. The background contours indicate the distribution of the parent sample, with contours enclosing 10%, 20%, ..., 90% of galaxies. Quenched galaxies were identified using the division line (black line in Figure 2) proposed by [43] and were excluded, leaving 13453 star-forming galaxies for analysis. We first perform single Sérsic fitting using the code IMFIT [56] and remove unresolved sources with effective radius in F444W band less than F444W PSF FWHM. We then perform bulge-disk decomposition, modeling each galaxy with an exponential disk and a round Sérsic bulge with the Sérsic index fixed at 4, in all available filters. The generation of point-spread functions (PSFs) used in the fitting is described in Section 3.1. Catastrophic fitting results are excluded. From this decomposition, we calculate the bulge-to-total light ratio (\(B/T\)) and exclude potential elliptical galaxies, either introduced by photometric measurement uncertainties or representing star-forming ellipticals, by adopting a threshold of \(B/T>0.5\). Using a higher threshold (e.g., \(B/T=0.7\)) does not affect our main results (see Section 4.3). We measure the axis ratio of the disk component (\(q = b/a\), where \(a\) and \(b\) are the semi-major and semi-minor axes, respectively) in each filter. Using interpolation, we derive the disk axis ratio at a consistent rest-frame wavelength of 1 \(\micron\) for all galaxies, which we define as the final disk axis ratio (\(q_d\)). Edge-on disks appear as flattened systems and can be efficiently selected using \(q_d\). This criterion is objective and reproducible, in contrast to visual classification. We therefore select our nearly edge-on disk sample by requiring \(q_d < 0.4\), following the criterion adopted by [27], which is slightly more permissive than the \(q_d < 0.3\) threshold used by [25]. As demonstrated in Section 4.2, our method for correcting biases introduced by projection effects is effective, such that adopting either \(q_d < 0.3\) or \(q_d < 0.4\) yields nearly identical bias-corrected scale heights. The above criteria reduce the 13453 galaxies to 2631 edge-on disks, the sample studied in this work. This sample is 15–30 times larger than in previous studies.

In contrast to our PSF-corrected measurements, applying an axis-ratio cut based on values derived without PSF correction (such as those from segmentation analysis), as done in some previous studies, can bias the sample toward intrinsically flatter systems at higher redshifts due to the smaller apparent sizes of galaxies. Our edge-on disk sample is shown in Figure 2 in the color-color diagram (blue dots), with median stellar mass values (\(\sim2.5\times10^{10}\,M_\odot\)) and scatter indicated at the bottom. The stellar masses of our sample lie well above the completeness limits at each redshift [45].

We note that there are apparent clustering of galaxies at some specific redshift is a known artifact of the spectral template fitting process for determining the photometric redshift. This can occur when prominent spectral features, such as the Balmer or 4000 Åbreak, transition between photometric bands. We have verified that the clustering is observed for both edge-on and face-on galaxies, indicating that the selection of edge-on disks does not introduce any unusual biases caused by the redshift determination. The overall quality of the photometric redshifts in our sample is high, with a catastrophic failure rate of only 1.44% and a normalized median absolute deviation of 0.011 compared to spectroscopic redshifts [45].

3 Measurements↩︎

3.1 Hybrid PSFs↩︎

Accurate PSFs are essential for robust morphological fitting of galaxies. Theoretical PSFs simulated with WebbPSF [57] tend to be narrower than those derived from real stars in drizzled mosaics [58], [59], as WebbPSF is designed to simulate the PSF on a single exposure, and should not be expected to agree with the PSF from drizzled mosaics, which include additional pixel-level convolutions. Since our analysis is based on drizzled mosaics, we construct empirical PSFs (ePSFs) from field-star images in each NIRCam filter using PSFEx [60], as shown in the top row of Figure 3. Still, these ePSFs suffer from relatively low signal-to-noise ratios (\(S/N\)) in their outer regions. A common strategy to improve PSF quality is to combine the high-\(S/N\) core of an empirical PSF with the noise-free wings of a theoretical one, forming a hybrid PSF [61], [62]. Following this approach, we generate theoretical PSFs for each NIRCam filter using WebbPSF³ and rotate them to match the orientation angle of the corresponding ePSFs by minimizing their difference. We then replace the noisy pixels (approximately the faintest 10% in flux) of the ePSF with values from the theoretical PSF. A smoothed step function and reversed step function are applied during this combination to ensure a seamless transition. Finally, we normalize the combined image to produce the hybrid PSF, as illustrated in the bottom row of Figure 3. These hybrid PSFs, which preserve the realistic core structure of observed stars while maintaining accurate outer profiles, are adopted for all model fittings throughout this work.

3.2 Scale-height measurement↩︎

The scale height provides a quantitative measure of disk thickness, with larger values corresponding to thicker disks. In a coordinate system aligned with the edge-on inclination, we define the radial distance along the disk plane (parallel to the major axis) as \(R\), and the vertical distance from the plane (parallel to the minor axis) as \(h\). The galaxy surface brightness at \((R, h)\) is denoted as \(I(R,h)\). Assuming an isothermal disk, the light distribution of an edge-on disk can be described by [63]: \[\label{Eq:general} I(R,h) = I_0 \cdot (R/R_d) \cdot K_1(R/R_d) \cdot {\rm sech}^{2/m}(mh/h_0),\tag{1}\]

where \(I_0\) is the central intensity, \(K_1\) is the modified Bessel function of the second kind of order one, \(R_d\) is the exponential scale length in the disk plane, \(h_0\) is the scale height, and \(m\) controls the overall shape of the vertical profile. Because \(m\) and \(h_0\) are degenerate, \(m\) is usually fixed to 1 or \(+\infty\) when determining the scale height, and Equation (1 ) then simplifies to the following two common forms: a sech\(^2\) function, \[\label{Eq:sech295} I \propto {\rm sech}^{2}(h/h_0),\tag{2}\] and an exponential function, \[\label{Eq:exp95} I \propto \exp(-2h/h_0 ).\tag{3}\] Due to the above functional forms, it is generally approximated that the sech\(^2\) scale height is twice the exponential scale height [13]. However, this approximation is not exact, as the vertical profile of a galaxy is independent of the functional form adopted. [5] suggest the sech\(^2\) scale height is 1.8 times the exponential scale height, whereas our analysis (see below) indicates that a factor of 1.37 provides a more accurate conversion for our sample.

Two-dimensional (2D) model fitting is an efficient approach for measuring scale heights [13], [40], [64]. We use IMFIT [56] to perform the 2D fittings, incorporating an edge-on disk model and a bulge model. The hybrid PSFs generated in Section 3.1 are used. The bulge is modeled with a round Sérsic function [65] with \(n=4\): \[\label{eq:bulge} I(R) = I_e \exp\left\{ -b_n \left[ \left( R/R_e \right)^{1/n} - 1 \right] \right\},\tag{4}\] where \(R_e\) is the half-light radius, \(I_e\) is the intensity at \(R_e\), and \(b_n\) is a function of \(n\) determined by the incomplete gamma function.

We adopt the sech\(^2\) vertical profile and the exponential radial profile for the edge-on disk model⁴: \[\label{Eq:sech2} I(R,h) = I_0 \cdot (R/R_d) \cdot K_1(R/R_d) \cdot {\rm sech}^{2}(h/h_0),\tag{5}\] where \(h_0\) is the sech\(^2\) scale height. The Sérsic-sech\(^2\) decomposition is carried out for each galaxy in all available filters. The galaxy center is a free parameter to fit, as this can mitigate the effects of dust lanes (see Section 4.1). The uncertainties of the best-fit parameters are provided by IMFIT as the statistical uncertainties obtained from evaluating the structure of the \(\chi^2\) minimum. The F115W and F150W images generally have lower \(S/N\) than the F277W and F444W images, as the stellar emission of galaxies typically peaks at the NIR. At high redshift, some galaxies observed in F115W or F150W are indistinguishable from, or only marginally above, the background noise. To ensure robust measurements, we exclude any filter whose mean flux within the galaxy segmentation falls below the background noise level and instead adopt filters with sufficiently high \(S/N\). The scale height derived from the Sérsic-sech\(^2\) fitting represents the scale height averaging thin and thick components when present. Figure 4 presents representative fitting results for galaxies at \(z=0.54\) and \(z=2.45\). We further interpolate the measurements to a fixed rest-frame wavelength of 1 to mitigate the effects of dust extinction and spatial variation in stellar populations, as detailed in Section 4.1. The correction for projection biases arising from deviations from perfect edge-on orientation is described in Section 4.2.

For comparison, we also measure the exponential scale height (\(h_{0,\rm{exp} }\)). To minimize degeneracy between components, we fix the bulge parameters to those obtained from the Sérsic-sech\(^2\) fit and then refit the image using an exponential vertical profile: \[\label{Eq:exp} I(R,h) = I_0 \cdot (R/R_d) \cdot K_1(R/R_d) \cdot \exp(-h/h_{0,\rm{exp}} ).\tag{6}\] We use \(h_{0,\rm{exp} }\) instead of \(h_{0,\rm{exp}}/2\) in Eq. (3 ), as \(h_{0,\rm{exp} }\) is more commonly adopted in the literature. Figure 5 presents the correlation between the two measurements, which follows a tight linear relation described by \[\label{exp95sech2} {\rm sech}^2~h_0 = 1.37\times {\rm exponential}~h_{0,\rm{exp}}.\tag{7}\] The Pearson correlation coefficient between sech\(^2\) \(h_0\) and exponential \(h_{0,\rm{exp}}\) is \(r_{\rm p}=0.99\), indicating a statistically significant correlation (\(p\)-value \(< 0.01\)). Therefore, our results remain unaffected by whether the vertical profile is described using a sech\(^2\) or an exponential function. In the rest of this work, the term “scale height” denotes the sech\(^2\) scale height unless the exponential scale height is specified.

In addition to the 2D method, one-dimensional (1D) fitting has also been used to measure scale heights [25], [27], [64]. The 1D method is able to obtain scale height as a function of radius. We evaluate the consistency between \(h_0\) derived using 1D and 2D approaches. For the 1D method, we first subtract the best-fit bulge component (if present) and rotate the galaxy so that its major axis aligns with the image \(x\)-axis. We then extract vertical intensity profiles along the galaxy major axis. Profiles with peak intensities below three times the background noise are excluded. Each vertical profile is fitted with a 1D \({\rm sech}^2\) function, convolved with the line-spread function (LSF), from which we derive the radial profile of scale height. The LSF is constructed by convolving a one-pixel line with a PSF following the strategy in [25]. The profiles are stacked to create a composite profile. A final \({\rm sech}^2\) fit to the stacked profile yields the average \(h_0\) from the 1D method, with weighting similar to that of the 2D fitting.

A comparison of the 1D and 2D measurements is shown in Figure 6. The median difference between the two is 0.04 kpc, with a standard deviation of 0.05 kpc. These small values indicate that the scale heights derived from the 1D and 2D methods are consistent when similar weighting schemes are applied. The radial variation of scale height will be discussed further in Section 4.

4 Results↩︎

4.1 Wavelength dependence of scale height and the correction↩︎

Since younger stellar populations dominate the thin disks of galaxies, observations at shorter rest-frame wavelengths without dust extinction are expected to reveal thinner disks with smaller scale heights. However, for edge-on systems, particularly the actively star-forming disks in our sample, dust lanes along the mid-plane substantially attenuate the flux especially at shorter wavelengths. This obscuration impacts the apparent vertical light profile and can overestimate the measured scale height [14]. As described by Eq. (5 ), the scale height is defined as the height at which the intensity falls to sech\(^2\)(1) of the central value, rather than as an absolute geometric thickness. When the mid-plane flux is suppressed by a dust lane, the observed profile appears artificially broadened, leading to an overestimation of the intrinsic scale height.

In Figure 7, we compare the vertical scale heights (\(h_0\)) measured in the F444W band with those measured in F277W, F150W, and F115W. The scale heights derived from the bluer filters are systematically larger than those from redder filters. Quantitatively, the \(h_0\) measured in F277W, F150W, and F115W exceeds those in F444W by \(0.12\pm0.11\) kpc, \(0.29\pm0.28\) kpc, and \(0.30\pm0.29\) kpc, respectively.

To visualize how dust attenuation overestimates the measured scale height, we plot in Figure 8 examples of the 1D vertical profile fitting (see Section 3.2 for the method) for a galaxy observed in F444W (top) and F115W (bottom). The first column shows the bulge-subtracted images, where the prominent dust lane in F115W is clearly visible. The second column presents the radial profile of \(h_0\), showing systematically larger values in F115W than in F444W, consistent with the 2D fitting results for the same galaxy shown in Figure 4. The third column shows the corresponding LSFs, while the final column compares the vertical light profiles. The F444W vertical light profile is symmetric and smooth, whereas the F115W profile exhibits a flux depression caused by the dust lane near the center. This depression shifts the apparent center and broadens the observed distribution, producing an overestimated scale height. When all fitting parameters are allowed to vary freely, we obtain \(h_0=0.71\) kpc in F444W and \(h_0=0.99\) kpc in F115W. If the galaxy center in F115W is fixed to that of F444W, the fitted value in F115W increases further to \(1.37\) kpc. Setting the center free for the fitting can mitigate the effect caused by dust lane, which motivates our choice of leaving the center unconstrained in the bulge-edge-on-disk decomposition (Section 3). In addition to dust extinction, spatial variations in the stellar population, such as younger stars concentrated near the mid-plane and older stars located at larger vertical distances, can in principle also affect the scale height measured at different wavelengths, although this effect is less significant than dust extinction.

Figure 9 presents the dependence of the measured \(h_0\)on rest-frame wavelength, with results shown separately for \(0<z\leq1\) (top) and \(1<z<3.5\) (bottom). The median values and their uncertainties reveal a clear trend: scale heights systematically decrease toward longer rest-frame wavelengths, where dust extinction becomes weaker. This dependence on rest-frame wavelength should have also incorporated the effects of spatial variations in the stellar population. The measured \(h_0\) at rest-frame 5000 Åis overestimated by roughly 30–40%, whereas at rest-frame \(1\,\micron\), the impact of extinction is small.

This wavelength dependence therefore acts as a confounding factor when a single filter is used to study the redshift evolution of scale height, if the filter traces rest-frame optical to UV wavelengths at high redshift but NIR wavelengths at low redshift. To isolate the true redshift evolution, we determine \(h_0\) for each galaxy at a fixed rest-frame wavelength of \(1\,\micron\) (\(h_{0,1\micron}\)) by fitting a power law to the \(h_0\) versus rest-frame wavelength, using the uncertainties from 2D fitting as weights. A similar wavelength dependence affects the disk scale length, \(R_d\), owing to the inside-out growth of galaxies [66], [67]. To ensure consistency, we likewise derive \(R_d\) at 1 \(\micron\), \(R_{d,1\micron}\). The effective radius of the disk component at 1 (\(R_{e,\rm{disk},1\micron}\)) is converted from \(R_{d,1\micron}\).

We adopt \(1\,\micron\) as the reference wavelength as extinction effects are largely mitigated at this wavelength, and it remains accessible at high redshift with the F444W filter. At longer rest-frame wavelengths (\(\gtrsim1.5\,\micron\)), \(h_0\) gradually converges to values \(\sim\)​10% lower than \(h_{0,1\micron}\) at both low and high redshifts. Nevertheless, by anchoring both \(h_0\) and \(R_d\) to the same rest-frame wavelength, our analysis mitigates wavelength-dependent effects, including those caused by dust extinction and stellar population effects, and enables a clean and uniform comparison of disk structure across redshift. The measurements of scale heights and scale lengths can be accessed at Zenodo [68].

4.2 Projection effects of deviations from a perfect edge-on inclination and the bias correction↩︎

A randomly oriented three-dimensional (3D) galaxy disk is unlikely to be perfectly edge-on. A deviation from perfect edge-on inclination (\(i = 90°\)) amplifies the scale height measured in a 2D projection image, which is done by assuming a perfect edge-on inclination. This constitutes a projection effect overestimating the scale height. Consequently, correcting for this effect is essential to derive the intrinsic scale heights and their intrinsic cosmic evolution, particularly for our sample, which is selected based on a cut in the observed disk axis ratio \(q_d<0.4\).

A constant correction factor of 0.78 has been proposed to address this bias [27]. However, this correction is not accurate enough, because the influence of the deviation from a perfect edge-on inclination on scale height measurement actually depends on the relative scale height of the disk, demonstrated below. The relative scale height (or relative thickness) is defined as \[{\rm relative~scale~height}=h_0/R_d.\] We generate simulated projected 2D images of disks at specific inclination by using IMFIT [56], which performs line-of-sight integration of 3D disk models with exponential radial and \({\rm sech}^2\) vertical profiles. For illustration purpose, we consider two disks with \(h_0/R_d=0.4\) and \(0.2\), respectively, and rotate them at \(i=90°\), \(85°\), \(80°\), and \(75°\). The projected 2D images with intrinsic and observed values of the disk parameters are shown in Figure 10. For the model with \(h_0/R_d=0.4\), when inclination is set from \(i=90°\) to \(i=75°\), \(h_0\) is amplified from 0.5 kpc to 0.8 kpc, and \(q_d\) is amplified from \(0.246\) to \(0.389\). But for \(h_0/R_d=0.2\), \(h_0\) is amplified from 0.5 kpc to 1.25 kpc and \(q_d\) is amplified from \(0.123\) to \(0.307\). We can learn two key points from this test: (1) the scale height of a disk with lower relative scale height (e.g., lower \(h_0\) at fixed \(R_d\), or larger \(R_d\) at fixed \(h_0\)) suffers more severely from the projection effect; (2) a disk with lower relative scale height is less likely to meet the cut in \(q_d\) when selecting a sample.

Next we quantify the impact of this projection effect and derive a new bias correction method. We use Imfit to construct a suite of disk models spanning a range of \(h_0/R_d\) and inclination angles. For each model, we perform fits to measure \(h_0\), \(R_d\), and \(q_d\). The measurement of \(R_d\) remains robust, but \(h_0\) is overestimated. To derive the correction for measured \(h_0\), we apply our sample selection criterion by excluding models with \(q_d \ge 0.4\) and compute the median inclination angle and median observed \(h_0\), weighted by the inclination probability distribution \(P(i) \propto \sin(i)\).

The results are presented in Figure 11. The top-left panel illustrates the median observed value of \(h_0\)of disks with larger intrinsic \(h_0/R_d\)are less overestimated. The top-middle panel demonstrates that disks with larger intrinsic \(h_0/R_d\)ratios can have smaller maximum deviations from \(i = 90°\) before exceeding the \(q_d\) threshold, simply because relative thinner disk has smaller axis ratio. The top-right panel present the median inclination angle as a function of intrinsic \(h_0/R_d\). The bottom-left panel show the conversion between intrinsic and median observed \(h_0/R_d\). The bottom-middle and bottom-right panel shows the derived statistical correction factor \(f_{\rm corr}\), defined as: \[f_{\rm corr} = \frac{{\rm median~intrinsic}~h_0}{{\rm median~observed}~h_0},\] plotted as a function of intrinsic and median observed \(h_0/R_d\), respectively. Disks with larger \(h_0/R_d\)ratios exhibit correction factors closer to 1, because (1) their \(h_0\)are less affected for a given inclination compared to those with lower intrinsic \(h_0/R_d\)(see Figure 10) and (2) their maximum deviation from a perfect edge-on inclination before meeting the cut in \(q_d\) is smaller (top-middle panel). The variation in \(f_{\rm corr}\) highlights the importance of the newly developed correction method, as a redshift evolution of \(h_0/R_d\) indeed exists, as shown in the next section. Table [tab] shows the parameters for deriving the correction factors in samples selecting using the axis ratio cut \(q_d<0.4\) as in this work. Table [tab9503] provides values for samples selecting using \(q_d<0.3\). If a smaller critical axis ratio is used, the \(f_{\rm corr}\) is closer to unity.

The application of our bias correction is straightforward: we first compute the median observed \(h_0/R_d\)and determine the corresponding \(f_{\rm corr}\), and then obtain the corrected median \(h_0\)by multiplying the observed value by \(f_{\rm corr}\). As described above, the correction curve is derived for a given intrinsic \(h_0/R_d\). In real observations, however, disks exhibit a distribution of \(h_0/R_d\) rather than a single value. To assess the effectiveness of our method under real observational conditions, we create four groups of models, each containing 10,000 models with random inclinations. The distributions of \(R_d\), \(h_0\), and \(h_0/R_d\)for these four groups are shown in Figure 12. A truncation is observed at \(h_0/R_d=0.65\) in the \(h_0/R_d\)distribution, as models with \(h_0/R_d\)above this threshold exhibit an axis ratio \(q_d \geq 0.4\) even in a perfectly edge-on inclination. For each inclined model, the observed \(h_0\) and \(h_0/R_d\)are measured using 2D \({\rm sech}^2\) fitting. For each group, we calculate the median observed \(h_0\) and \(h_0/R_d\), and apply our new correction method to obtained the bias-corrected \(h_0\).

Figure 13 presents a comparison of the median measured \(h_0\), with and without correction, against the intrinsic values. The results without correction are shown as crosses. As expected, the observed \(h_0\) values are systematically overestimated compared to their intrinsic values when no correction is applied. Furthermore, the degree of overestimation varying across different model groups, suggesting there is no constant factor that can entirely correct for the effects. If a constant correction factor of 0.78, as suggested by [27], is applied, the corrected \(h_0\) values (shown as diamonds) do not reproduce the intrinsic values. In contrast, our correction method, represented by solid points, yields \(h_0\) values that align closely with the one-to-one relation, demonstrating the effectiveness of our approach in accounting for the projection effect caused by deviations from a perfect edge-on inclination. We emphasize that our method is statistical in nature and can be applied only to a sample, not to an individual source.

4.3 Redshift evolution of bias-corrected scale heights↩︎

Before examining the bias-corrected scale heights, we first probe the uncorrected measurements derived from the F444W, F277W, F150W, and F115W filters as a function of redshift (Figure 14). Only a weak trend is seen in the F444W filter, where the scale height increases slightly toward lower redshift and then remains roughly constant. At high redshift, the number of galaxies included in the F115W and F150W measurements decreases because in these filter the galaxies become too faint in surface brightness to yield reliable measurements and are therefore excluded (see Section 3).

To remove the bias introduced by variations in rest-frame wavelength, we have computed the scale height (\(h_{0,1\micron}\)), scale length (\(R_{d,1\micron}\)), and relative scale height (\(h_{0,1\micron}/R_{d,1\micron}\)) at a fixed rest-frame wavelength of 1 \(\micron\) (see Section 4.1). Their redshift dependence is shown in Figure 15. The relation between \(h_{0,1\micron}\) and redshift is more clearly defined than that derived from single-band measurements, as biases caused by variations in rest-frame wavelength, such as dust extinction and stellar population effects, have been removed. The \(R_{d,1\micron}\) exhibits a general increase toward lower redshifts, consistent with the well-established growth of galaxy sizes over cosmic time [66], [69], [70]. The observed \(h_{0,1\micron}/R_{d,1\micron}\) decreases with lower redshift, resulting in smaller values of \(f_{\rm corr}\) at low redshift. This indicates that the measured scale heights is more overestimated at lower redshifts than at higher redshifts.

To remove the bias introduced by the projection effect, we divide the sample into 7 redshift bins and compute the median observed value of \(h_{0,1\micron}/R_{d,1\micron}\) in each bin. These values are used as input to determine the \(f_{\rm corr}\) (see Figure 11 or Table [tab]), which ranges from 0.6 to 0.8. \(f_{\rm corr}\) is then applied to \(h_{0,1\micron}\) to obtain the bias-corrected scale height. The bias-corrected \(h_{0,1\micron}\) as a function of redshift is presented in Figure 15 (bottom-right) and are listed in Table [tab:results]. Therefore, after correcting for the effects of variation in rest-frame wavelength and the projection effect, the median and scatter of the scale heights are \(0.56\pm0.03\), \(0.61\pm0.04\), \(0.72\pm0.03\), \(0.78\pm0.03\), \(0.84\pm0.04\), \(0.82\pm0.03\), and \(0.67\pm0.06\) kpc at \(z=3.25\), \(2.75\), \(2.25\), \(1.75\), \(1.25\), \(0.75\), and \(0.25\), respectively. The scatters among individual measurements are 0.2–0.3 kpc. The overall trend shows an increase in scale height from \(z\sim3\) to \(z\sim1\), followed by a decline toward lower redshift. The bias-corrected ratios of \(h_{\rm vert, 1\micron}\) to \(R_{d,1\micron}\) and their inverse ratios are also presented in Table [tab:results].

Furthermore, to examine whether, during sample selection, the adopted cut in \(B/T\) used to remove potential elliptical galaxies, which may fall within the star-forming region of the color-color diagram due to measurement uncertainties or are intrinsically star-forming, affects our results, we plot the redshift evolution of the bias-corrected \(h_{0,1\micron}\) using a higher \(B/T\) threshold of 0.7, shown as the open circles in Figure 15. The result does not change significantly, suggesting that our conclusions are insensitive to adopting a higher \(B/T\) threshold.

4.4 Artificial radial increase of measured scale heights↩︎

The projection effect of a deviation from a perfectly edge-on inclination not only causes a systematic overestimation of the global vertical scale height but also produces an artificial radial increase of the measured height toward larger radii, as previously noted by [15]. In addition, the LSF cannot fully account for the light from the bright central region scattered into the outskirts of the disk, leaving a residual PSF effect that introduces further artificial flaring. Figure 16 illustrates these effects using 3D model disks with an intrinsic constant scale height. The measured radial profiles are obtained by fitting 1D vertical light profiles, convolved with a LSF if need, at each radius in the 2D projected image. For a perfectly edge-on disk without PSF convolution, the measured profile normalized by central scale height remains flat at 1. When a PSF is convolved, however, the measured scale height becomes artificially larger even though a LSF is applied when doing the 1D fitting (see Section 3). For a disk at \(i=70°\) without PSF convolution, the measured scale height artificially increases with radius due to projection effects; when the PSF is convolved, this artificial increase becomes slightly stronger owing to the combined impact of projection and residual PSF effects.

To examine whether the disk galaxies in COSMOS-Web exhibit genuine or artificial radial gradients, we measured the vertical scale height at each radius using the 1D method described in Section 3. We refrain from interpolating to rest-frame 1 \(\micron\) because the 1D profiles at individual radii, extracted from 2D flux maps, have relatively lower \(S/N\) than global measurements, which would propagate larger uncertainties through interpolation. Instead, we adopted results from the filter probing rest-frame wavelengths just above 1 , where the effects of dust extinction are small. The galaxies were grouped into two redshift bins: \(0 < z \leq 1\) and \(1 < z \leq 3.5\). For each galaxy, the radial profile of the scale height was normalized by its central value within one FWHM of the PSF, and the radius was normalized by the disk scale length \(R_d\). The resulting normalized profiles are shown in Figure 17, where large points and error bars indicate the median and scatter, respectively. In both redshift ranges, the measured scale height systematically increases with radius, giving the appearance of disk flaring.

However, as shown above, such a trend can arise purely from the projection effect and residual PSF effect. To understand whether this observed increase is physical or artificial, we estimated the median observed \(h_{0,1\micron}/R_{d,1\micron}\) for each redshift bin and derived the corresponding intrinsic \(h_{0,1\micron}/R_{d,1\micron}\) and median inclination. Using these parameters, we generated PSF-convolved 2D projected models from 3D model disks with constant intrinsic scale height. The 1D fitting method yields the predicted artificial radial trends shown as blue solid curves in Figure 17. The predicted trends closely reproduce the observed radial profiles at both low and high redshift. This agreement indicates that the apparent flaring is primarily a projection artifact, with a minor contribution from residual PSF effects, rather than a genuine structural feature. The implications of this result for the origin of disk scale heights are discussed in Section 5.

4.5 Comparison with scale height measurements in previous works↩︎

There are several studies on the redshift evolution of disk scale height based on HST or JWST images [25], [27], [40], [71]. We compare our results with their measurements in Figure 18. Our measurements are plotted as black points, with the black curve showing the best-fit relation and the gray curves indicating its uncertainties. The fitting method is described later in Section 5.

At the Solar Circle of the Milky Way, [9] reported exponential scale heights of 0.9 kpc and 0.3 kpc for the thick and thin disks, respectively, with uncertainties of about 20%. If these values are multiplied by the commonly used, although inappropriate, factor of 2 to convert to a sech\(^2\) scale height, one obtains 1.8 kpc and 0.6 kpc. Under this assumption, the Milky Way’s thick disk would appear substantially thicker than the extragalactic values with median less than \(\sim\,\)​1 kpc, whereas its thin disk would be comparable to the global scale height at \(z=0\), suggesting that the Milky Way might be an unusually thick system. However, as shown in Eq. (7 ), the appropriate conversion factor is 1.37 rather than 2. Using this conversion, the sech\(^2\) scale heights of the Milky Way’s thick and thin disks become 1.23 kpc and 0.41 kpc, respectively, as indicated by arrows in Figure 18. These values are a few hundred parsecs higher and lower, respectively, than the global scale heights at \(z\approx0\) in this work and in most of the previous studies. This is expected, because the global measurement represents a flux-weighted average of the thin and thick components.

Based on HST/ACS imaging, [25] measured the scale heights of 107 visually selected edge-on disks with an average stellar mass of \(10^{9.5\textendash10.5}\,M_\odot\) at \(0.5\leq z\leq 3.5\). They reported a mean scale height of \(0.63\) kpc with a scatter of 0.24 kpc in the F814W band, which probes rest-frame wavelengths of 0.20–0.54 and is therefore potentially affected by dust extinction. We regrouped their sample into three redshift bins and plot the results as light-gray diamonds in Figure 18. The data reveal an increase in scale height from \(\sim0.32\) kpc at \(z\approx3\) to \(\sim1.03\) kpc at \(z\approx1\).

[27] measured the vertical scale heights of 491 disk galaxies with stellar masses of \(M_* = 10^{9\textendash11}\,M_\odot\) (the majority below \(10^{10}\,M_\odot\)) over the redshift range \(0.4 \leq z \leq 2.5\) using HST/ACS F850LP, HST/WFC3 F125W, and F160W imaging. The images are tracing rest-frame optical wavelength range (0.46–0.66 ). The galaxies were selected with an axis ratio cut of \(q < 0.4\), and the scale heights were derived from 1D vertical profile fitting. A constant correction factor of 0.78 was applied to account for the projection effect, yielding a median scale height of \(0.74\) kpc with a scatter of 0.35 kpc. The median values of three redshift bins are shown as gray squares in Figure 18. They found little to no redshift evolution in the scale height. Nevertheless, the relatively shallow HST observations and their lower spatial resolution may introduce larger measurement uncertainties.

Using F115W imaging, which traces the rest-frame NIR at \(z\sim0\) but shifts to the rest-frame UV at \(z>3\), [39] measured \(\operatorname{sech}^2\) scale heights for 191 disk galaxies with a median stellar mass of \(\sim10^9\,M_\odot\) over \(0.2\le z\le5\), reporting a global median of \(0.38\) kpc with a scatter of 0.13 kpc. The red crosses in Figure 18 show their median values in each redshift bin (plotted only for \(z<3.5\)). Their measured global scale height exhibits only mild evolution at \(z>1.5\) but decreases toward lower redshift at \(z<1.5\), reaching a surprisingly small value of \(\sim0.18\) kpc at \(z=0.35\). Perhaps due to uncertainties caused by small sample statistics, this value is roughly half of the global scale heights reported by [40] for JWST galaxies of similar stellar mass at similar redshift. [39] further argued that scale-height measurements are consistent across NIRCam filters from F444W to F115W. However, as demonstrated in Section 4.1 and in [14], measurements at shorter rest-frame wavelengths are systematically biased toward larger values due to dust extinction across the mid-plane.

[40] studied 111 nearly edge-on disk galaxies at \(z \leq 3\) using F277W, F356W, and F444W images, tracing rest-frame wavelengths of \(1\)–\(2\,\micron\), where dust extinction effect is small. They identified 67 galaxies better fitted with a single-disk component and 44 galaxies requiring two disk components (thin and thick disks), with the earliest two-disk system found at \(z = 1.96\). After correcting for mass dependence, the median scale heights at \(M_*=10^{9.2}\,M_\odot\) for single disks, thin disks, and thick disks are shown in Figure 18. Consistent with our results, the scale height of single-disk galaxies increases toward lower redshift, with a median value of about \(0.66\pm0.07\,\)kpc, slightly smaller than our measurements due to their lower stellar masses. In contrast, the thin and thick components of two-disk galaxies show nearly constant, or perhaps slightly decreasing, scale heights with decreasing redshift, with median values of \(0.33\pm0.04\) and \(0.94\pm0.07\,\)kpc, respectively.

5 Implications for Galaxy Evolution↩︎

It is important to emphasize that the redshift evolution of disk scale height reported in this work refers to galaxies with a fixed stellar mass of approximately \(2.5\times10^{10}\,M_{\odot}\). In reality, galaxies increase their stellar mass over cosmic time; therefore, our analysis does not trace the exact progenitor–descendant relationship of individual systems across redshift. Nevertheless, the results offer valuable clues to the overall cosmic evolution of disk galaxies.

The formation mechanism of thin disks is relatively well understood. Thin disks are composed of stars born from a dynamically cold, rotationally supported gaseous disk. The newborn stars inherit the kinematics of the gas, exhibiting low ratios of velocity dispersion to rotation velocity, and consequently form geometrically thin stellar structures [10], [11]. However, open questions remain regarding when thin disks first emerged in the Universe and how they evolved over time.

In contrast, the dominant formation mechanism of thick disks remains debated despite decades of research. This uncertainty partly stems from observational limitations, particularly the difficulty in resolving the vertical structure of distant galaxies prior to the JWST era. Broadly speaking, proposed thick-disk formation mechanisms fall into two categories:

(1) Thick disks are born thick. In this picture, stars form with large scale heights of \(\sim\)​1 kpc either at birth or shortly thereafter. This may occur through rapid stellar scattering driven by violent gravitational instabilities in gas-rich, clumpy high-redshift disks [20], [72], [73], combined with the dissolution of massive clumps [25], [26]. It may also arise from star formation within turbulent high-redshift gaseous disks [74], [75], where giant molecular clouds (GMCs) have sizes comparable to clumps, typically \(\sim\)​300 pc [76]–[81]. Together with the effect that the GMCs may not in coplanar, these conditions may directly produce stars with high velocity dispersion, naturally giving rise to thick-disk populations and providing a plausible explanation for the thick disks found in hydrodynamical simulations [11], [12]. Furthermore, the sloshing effect that perturbation from non-axisymmetric structures and the stochastic star formation in the disk induce bulk motion in the disk, increasing velocity dispersion of the stars [82], [83].

(2) Thick disks were born thin but later thickened. The potential mechanisms include dynamical heating by minor mergers [21], [23], secular evolution by stellar structures [84], [85], and decreasing in surface mass density due to galaxy size growth.

5.1 Thick disks begin as intermediate-scale-height, dynamically hot, dense structures↩︎

The detection of relatively small disk scale heights at high redshifts was first reported by [25], who found that the median scale height of disks at \(z=3\) is \(\sim0.3\) kpc, which is although based on rest-frame \(\sim0.2\,\micron\) observations with HST. This result was recently reinforced by [40], who analyzed JWST NIR rest-frame images of 111 edge-on disks. After correcting their measurements to a median stellar mass of \(10^{9.2}\,M_\odot\), they obtained a median global disk scale height of 0.51 kpc at \(z=2.6\), compared to 0.77 kpc at \(z=0.26\). Consistently, using a JWST sample that is 15–30 times larger than in previous studies, we find a median scale height of \(0.56\) kpc at \(z=3.25\). Together, these findings suggest that high-redshift disks are neither thin nor thick, but instead occupy an intermediate scale-height regime between the typical thin disks (0.41 kpc for the Milky Way) and thick disks (1.23 kpc for the Milky Way), and experienced vertical growth over Gyr timescales. Consequently, the category (1) mechanisms may involve physical conditions in their galaxy models that do not fully match real observations, and may therefore require further improvement.

This, however, does not exclude the role of gravitational-instability-driven scattering [20] or the sloshing effect [83]. On the contrary, such processes likely play a crucial role in producing dynamically hot stellar populations within thick disks. In the same vein, the “born thick” disk stars formed from large-size GMCs or clumps, from bursty star formation seen in the simulations of [11] also contribute to form hot stars. Our results are not necessarily inconsistent with these simulations, because the bridge connecting simulations and observations, namely, realistic mock observations derived from simulations, is still lacking.

Although disks at \(z\approx3.5\) are geometrically intermediate-thickness rather than thin or thick, they are not scaled-up versions of local cold thin disks. Instead, they more closely resemble thick disks, as they are dynamically hot and exhibit high stellar mass surface densities, as discussed below. In Figure 19, we show that the stellar surface mass density, \(\Sigma_{*, \rm disk}=M_{*,\rm disk}/(2\pi R_{e,\rm{disk},1\micron}^2)\), increases toward higher redshift: the median value is \(\Sigma_{*, \rm disk}=10^{8.13}\,M_\odot\,\mathrm{kpc}^{-2}\) at \(z=0.25\), and \(\Sigma_{*, \rm disk}=10^{8.67}\,M_\odot\,\mathrm{kpc}^{-2}\) at \(z=3.25\). A best-fit relation to the median values in the seven redshift bins yields \[\label{SigmaMass} \Sigma_{*, \rm disk}(z) = 10^{8.07\pm0.02} \exp\!\left((0.34\pm0.02)\,z\right) ~ M_\odot\,\mathrm{kpc}^{-2}.\tag{8}\] Assuming an isothermal disk [63], the vertical stellar velocity dispersion is related to \(h_0\) and \(\Sigma_{*, \rm disk}\) through (also see derivation in the Appendix):

\[\label{Eq:h} h_0 = \frac{\sigma^2_{\rm *, vert}}{\pi\,G\,\Sigma_{*, \rm disk}},\tag{9}\]

or equivalently,

\[\sigma_{\rm *, vert} = \sqrt{\pi\,G\,\Sigma_{*, \rm disk}\,h_0}.\]

Given \(\Sigma_{*, \rm disk}=10^{8.67}\,M_\odot\,\mathrm{kpc}^{-2}\) and a bias-corrected \(h_{0,1\micron}=0.56\) kpc at \(z=3.25\), we obtain \(\sigma_{\rm *,vert}=59\,{\rm km\,s^{-1}}\). This value is consistent with \(\sigma_{\rm *,vert}=50\pm5\,{\rm km\,s^{-1}}\) for stars at the solar neighborhood in the Milky Way’s thick disk [5]. If including the contribution of the gaseous component, assumed to account for 40% of the total mass [86], the inferred velocity dispersion increases to \(\sim78\,{\rm km\,s^{-1}}\). Therefore, although the disks at \(z=3.25\) appear geometrically thin, their vertical velocity dispersions of 59–78 km s\(^{-1}\) indicate that they are dynamically hot. These hot stellar populations are likely produced by star scattering [20], sloshing effect [83], or formed inherently hot [11]. The strong gravitational potential of the dense mid-plane confines stars with high vertical velocity dispersions within a geometrically thin structure.

We may further gain insights from gas kinematics. Observations show that the velocity dispersion of the gaseous component, including the ionized, atomic, and molecular phases, increases systematically with redshift [87]. By compiling and analyzing measurements from the literature, [88] derived a redshift-dependent relation for the gas velocity dispersion: \(\sigma_{\rm g}/{\rm km\,s}^{-1} = 10.9 + 11.0\,z\). Applying this relation yields \(\sigma_{\rm g,vert}\sim50\,{\rm km\,s^{-1}}\) at \(z=3.25\), consistent with the inferred \(\sigma_{\rm *,vert}\) above. Such turbulent, gas-rich disks facilitate the formation of stars with high velocity dispersion [11]. Consequently, even though the stellar disks at \(z\approx3.5\) are geometrically intermediate-thickness rather than thick, with median \(h_{0,1\micron}=0.56\) kpc, they are dynamically hot. Furthermore, these \(z\approx3.5\) disks exhibit a single-component vertical structure and occupy a structural-parameter space comparable to that of thick disks identified at lower redshifts, while a distinct thin component only begins to emerge at \(z<2\) [40]. This may suggest that most galaxies initially form a thick disk, observed as a dense, dynamically hot, intermediate-thickness single component, followed by the subsequent development of a thin disk.

Because our sample is selected to have an approximately constant median stellar mass across redshift, it does not trace the evolutionary paths of individual galaxies. However, [40] similarly find relatively small scale heights (\(\sim\,0.51\) kpc) in high-redshift disks at a lower median mass of \(10^{9.2}\,M_\odot\). Combined, these findings point toward a picture in which present-day thick disks originate from dense, dynamically hot, intermediate-thickness disks formed in the turbulent, gas-rich conditions of early cosmic epochs.

5.2 Minor mergers and secular evolution do not drive present-day thick disks↩︎

Dynamical heating by minor mergers [21], [23] has long been proposed as a possible mechanism for thick disk formation, especially at early cosmic times when merger rates were significantly higher [89], [90]. This process, however, is predicted to induce strong vertical flaring of stellar disks [20], [91], i.e., an intrinsic increase of scale height with radius. Yet, observations suggest otherwise. [20] showed that high-redshift disks do not exhibit pronounced flaring, at least based on visual inspection of edge-on images. [26] did not find signature of minor mergers in their sample. Moreover, [15] pointed out that even a slight deviation from a perfectly edge-on orientation can produce artificial flaring, which may mimic the expected intrinsic trend.

As shown in Figure 17, the observed radial gradient of the measured scale height in our sample is consistent with the amount of artificial flaring expected from the projection effect and the residual PSF effect (see Section 4.4 for detailed calculation). This agreement implies that the apparent flaring in our sample is likely due to observational effects rather than intrinsic structural variations. Therefore, we find no evidence for intrinsic vertical flaring across our edge-on disk sample over the redshift range \(0<z<3.5\), effectively ruling out minor mergers as the dominant mechanism responsible for the formation of the thick disk component.

Nevertheless, mergers likely played an important role in specific cases such as the Milky Way, where they contributed to disk heating and the formation of in-situ halo stars with chemical abundances similar to those of thick disk stars [50], [92], [93].

Our results also do not support scenarios in which thick disks formed primarily through slow heating of a pre-existing thin disk by bar buckling [84] or spiral scattering [85]. These mechanisms are potentially important, as bars and spirals are common in nearby galaxies [94]–[97] and have been shown to play a key role in driving galaxy secular evolution [98]–[100]. If bar buckling or spiral scattering played a dominant role in disk thickening, one would expect disk thickness to increase toward lower redshifts at \(z < 1\), when bars and spiral arms are well developed and dynamically active. In contrast, we observe a clear decrease in disk thickness over this redshift range (Fig. 18). Furthermore, these secular mechanisms are also disfavored by the observed high-[\(\alpha\)/Fe] stars in the thick disk of the Milky Way. If they were responsible for building the thick disk over the last 8 Gyr of the Universe, the gas reservoir would inevitably produce a substantial number of low-[\(\alpha\)/Fe] stars that would be heated into the thick disk, which is contrary to observations. Consistently, hydrodynamical simulations also show no significant thick disk formation at \(z\lesssim1\), after the end of the bursty star-formation phase [11].

5.3 Scale heights from \(z=3\) to 1↩︎

From \(z=3.25\) to \(z=1.25\), the disk scale height increases from \(0.56\) to \(0.84\) kpc and \(R_d\) grows from \(1.52\) to \(2.56\) kpc, while \(\Sigma_{*, \rm disk}\) declines from \(10^{8.67}\) to \(10^{8.24}\,M_\odot\,\mathrm{kpc}^{-2}\). The errors of these values are approximately \(0.015\) kpc, \(0.05\) kpc, and \(10^{7.3}\,M_\odot\,\mathrm{kpc}^{-2}\), respectively. For comparison, \(\Sigma_{*, \rm disk}=10^{8.07}\,M_\odot\,\mathrm{kpc}^{-2}\) at \(z=0\) in our sample. As discussed below, these trends reflect a combination of effects, including the decreasing of \(\Sigma_{*, \rm disk}\) due to size growth, and the generation of dynamically hot stars formed through bursty star formation, stellar scattering, and the sloshing effect.

We first examine the impact of size growth. It is well established that, at fixed stellar mass, galaxy size increases toward lower redshift [66], [69], [70], likely as a result of dark-matter halo growth and radial stellar migration. Because our galaxies have a nearly constant stellar mass of \(10^{10.4}\,M_\odot\) across redshift, size growth naturally leads to a lower \(\Sigma_{*, \rm disk}\), which plays a key role in determining disk scale height. To test whether the decreasing of \(\Sigma_{*, \rm disk}\) due to size growth alone can account for the observed increase in \(h_0\), we assume that there are no exchanges of mass, energy, or angular momentum between the system and its environment to isolate the pure effect of size growth. The size growth is therefore assumed to be adiabatic. The scale height varies with \(\Sigma_{*, \rm disk}\) as follows (See Appendix for formula derivation): \[h_0 \propto \Sigma_{*, \rm disk}^{-0.2}.\] This relation predicts an increase in scale height by a factor of \(1.16\) from \(z=3.25\) to \(1.75\), smaller than the observed factor of \(1.4\). Hence, additional energy sources are required to raise the stellar velocity dispersion.

Both observations and simulations suggest that thick disks were established rapidly during a turbulent phase rather than through a gradual buildup at early cosmic times (\(z\gtrsim1\)) [6]–[8], [11], [12], [50]. During this phase, fueled by intense gas inflows along cosmological cold streams [101], high-redshift disks exhibited elevated star formation rates [102]–[105] and clumpy, irregular morphologies [28], [29], [106]–[108], driven by violent gravitational instabilities [74], [75]. In such physical conditions, dynamically hot stars that form the thick disk may form, via star scattering [20] or the sloshing effect that perturbation from non-axisymmetric structures and stochastic star formation in the disk induce bulk motion in the disk respect to the halo potential, increasing the velocity dispersion [82], [83].

If we consider an extreme case where bursty star formation and gravitational instabilities are strong enough to keep \(\sigma_{\rm *,vert}\) unchanged during size growth (isothermal thickening), the scale height follows (see Eq. (9 )) \[h_0 \propto \Sigma_{*, \rm disk}^{-1}.\] This predicts an increase in scale height by a factor of \(2.1\) from \(z=3.25\) to \(1.75\), larger than observed. Thus, the actual evolution lies between the adiabatic and isothermal limits. The observed rise in scale height at high redshift therefore reflects the combined effects of size growth, bursty star formation, and violent disk instabilities. These disks are single-component systems that serve as progenitors of present-day thick disks, followed by the subsequent development of a thin disk.

After \(z\approx1\)–2, the penetration of cold streams into massive halos becomes inefficient as halos develop shock-heated gas atmospheres, leading to a decline in gas accretion rates [101], [109], [110]. Consequently, galaxies become less clumpy, reducing stellar scattering and the rate of thickening. Simulations have demonstrated that at this epoch some stars with low vertical velocity dispersion begin to form [11] facilitating the formation of old thin disks.

Observations have also confirmed the presence of old thin disks [17], [25], [40]. [17] found old stellar populations in thin disks of S0 galaxies, and [25] inferred thin disks at \(z\sim2\) from vertical color gradients. Using JWST imaging, [40] decomposed disks at \(z>1\) and identified the first thin disk at \(z=1.96\), finding that the thin-disk fraction increases toward lower redshift. The detection of spiral arms and bars in high-redshift galaxies [28], [29], [35]–[38], [111]–[114] further supports the presence of dynamically cold, settled disks [115], [116]. [40] show that the thin disk forms in a downsizing way: thin disks appear to form earlier in more massive galaxies, and on average thin disks appear to dominate over the past 8–9 Gyr for the stellar mass bin \(2.5\times10^{10}\,M_\odot\), corresponding to the bend in Figure 18.

5.4 Scale heights from \(z=1\) to 0↩︎

At \(z\approx1\), the virialization of the inner circumgalactic medium begins to regulate cosmological cold streams, reducing gas inflow to the disk and marking the transition from the bursty phase to a steady star-forming phase that enables the formation of dynamically cold stellar disks [11], [101]. Although minor mergers may still occur during this later phase at \(z\lesssim1\), simulations suggests that they do not account for the majority of thick-disk stars [11], consistent with our result that no intrinsic flaring in the disks is observed. Instead, minor mergers mildly heat some thin-disk stars and contribute to the young-star tail of the thick-disk population [11]. Fueled by cold gas accretion into the mid-plane and leftover gas from previous phase, the thin disk gradually builds up through the steady star formation.

Toward \(z=0\), galaxy sizes at fixed stellar mass continue to increase [66], reducing \(\Sigma_{*, \rm disk}\) and the mid-plane gravitational potential, thereby tending to increase the scale height of the disk. This effect competes with thin-disk formation, which drives a reduction in the global scale height through two mechanisms. First, the deepening of the mid-plane gravitational potential by cold gas accretion contracts the thick disk vertically [71]. These effects become more pronounced with the progressive growth of the thin disk over cosmic time, driving the overall decline in global thickness. Second, the superposition of thin and thick components introduces a geometric effect that reduces the measured global scale height in single-disk sech\(^2\) fitting of the combined vertical profile. During this period, Type Ia supernovae enrich the interstellar medium with Fe, lowering the [\(\alpha\)/Fe] ratios characteristic of thin-disk stars [5].

The ratio of scale length to scale height (\(R_{d,1\micron}/h_{0,1\micron}\)), i.e., the inverse of the relative scale height, serves as a useful indicator of the galaxy’s dynamical state, roughly tracing the \(v/\sigma\) of the system [117]. The bias-corrected values are summarized in Table [tab:results] and shown in Figure 20. The ratio remains approximately constant at \(2.7\pm0.2\) for \(z>1.5\), but increases to \(4.0\pm0.4\) at \(z=0.25\), indicating that galaxies become increasingly rotationally supported toward lower redshift, which is consistent with the picture in which thin-disk growth dominates galaxy evolution in this epoch.

5.5 Simplified Model for Evolution and Its Fitting↩︎

Based on the discussion above, we construct a simplified model for the redshift evolution of scale height, aiming to gain further insight on the galaxy disk evolution (Figure 21). Since the mid-plane mass density is the key factor determining disk thickness, the model is formulated as a function of \(\Sigma_{*, \rm disk}\), which evolves with redshift and is connected back to the redshift according to the best-fit relation in Eq. (8 ). We begin by assuming the dense, dynamically hot, intermediate-thicness disk observed at \(z\approx3.5\) represent the progenitors, or “seed galaxies”, of present-day thick disks. An example of its vertical density profile is shown in the bottom-right inset panel of Figure 21. As previously demonstrated, the scale-height increasing, due to the decreasing in \(\Sigma_{*, \rm disk}\), violent gravitational instabilities and bursty star formation, between \(z=3.5\) and \(z=1\), is a intermediate process between the adiabatic process and isothermal regimes, we parameterize the relation between scale height and \(\Sigma_{*, \rm disk}\) as \[\label{medadia} h_0 = D\,\Sigma_{*, \rm disk}^{k},\tag{10}\] where \(D\) is a normalization constant and \(k\) lies between \(-0.2\) and \(-1\). Here we adopt kpc and \(M_{\odot}\,{\rm kpc}^{-2}\) as the unit of scale height and stellar mass surface density, respectively. For simplicity, the increase in the mass of gas component is absorbed into the factor \(D\), given that both the stellar surface density \(\Sigma_{*, \rm disk}\) and the gas fraction increase with redshift.

For \(z\lesssim 1\), the formation of thick disks are expected to significantly reduced. We therefore assume pure adiabatic thickening for the thick disk at \(z<1\), and use a step function \(\zeta\) to connect \(z>1\) and \(z<1\) formulae (before addition of the thin disk): \[\label{thick95ori} \begin{align} h_{\rm thick}^{0}(z) = D\biggl[ \zeta(z)\,\Sigma_{*,{\rm disk}}^{\,k+0.2}(z=1)\, \Sigma_{*,{\rm disk}}^{-0.2}(z) \\ + \bigl(1 - \zeta(z)\bigr)\,\Sigma_{*,{\rm disk}}^{\,k}(z) \biggr]. \end{align}\tag{11}\] where \(\Sigma_{*,{\rm disk}}^{k+0.2}(z=1)\) is to ensure the continuity at \(z=1\). Motivated by the detection of a thin disk as early as \(z\sim1.96\) and the increasing number fraction of thin disks toward lower redshift reported in [40], we adopt the following functional form for the mass fraction of thin disk: \[\label{fthin} f_{\rm thin} = \left\{ \begin{array}{ll} f_0\, (2-z)^{t}, & \text{for } z \leq 2 \\ 0, & \text{for } z > 2 \end{array} \right.\tag{12}\] where \(f_0\) and \(t\) characterize the amplitude and growth rate. We assume the mass-to-light ratio for the thick and thin disks is 1.2/1 [118]. [40] find that scale heights for thin disk are approximately unchanged with \(z\) with a median value of \(\sim\,0.33\,\)kpc. We therefore adopt a constant scale height for the thin-disk model: \[\label{thin95ori} h_{\rm thin}(z) = 0.43\,{\rm kpc}\tag{13}\] scaled upward by a factor of 1.3, difference between their and our global scale heights, to account for stellar mass differences between their sample and ours. This value is almost the same with the sech\(^2\) scale height of the Milky Way’s thin disk (0.41 kpc). An example vertical light profile showing a thick disk overlapped with a newly formed thin disk is displayed in the inset at \(z=1.5\), while a later stage with a more developed thin disk is illustrated at \(z=1\) in Figure 21.

Because our sample has constant stellar mass across redshift, the growth of the thin disk must be accompanied by a reduction of stellar mass in the thick component. To model this redistribution, we assume that the vertical velocity dispersion of thick-disk stars remains unchanged, with a modest increase caused by adiabatic contraction during thin-disk buildup. This is reasonable because the dominant heating mechanisms responsible for thick-disk formation operate primarily at \(z\gtrsim1\) and are largely ineffective at later times. The gradual buildup of the thin disk, with surface density \(f_{\rm thin}\Sigma_{*, \rm disk}(z)\), is assumed to proceed adiabatically, contracting the vertical structure and slightly increasing the velocity dispersion of stars of the pre-existing thick disk [25], [26]. We assume the approximation that the thin and thick components each maintain independent pressure equilibria, while their gravitational potentials are coupled [26], [119]. The connection between the scale heights of the thin and thick components in a thick-thin system can then be approximated as (see Appendix and [10]): \[\label{thick95response} h_{\rm thick} = \frac{\sigma_{\rm thick}^2}{\pi G \left( (1-f_{\rm thin})\Sigma_{*, \rm disk} + f_{\rm thin}\Sigma_{*, \rm disk} \dfrac{h_{\rm thick}}{h_{\rm thin}} \right)},\tag{14}\] where \(\sigma_{\rm thick}\) is linked to \(h_{\rm thick}^0\) and \(\Sigma_{*, \rm disk}\) through the adiabatic invariant of the thick disk (Eq. 19 ). The quantity \(h_{\rm thick}(z)\) is then solved numerically.

The corresponding vertical light profiles are given by \[I_{\rm thick}(h,z) = \frac{(1-f_{\rm thin})\,\Sigma_{*, \rm disk}(z)}{2\,h_{\rm thick}(z)}\,{\rm sech}^2\!\left( \frac{h}{h_{\rm thick}(z)}\right)\] for the thick disk and \[I_{\rm thin}(h, z) = \frac{f_{\rm thin}\,\Sigma_{*,\rm disk}(z)}{2\,h_{\rm thin}(z)} \,{\rm sech}^2\!\left(\frac{h}{h_{\rm thin}(z)}\right)\] for the thin disk. An example of the combined thin+thick profile at \(z=0\), showing a dominant thin disk embedded within a thicker component, is presented in Figure 21. To obtain the disk scale height averaging the thin and thick disk, we fit a single sech\(^2\) function to the combined profile, mimicking our observational procedure applied to edge-on disks in COSMOS-Web: \[\label{final95h} h_0(z) = \Gamma\!\left(\,I_{\rm thin}(h, z) + I_{\rm thick}(h, z)\,\right),\tag{15}\] where \(\Gamma\) denotes the fitting operation.

We fit Eq. (15 ) to the bias-corrected \(h_{0,1\micron}(z)\) data using four free parameters: \(D\), \(k\), \(f_0\), and \(t\). The best-fit function is shown as the black curve in Figure 18 and 21. To estimate uncertainties, we perform 1000 bootstrap resamplings with replacement and refit the model each time; the resulting best-fit curves are shown as gray lines in Figure 18 and 21. The resulting parameters yield \(D=10^{6.0}\pm10^{4.9}\), \(k=-0.73\pm0.02\), \(f_0=0.11\pm0.03\), and \(t=2.52\pm0.52\). The best-fit evolutionary curve for the thick-disk component is plotted as a red dashed line in Figure 21. The continued increasing of the thick-disk scale height below \(z=1\) is primarily due to the decline in \(\Sigma_{*, \rm disk}\) caused by size growth and mass reduction. The redshift corresponding to the peak of the best-fit curve is \(z_{\rm peak}=1.12\pm0.04\).

We assess the robustness of the fit by comparing the Bayesian Information Criterion (BIC) between Eq. (15 ) and a simple linear model. The BIC for Eq. 15 is lower by 20, indicating a strong statistical preference for our simplified model and confirming the presence of a peak in the evolution of disk scale height.

The best-fit thin-disk mass fraction, \(f_{\rm thin}\), is shown in Figure 22. It increases from 0 at \(z=2\) to \(0.10\pm0.03\) at \(z=1\) and \(0.60\pm0.10\) at \(z=0\). Our inferred \(f_{\rm thin}\) values are systematically higher than those reported by [40] for a sample with a median mass of \(10^{9.2}\,M_\odot\), likely because galaxies with higher stellar masses tend to host a larger fraction of thin disks ([18] and see Figure 10 of [40]). Another plausible explanation is that the decomposition of thin and thick disk components may be systematically incomplete relative to the idealized model expectations. Detecting two different components requires sufficient contrast between their vertical profiles, so galaxies with comparable thin and thick disk masses are more likely identified as having a two-component structure. In contrast, when either component dominates the flux, the secondary component becomes difficult to detect observationally and the galaxy is more likely classified as a single-component system. This may lead to the absence of very low thin-disk fractions in massive galaxies reported by [40] and [18]. Nevertheless, the inferred fraction at \(z=0\) agrees well with the measurements by [18], which gives \(f_{\rm thin}=0.66\pm0.09\) for nearby galaxies with \(M_*\sim10^{10}\)–\(10^{10.8}\,M_\odot\). This agreement suggests that, despite its simplicity, our model captures the essential physical processes governing the coevolution of thin and thick disks.

In summary, the observed redshift evolution of disk scale height from COSMOS-Web can be explained by the interplay between thick- and thin-disk formation and evolution. The rise in \(h_0\) from \(z=3.5\) to \(1\) results from size-driven decreasing of \(\Sigma_{*, \rm disk}\), violent gravitational scattering, and bursty star formation. The subsequent decline from \(z=1\) to \(0\) reflects geometric overlap with thin disks, continuing \(\Sigma_{*, \rm disk}\) decreasing, and contraction driven by thin-disk growth.

6 Summary and Conclusion↩︎

We conducted a comprehensive analysis of the global disk scale heights for a sample of 2631 clean, nearly edge-on disk galaxies at \(0<z<3.5\), observed by JWST/NIRCam in the COSMOS-Web survey with the F115W, F150W, F277W, and F444W filters. The sample has a median stellar mass of \(2.5\times10^{10}\,M_{\odot}\), and a disk-axis ratio cut of \(q_d < 0.4\) was applied to ensure edge-on orientations. We measured the vertical scale height, \(h_0\), by fitting a 2D model consisting of a bulge and an edge-on disk described by a single sech\(^2\) vertical profile. The disk size is characterized by the radial exponential scale length \(R_d\). The resulting \(h_0\) represents an average scale height that combines the contributions of both thin and thick disks, where such components coexist. For comparison, exponential vertical profiles were also tested. In addition, 1D vertical light profiles were fitted with a line spread function (LSF), constructed using a PSF, to derive scale heights as a function of radius. The measurements of scale heights and scale lengths, as well as the correction curves used to account for projection effects, can be accessed at Zenodo [68]. Our main findings are summarized as follows:

1. \(h_0\) derived from sech\(^2\) vertical profiles are 1.37 times as large as those obtained from exponential profiles, with negligible scatter.

2. \(h_0\) measured at shorter rest-frame wavelengths are systematically overestimated relative to those measured at longer wavelengths, primarily due to dust extinction across the mid-plane. This effect is consistent with [14]. We adopt a fixed rest-frame 1 for our measurements to mitigate the dust extinction and spatial variation in stellar population.

3. Deviations from a perfectly edge-on inclination systematically bias \(h_0\) toward larger values, with this projection effect becoming progressively more severe in relatively thinner disks (i.e., those with smaller \(h_0/R_d\)). We developed a new statistical correction method for this inclination bias, which is more accurate than the constant correction factor proposed by [27].

4. The observed radial variation of scale heights across our sample is consistent with artificial flaring primarily caused by projection effects, with a minor contribution from imperfect PSF treatment in the 1D fitting. This disfavors minor mergers as the dominant mechanism for the disk thickening and hence thick-disk formation, which are expected to produce genuine flaring.

5. After correcting for the biases caused by variation in observed rest-frame wavelength and by the projection effect, we find an unambiguous trend in scale height as a function of redshift: it increases from \(0.56\pm0.03\) kpc at \(z=3.25\) to \(0.84\pm0.04\) kpc at \(z=1.25\), and then decreases to \(0.67\pm0.06\) kpc at \(z=0.25\). The scatters among individual measurements are 0.2–0.3 kpc. The disks at \(z\sim3.5\) were geometrically intermediate-thickness, yet dynamically hot and highly dense, likely representing the progenitors of present-day thick disks.

6. The bias-corrected ratio of scale length \(R_d\) to \(h_0\) (\(R_d/h_0\)) remains approximately constant at \(2.7\pm0.2\) at \(z>1.5\), while increases to \(4.0\pm0.4\) at \(z=0.25\).

7. Using a simplified evolutionary model, we suggest that the early dense, dynamically hot, intermediate-thickness, single-disk galaxies at \(z\approx3.5\) thickened as a result of decreasing surface mass density (as galaxies grew in size), violent gravitational instabilities, and bursty star formation, leading to increasing scale height toward lower redshift at \(z>1\). Below \(z\sim1\), thin-disk growth became dominant, producing a vertically more compact structure with smaller scale heights toward \(z=0\).

8. The model fitting further indicates that the thin-disk mass fraction \(f_{\rm thin}\) increases from \(f_{\rm thin}=0\) at \(z\approx3.5\), to \(f_{\rm thin}=0.10\pm0.03\) at \(z=1\), and to \(f_{\rm thin}=0.60\pm0.10\) at \(z=0\). The best-fit value at \(z=0\) agrees well with the near-infrared observational constraints for nearby galaxies with comparable stellar mass reported by [18].

A detailed investigation of disk thickness as a function of galaxy properties, such as stellar mass, will be pursued in future work. Future studies focusing on large-sample thick-thin disk decompositions, as well as direct comparisons between realistic mock images from high-resolution hydrodynamical simulations and observations, will be essential for disentangling the physical processes that govern the evolution of disk thickness and the transition from high-redshift dense hot disks with intermediate-scale-height to present-day thin–thick disk systems.

7 Derivation of Key Equations for a Self-Gravitating Single-Disk Model↩︎

In the main text, we use \(h\) to denote the vertical distance instead of the commonly used \(z\), to avoid confusion with the symbol for the redshift. In the Appendix, however, we revert to \(z\) for vertical distance, as it is the convention most widely adopted in the literature.

For a stellar disk in vertical equilibrium, the balance between the vertical component of the gravitational force and the dynamical pressure (\(P=\rho\,\sigma^2\)) associated with the random stellar motions can be expressed as \[\frac{1}{\rho}\frac{d(\rho\,\sigma_z^2)}{dz} = -\frac{d\Phi}{dz},\] where \(\rho(z)\) is the mass density, \(\sigma_z\) is the vertical velocity dispersion, and \(\Phi(z)\) is the gravitational potential. This is the vertical component of the steady-state Jeans equation for a collisionless stellar system. If the system is isothermal in the sense that \(\sigma_z\) is constant with height, the equation reduces to \[\sigma_z^2\,\frac{d\rho}{dz} = -\rho\,\frac{d\Phi}{dz}.\] The gravitational potential obeys the Poisson equation \[\frac{d^2\Phi}{dz^2} = 4\pi G\rho.\] Combining these two equations and eliminating \(\Phi\) gives \[\label{A4} \frac{d^2}{dz^2}\ln\rho = -\frac{4\pi G}{\sigma_z^2}\rho,\tag{16}\]

whose analytic solution is

\[\rho(z) = \rho_0\,\mathrm{sech}^2\!\left(\frac{z}{z_0}\right),\] where \(\rho_0\) is the midplane density and \(z_0\) is the vertical scale height. Substituting this form into Eq. (16 ) yields \[\label{A6} z_0^2 = \frac{\sigma_z^2}{ 2\pi G\rho_0}.\tag{17}\]

The total surface density of the disk is \[\Sigma = \int_{-\infty}^{+\infty} \rho(z)\,dz = 2\rho_0 z_0.\] Substituting \(\rho_0 = \Sigma / (2z_0)\) into Eq. (17 ) gives \[\label{eq:sigma95relation} z_0 = \frac{\sigma_z^2}{\pi G \Sigma}.\tag{18}\]

If the disk evolves slowly in an adiabatic manner, the vertical pressure and density satisfy \[P = K\rho^{\gamma},\] where \(K\) is the adiabatic constant and \(\gamma=5/3\) for a collisionless stellar system. With \(P_0=\rho_0\sigma_z^2\) and \(\rho_0=\Sigma/(2z_0)\), we obtain \[\label{adiabatic95invar} \sigma_z^2 = K\left(\frac{\Sigma}{2z_0}\right)^{\gamma-1}.\tag{19}\] Combining this with Equation 18 yields \[\Sigma^{\,2-\gamma}z_0^{\,\gamma}=\mathrm{const.}\] For \(\gamma=5/3\), the scale height and velocity dispersion vary as \[z_0 \propto \Sigma^{-1/5}, \qquad \sigma_z \propto \Sigma^{2/5}.\] Thus, during a slow adiabatic expansion, the disk becomes slightly thicker and dynamically cooler as its surface density decreases.

8 Derivation of Key Equations for a Self-Gravitating Double-Disk Model↩︎

To derive an approximate formula that connects the scale heights of the two disk components, we assume that the thin and thick components each maintain independent pressure equilibria, while their gravitational potentials are coupled. We emphasize that, within this approximation, the two resulting sech\(^2\) profiles do not constitute a strictly self-consistent solution to the Poisson equation together with vertical hydrostatic balance. Even so, this approximation and the adoption of two sech\(^2\) profiles offer a practical and decent representation of the observed scale heights. Each disk follows a sech\(^2\) density profile:

\[\rho_i(z) = \rho_{0i} \operatorname{sech}^2\left(\frac{z}{z_{0i}}\right), \quad i = 1,2,\]

with corresponding surface densities

\[\Sigma_i = 2\rho_{0i} z_{0i}.\]

The total gravitational potential satisfies the Poisson equation:

\[\frac{d^2 \Phi}{dz^2} = 4\pi G \left[ \rho_1(z) + \rho_2(z) \right].\]

Each disk is assumed to be isothermal in the vertical direction, so that vertical hydrostatic equilibrium for each component gives

\[\sigma_i^2 \frac{d \ln \rho_i}{dz} = - \frac{d\Phi}{dz},\]

where \(\sigma_i\) is the vertical velocity dispersion of disk \(i\).

Taking the logarithmic derivative of the sech\(^2\) profile, we obtain

\[\frac{d \ln \rho_i}{dz} = - \frac{2}{z_{0i}} \tanh\left( \frac{z}{z_{0i}} \right),\]

so that hydrostatic equilibrium becomes

\[\frac{d\Phi}{dz} = \frac{2\sigma_i^2}{z_{0i}} \tanh\left( \frac{z}{z_{0i}} \right).\]

Differentiating once more with respect to \(z\) gives

\[\frac{d^2\Phi}{dz^2} = \frac{2\sigma_i^2}{z_{0i}^2} \operatorname{sech}^2\left(\frac{z}{z_{0i}}\right).\]

Evaluating at the mid-plane \(z=0\) where sech\(^2(0)=1\) and comparing with the Poisson equation yields

\[\frac{2\sigma_1^2}{z_{01}^2} = \frac{2\sigma_2^2}{z_{02}^2} = 4 \pi G (\rho_{01} + \rho_{02}).\]

Substituting \(\rho_{0i} = \Sigma_i / (2 z_{0i})\) gives the coupled algebraic equations for the scale heights:

\[\sigma_1^2 = \pi G \left( \Sigma_1 z_{01} + \Sigma_2 \frac{z_{01}^2}{z_{02}} \right), \quad \sigma_2^2 = \pi G \left( \Sigma_2 z_{02} + \Sigma_1 \frac{z_{02}^2}{z_{01}} \right).\]

Dividing the two equations provides a simple relation between the scale heights:

\[\frac{z_{01}}{z_{02}} = \frac{\sigma_1}{\sigma_2}.\]

The explicit expressions for \(z_{01}\) and \(z_{02}\) are then

\[\label{eq:z01z02} z_{01} = \frac{\sigma_1^2}{\pi G \left( \Sigma_1 + \Sigma_2 \dfrac{z_{01}}{z_{02}} \right)}, \quad z_{02} = \frac{\sigma_2^2}{\pi G \left( \Sigma_2 + \Sigma_1 \dfrac{z_{02}}{z_{01}} \right)}.\tag{20}\]

References↩︎