Bars and boxy/peanut bulges in thin and thick discs III. Boxy/peanut bulge formation and evolution in presence of thick discs

Boxy/peanut (b/p) bulges, the vertically extended inner parts of bars, are ubiquitous in barred galaxies in the local Universe, including our own Milky Way. At the same time, a majority of external galaxies and the Milky Way also possess a thick-disc. However, the dynamical effect of thick-discs in the b/p formation and evolution is not fully understood. Here, we investigate the effect of thick-discs in the formation and evolution of b/ps by using a suite of N-body models of (kinematically cold) thin and (kinematically hot) thick discs. Within the suite of models, we systematically vary the mass fraction of the thick disc, and the thin-to-thick disc scale length ratio. The b/ps form in almost all our models via a vertical buckling instability, even in the presence of a massive thick disc. The thin disc b/p is much stronger than the thick disc b/p. With increasing thick disc mass fraction, the final b/p structure gets progressively weaker in strength and larger in extent. Furthermore, the time-interval between the bar formation and the onset of buckling instability gets progressively shorter with increasing thick-disc mass fraction. The breaking and restoration of the vertical symmetry (during and after the b/p formation) show a spatial variation -- the inner bar region restores vertical symmetry rather quickly (after the buckling) while in the outer bar region, the vertical asymmetry persists long after the buckling happens. Our findings also predict that at higher redshifts, when discs are thought to be thicker, b/ps would have more 'boxy-shaped' appearance than more 'X-shaped' appearance. This remains to be tested from future observations at higher redshifts.


Introduction
A number of observational studies find that nearly half of all edge-on disc galaxies in the local Universe exhibit a prominent boxy or peanut-shaped structure (hereafter b/p structure; e.g.see Bureau & Freeman 1999;Lütticke et al. 2000;Erwin & Debattista 2017).A wide variety of observational and theoretical evidence indicates that many bars are vertically thickened in their inner regions, appearing as "boxy" or "peanutshaped" bulges when seen in edge-on configuration (e.g.Combes & Sanders 1981;Raha et al. 1991;Erwin & Debattista 2016).The presence of b/p structure has also been detected for galaxies in the face-on configurations, for example, in IC 5240 (Buta 1995), IC 290 (Buta & Crocker 1991), and several others (see examples in Quillen et al. 1997;McWilliam & Zoccali 2010;Laurikainen et al. 2011;Erwin & Debattista 2013).Several photometric and spectroscopic studies of the Milky Way bulge revealed that the Milky Way also has an inner b/p structure (e.g.see Nataf et al. 2010;Shen et al. 2010;Ness et al. 2012;Wegg & Gerhard 2013;Wegg et al. 2015).The occurrence of b/p bulges is observationally shown to depend strongly on the stellar mass of the galaxy, and a majority of barred galaxies above stellar mass log(M * /M ) ≥ 10.4 host b/p bulges (e.g.see Yoshino & Yamauchi 2015;Erwin & Debattista 2017;Marchuk et al. 2022).A similar (strong) stellar-mass dependence of the b/p bulge occurrence, at redshift z = 0, is shown to exist for the TNG50 suite of cosmological zoom-in simulations (Anderson et al. 2024).
Much of our current understanding of the b/p formation and its growth in barred galaxies is gleaned from numerical simulations.Studies using N-body simulations often find that, soon after the formation of a stellar bar it undergoes a vertical buckling instability, which subsequently gives rise to a prominent b/p bulge (e.g.see Combes et al. 1990;Raha et al. 1991;Merritt & Sellwood 1994;Debattista et al. 2004;Martinez-Valpuesta et al. 2006;Martinez-Valpuesta & Athanassoula 2008;Saha et al. 2013).Indeed, Erwin & Debattista (2016) detected two such local barred-spiral galaxies that are undergoing such a buckling phase.If a barred N-body model is evolved for a long enough time, it might go through a second and prolonged buckling phase, thereby producing a prominent X-shape feature (e.g.Martinez-Valpuesta et al. 2006;Martinez-Valpuesta & Athanassoula 2008).Furthermore, Saha et al. (2013) showed that a bar buckling instability is closely linked with the maximum meridional tilt of the stellar velocity ellipsoid (denoting the meridional shear stress of stars).Alternatively, a b/p bulge can be formed via the trapping of disc stars at vertical resonances during the secular growth of the bar (e.g.see Combes & Sanders 1981;Combes et al. 1990;Quillen 2002;Debattista et al. 2006;Quillen et al. 2014;Li et al. 2023) or by gradually increasing the fraction of bar orbits trapped into this resonance (e.g.see Sellwood & Gerhard 2020).The main difference between these two scenarios of b/p formation is that when the bar undergoes the buckling instability phase, the symmetry about the mid-plane is no longer preserved for a period of time (see discussion in Cuomo et al. 2023).
Regardless of the formation scenario, the b/p bulges are shown to have a significant effect on the evolution of disc galaxies by reducing the bar-driven gas inflow (e.g.see Fragkoudi et al. 2015Fragkoudi et al. , 2016;;Athanassoula 2016).The formation of b/p bulges can affect metallicity gradients in the inner galaxy (e.g.Di Matteo et al. 2014;Fragkoudi et al. 2017a) and can also lead to bursts in star formation history (e.g.Pérez et al. 2017).In addition, Saha et al. (2018) showed that for a 3D b/p structure (i.e.b/p seen in both face-on and edge-on configurations), it introduces a kinematic pinch in the velocity map along the bar minor axis.Furthermore, Vynatheya et al. (2021) demonstrated that for such a 3D b/p structure, the inner bar region rotates slower than the outer bar region.
On the other hand, a thick-disc component is now known to be ubiquitous in majority of external galaxies as well as in the Milky Way (e.g.see Tsikoudi 1979;Burstein 1979;Yoachim & Dalcanton 2006;Comerón et al. 2011aComerón et al. ,b, 2018)).The existence of this thick-disc component covers the whole range of the Hubble classification scheme, from early-type S0 galaxies to late-type galaxies (Pohlen et al. 2004;Yoachim & Dalcanton 2006;Comerón et al. 2016Comerón et al. , 2019;;Kasparova et al. 2016;Pinna et al. 2019a,b;Martig et al. 2021;Scott et al. 2021).The thick-disc component is vertically more extended and kinematically hotter as compared to the thindisc component.The dynamical role of a thick-disc on the formation and growth of non-axisymmetric structures has been been studied for bars (e.g., Klypin et al. 2009;Aumer & Binney 2017;Ghosh et al. 2023) and spirals (Ghosh & Jog 2018, 2022).Past studies demonstrated that the (cold) thin and (hot) thick discs are mapped differently in the bar and boxy/peanut bulge (e.g.Athanassoula et al. 2017;Fragkoudi et al. 2017b;Debattista et al. 2017;Buck et al. 2019).Since the presence of a thick disc can significantly affect the formation, evolution, and properties of bars (Ghosh et al. 2023), we need to explore how it will affect the b/ps, since b/ps are essentially the vertical extended part of the bar.
Stellar bars are known to be present in high-redshift (z ∼ 1) galaxies (e.g.see Sheth et al. 2008;Elmegreen et al. 2004;Jogee et al. 2004;Guo et al. 2023;Le Conte et al. 2023).Furthermore, a recent study by Kruk et al. (2019) showed the existence of b/p structure in high redshift (z ∼ 1) galaxies as well.At high redshift, discs are known to be thick, kinematically hot (and turbulent), and more gas rich.So, the question remains as to how efficiently the b/p structures can form in such thick discs at such high redshifts.Fragkoudi et al. (2017b) studied the effect of such a thick disc component on the b/p formation using a fiducial two-component thin+thick disc model where the thick disc constitutes 30% of the total stellar mass.The formation and properties of b/p bulges in multi-component discs (i.e. with a number of disc populations greater than two) was also studied in Di Matteo (2016), Debattista et al. (2017), Fragkoudi et al. (2018a,b), and Di Matteo et al. (2019).However, a systematic study of b/p formation in discs with different thin and thick discs, as well as composite thin and thick discs is still missing.We aim to pursue this here.
In this work, we systematically investigated the dynamical role of the thick-disc component in b/p formation and growth using a suite of N-body models with (kinematically hot) thick and (kinematically cold) thin discs.We varied the thick-disc mass fraction and considered different geometric configurations (varying ratio of thin-and thick-disc scale lengths) within the suite of N-body models.We quantified the strength and growth of the b/p in both the thin-and thick-disc stars and studied the vertical asymmetry associated with the vertical buckling instability.In addition, we investigated the kinematic phenomena (i.e.change in the velocity dispersion, meridional tilt angle) associated with the b/p formation and its subsequent growth.
The rest of the paper is organised as follows.Section 2 provides a brief description of the simulation setup and the initial equilibrium models.Section 3 quantifies the properties of the b/p structure, their temporal evolution, and the vertical asymmetry in different models and the associated temporal evolution.Section 4 provides the details of kinematic phenomena related to the b/p formation and its growth, while Sect. 5 provides the details of the relative contribution of the thin disc in supporting the X-shape structure.Section 6 contains the discussion, while Sect.7 summarises the main findings of this work.

Simulation setup, and N-body models
To motivate our study, we used a suite of N-body models consisting of a thin and a thick stellar disc, and the whole system is embedded in a live, dark-matter halo.One such model was already presented in Fragkoudi et al. (2017b).In addition, these models have been thoroughly studied in a recent work of Ghosh et al. (2023) in connection with a bar-formation scenario under varying thick-disc mass fractions.Here, we used the same suite of thin+thick models to investigate b/p formation and evolution with varying thick-disc mass fractions.
The details of the initial equilibrium models and how they are generated are already given in Fragkoudi et al. (2017b) and Ghosh et al. (2023).Here, for the sake of completeness, we briefly mention the equilibrium models.Each of the thin and thick discs is modelled with a Miyamoto-Nagai profile (Miyamoto & Nagai 1975), with R d , z d , and M d being the characteristic disc scale length, the scale height, and the total mass of the disc, respectively.The dark-matter halo is modelled with a Plummer sphere (Plummer 1911), with R H and M dm being the characteristic scale length and the total halo mass, respectively.The values of the key structural parameters for the thin and thick discs, dark-matter halo, and the total number of particles used to model each of these structural components are mentioned in Table 1.For this work, we analysed a total of 19 N-body models (including one pure thin-disc-only and three pure thick-disc-only models) of such thin+thick discs.
The initial conditions of the discs are obtained using the iterative method algorithm (see Rodionov et al. 2009).For further details, the reader is referred to Fragkoudi et al. (2017b) and Ghosh et al. (2023).The simulations are run using a TreeSPH code by Semelin & Combes (2002).A hierarchical tree method (Barnes & Hut 1986) with opening angle θ = 0.7 is used to calculate the gravitational force, which includes terms up to the quadrupole order in the multipole expansion.A Plummer potential was employed for softening the gravitational forces with a softening length = 150 pc.We evolved all the models for a total time of 9 Gyr.
Within the suite of thin+thick disc models, we considered three different scenarios for the scale lengths of the two disc (thin and thick) components.In rthickE models, both the scale lengths of thin and thick discs are kept the same (R d,thick = R d,thin ).In rthickS models, the scale length of the thick-disc component is shorter than that for the thin-disc one (R d,thick < R d,thin ), and in rthickG models, the scale length of the thick-disc component is Before we present the results, we mention that in our thin+thick models, we can identify and separate, by construction, which stars are members of the thin-disc component at initial time (t = 0) and which stars are members of the thick-disc component at t = 0, and we can track them as the system evolves self-consistently.Thus, throughout this paper, we refer to the b/p as seen exclusively in particles initially belonging to the thindisc population as the "thin-disc b/p" and that seen exclusively in particles initially belonging to the thick-disc population as the "thick-disc b/p".

Boxy/peanut formation and evolution for different mass fraction of thick-disc population
3.1.Quantifying the b/p properties Figure 1 shows the distribution of all stars (thin+thick) in the edge-on projection, calculated at the end of the simulation run (t = 9 Gyr), for all the thin+thick models considered here.In each case, the bar is placed in the side-on configuration (i.e.along the x-axis).A prominent b/p structure is seen in most of these thin+thick models.We further checked the same edge-on stellar density distribution, calculated for the thin-and thick-disc stars separately.Both of them show a prominent b/p structure in most of the thin+thick models.For the sake of brevity, we do not show it here (however, see Fig. 2 in Fragkoudi et al. 2017b).
3.1.1.Quantifying the b/p strength and its temporal evolution Here, we quantify the strength of the b/p structure and study its variation (if any) with thick-disc mass fraction.Following Martinez-Valpuesta & Athanassoula (2008) and Fragkoudi et al. (2017b), in a given radial bin of size ∆R (=0.5 kpc), we calculate the median of the absolute value of the distribution of particles in the vertical (z) direction, |z| (normalised by the initial value, z0 ) of a snapshot seen edge-on and with the bar placed side-on (along the x-axis).In Fig. 2, we show one example of the corresponding radial profiles of the |z|/z 0i (i = thin, thick, thin+thick) computed separately for the thin-and thick-disc particles, as well as for the thin+thick disc particles, as a function of time.As seen in Fig. 2, a prominent peak in the radial profiles (at later times) of |z|/z 0i denotes the formation of a b/p structure in the thin+thick model.Here, we mention that as the vertical scale height, and hence the vertical extent of thick disc is larger (by a factor of three) than that for the thin disc (see Table 1); the normalisation by z0 is necessary to unveil the intrinsic vertical growth due to the b/p formation.When only the absolute values of |z| are considered, the thick-disc stars always produce a larger value of |z| than the thin-disc stars.This happens due to the construction of equilibrium thin+thick models, and not due to the b/p formation.The normalised peak for the thin-disc is much larger than that for the thick-disc, in concordance with the previous results (Fragkoudi et al. 2017b).Furthermore, at later times, the peak in |z|/z 0i profiles shifts towards outer radial extent (more prominent for the thin disc b/p), indicating the growth of the b/p structure towards outer radial extent.These trends are also seen to hold for other thin+thick models that form a b/p structure during their evolutionary pathway.To quantify the temporal evolution of the b/p strength, we define the b/p strength at time t, S b/p (t) as the maximum of the peak value of the |z|/z 0i , i.e., In Martinez-Valpuesta et al. (2006) and Martinez-Valpuesta & Athanassoula (2008), a method based on the Fourier decomposition was formulated to quantify the strength of a b/p structure.
In Appendix A, we compare this method with Eq. ( 1) for the thin+thick model rthickE0.5.
Figure 3 shows the corresponding temporal evolution of the b/p strength, calculated separately for the thin and thick discs and for the composite thin+thick disc particles for all thin+thick models considered here.The thin-disc b/p remains much stronger than the thick-disc b/p structure, and this trend holds true for all thin+thick models with three different configurations (i.e.rthickS, rthickE, and rthickG) that develop a prominent b/p structure during their course of temporal evolution.
The temporal evolution of the b/p strength in three thickdisc-only models merits some discussion.For the rthickS1.0model, the values of S b/p show a monotonic increase with time, denoting the formation of a prominent b/p structure (also see the bottom row of Fig. 1).However, the final b/p strength for the rthickS1.0model remains lowest when compared to other rthickS models with different f thick values.For the rthickE1.01), for thin-disc (upper panels), thick-disc (middle panels), and total (thin+thick) disc stars (lower panels) for all thin+thick disc models with varying f thick values (see the colour bar).Left panels show the b/p strength evolution for the rthickS models, whereas middle panels and right panels show the b/p strength evolution for the rthickE and rthickG models, respectively.The thick disc fraction ( f thick ) varies from 0.1 to 0.9 (with a step-size of 0.2), as indicated in the colour bar.The blue solid lines in the middle and the bottom rows denote the three thick-disc-only models ( f thick = 1), whereas the red solid line in the top middle panel denotes the thin-disc-only model ( f thick = 0; for details, see text).model, the temporal evolution of S b/p shows a sudden jump around t = 6.5 Gyr and then remains constant.By the end of the simulation, this model forms a b/p structure which appears boxier than a peanut or X shape.For the rthickG1.0model, the temporal evolution of S b/p does not show much increment, and the model does not form a prominent b/p structure at the end of the simulation.
In addition, for a fixed value of f thick , we calculated the gradient of S b/p (t) with respect to time t (dS b/p (t)/dt) for three different geometric configurations considered here.One such example is shown in Fig. 4 for f thick = 0.5.A prominent (positive) peak in the dS b/p /dt profile denotes the onset of the b/p formation.As seen clearly, for a fixed f thick value, the peak in the dS b/p /dt profile occurs at an earlier epoch when compared with the other two geometric configurations.This confirms that the b/p forms at an earlier time in the rthickS0.5model compared to the rthickE0.5 and rthickG0.5models.These trends are in tandem with the fact that the bars in rthickS models form at earlier times and grow faster compared to the other two disc configurations (for details, see Ghosh et al. 2023).
Lastly, in Fig. 5 (top panel), we show the final b/p strength (i.e.calculated at t = 9 Gyr using the thin+thick disc particles) for all thin+thick models considered here.We point out that, in some cases, the maximum of the |z|/z 0i profile for the thick-disc is not always easy to locate; sometimes they display a plateau rather than a clear maximum (see Fig. 2).This, in turn, might have an impact in the estimation of the b/p strength.In order to derive an estimate of the uncertainty on the b/p strength measurement, we constructed a total of 5000 realisations by resampling the entire population using a bootstrapping technique (Press et al. 1986), and for each realisation we computed the radial profiles of |z|/z 0i as well as its peak value (denoting the b/p strength).The resulting error estimates are shown in Fig. 5.The final b/p strength shows a wide variation with the f thick values as well as with the thin-thick-disc configuration.To illustrate, for the rthickS models, the final b/p strength decreases monotonically as the f thick value increases.For the rthickE models, the final b/p strength increases from f thick = 0.1−0.3, and then decreases monotonically as the f thick value increases.A [kpc] similar trend is also seen for the rthickG models.Nevertheless, the strength of the b/p shows an overall decreasing trend with increasing f thick values, and this remains true for all three geometric configurations considered here.increases significantly (by a factor of ∼2) over the entire evolutionary phase.In addition, towards the end of the simulation run, the thick disc b/p is larger (by ∼10−15%) than the thin disc b/p.We found a similar trend in temporal variation of the b/p extent for other thin+thick models, and therefore they are not shown here.
We further checked how the extent of the b/p structure, by the end of the simulation run, varied with the thin-to-thick disc mass fraction ( f thick ).In Fig. 5 (bottom panel), we show the corresponding extents of the b/p structure computed at t = 9 Gyr using the thin+thick stellar particles for all thin+thick models considered here.The extent of the b/p structure increases steadily as the f thick value increases, and this trend holds true for all three different configurations considered here.Furthermore, at a fixed f thick value, the rthickG models show a higher value for the R b/p when compared to other two configurations, thereby denoting that rthickG models form a larger b/p structure, by the end of the simulation run, as compared to rthickE and rthickS models.Lastly, in Appendix C, we show how the extent of the thin disc b/p and thick disc b/p, at the end of the simulation (t = 9 Gyr), vary across different f thick values and different disc configurations.

Vertical asymmetry and buckling instability
To quantify the vertical asymmetry and the time at which it occurs, we first calculated the amplitude of the first coefficient in the Fourier decomposition (A 1z ), which provides a measure of the asymmetry (e.g.see Martinez-Valpuesta et al. 2006;Martinez-Valpuesta & Athanassoula 2008;Saha et al. 2013).The first Fourier coefficient (A 1z ) is defined as (Martinez-Valpuesta & Athanassoula 2008) where m i is the mass of the ith particle, and ϕ i is the angle of ith particle measured in the (x, z)-plane with the bar placed along the x-axis (side-on configuration).M tot is denotes the total mass of the particles considered in this summation.Following Martinez-Valpuesta & Athanassoula (2008), to make this coefficient more sensitive to a buckling, we only included the stars (see Eq. ( 2)) that are momentarily within the extent of the b/p (R b/p ) in the summation.The corresponding temporal evolution of the buckling amplitude (A 1z ), calculated separately for thin, thick, and thin+thick particles, for the model rthickE0.5 is shown in Fig. 7 (left panel).A prominent peak in the A 1z profile denotes the vertical buckling event.We further checked that for all the thin+thick models, a peak in the A 1z is associated with the dip/decrease in the bar strength.This is expected since it is well known that the bar strength decreases as it goes through the buckling phase.The A 1z amplitude is larger for the thin-disc stars when compared with that of the thick-disc stars.This is consistent is with the scenario that the thin disc b/p is stronger than the thick disc b/p.Another way of quantifying the buckling instability is by measuring the buckling amplitude, A buck , which is defined as (for details, see Sellwood & Athanassoula 1986;Debattista et al. 2006Debattista et al. , 2020) where m j , z j , and φ j denote the mass, vertical position, and azimuthal angle of the jth particle, respectively, and the summation runs over all star particles (thin, thick, thin+thickwhichever is applicable) within the b/p extent.The quantity A buck denotes the m = 2 vertical bending amplitude (for further details, see Debattista et al. 2006Debattista et al. , 2020)).The corresponding temporal evolution of the buckling amplitude (A buck ), calculated separately for thin, thick, and thin+thick particles for the model rthickE0.5 is shown in Fig. 7 (right panel).A prominent peak in the A buck profile denotes the onset of the vertical buckling instability.Furthermore, the peak of value of A buck is higher for the thin-disc than that for the thick-disc.This is again consistent with the thin disc b/p being stronger than the thick disc b/p.This trend holds true for all the thin+thick models considered here.Next, we define τ buck as the epoch for the onset of the buckling event when the peak in A buck occurs.As seen from Fig. 7, the epoch at which the peak in A 1z occurs coincides with τ buck .This is not surprising as both the quantities denote the same physical phenomenon of vertical buckling instability.In Sect.3.3, we further investigate the variation of τ buck with f thick and its connection with bar-formation epoch.While Fig. 7 clearly demonstrates the temporal evolution of the vertical asymmetry associated with the b/p structure formation, it should be borne in mind that A 1z (quantifying vertical asymmetry) or A buck (quantifying the m = 2 vertical bending amplitude) only informs us about the buckling instability in an average sense, and hence it lacks any information about the 2D distribution of the vertical asymmetry.To investigate that, we computed the 2D distribution of the mid-plane asymmetry.Following Cuomo et al. (2023), we define where Σ(x, z) denotes the projected surface number density of the particles at each position of the image of the edge-on view of the model.The resulting distribution of A Σ (x, z), computed separately for the thin-disc, thick-disc, and thin+thick discs at six different times (before and after the buckling happens) are shown in Fig. 8 for the model rthickE0.5.At the initial rapid  2)) denoting vertical asymmetry in bar region, calculated using thin-disc (blue lines), thick-disc (red lines), and thin+thick disc (black lines) for the model rthickE0.5.Right panel: temporal evolution of buckling amplitude, A buck (Eq.( 3)), calculated using thin-disc (blue lines), thick-disc (red lines), and thin+thick disc (black lines) for the model rthickE0.5.The vertical magenta dotted line denotes the onset of buckling instability (τ buck ), calculated from the peak of A buck profile (for details, see text).Furthermore, for reference, we indicate the onset of bar formation (τ bar , vertical maroon dotted line), calculated from the amplitude of the m = 2 Fourier moment.
bar growth phase (t ∼ 1 Gyr), the A Σ (x, z) values remain close to zero, indicating no breaking of vertical symmetry in that evolutionary phase of the model.Around t ∼ 2.7 Gyr, the model undergoes a strong buckling event (see the peak in Fig. 7).As a result, the distribution A Σ (x, z) shows large positive and negative values at t = 2.75 Gyr, thereby demonstrating that the vertical symmetry is broken around the mid-plane.At a later time (t = 5 Gyr), the vertical symmetry is restored in the inner region (as indicated by A Σ (x, z) ∼ 0).However, in the outer region (close to the ansae or handle of the bar), A Σ (x, z) still displays non-zero values, thereby indicating that the vertical asymmetry still persists in the outer region.Around t ∼ 6.85 Gyr, the model undergoes a second buckling event (see the second peak in A 1z , albeit with smaller values in Fig. 7).As a result, at a later time (t = 6.95Gyr), the model still shows non-zero values for A Σ (x, z) in the outer region.By the end of the simulation run (t = 9 Gyr), the values of A Σ (x, z) become close to zero throughout the entire region, thereby demonstrating that the vertical symmetry is finally restored.As Fig. 8 clearly reveals that the thin-disc stars show a larger degree of vertical asymmetry (or equivalently larger values of A Σ (x, z)) when compared with the thick-disc stars, we further checked the distribution of A Σ (x, z) in the (x, z)-plane at different times for other thin+thick models that host a prominent b/p structure.We found an overall similar trend of spatio-temporal evolution of the A Σ (x, z) as seen for the model rthickE0.5.

Correlation between bar and b/p properties
Past theoretical studies of the b/p formation and its subsequent growth have revealed a strong correlation between the (maximum) bar strength and the resulting (maximum) b/p strength (e.g.see Martinez-Valpuesta & Athanassoula 2008).We tested this correlation for the suite of thin+thick models considered here.The maximum bar strengths for the models are obtained from Ghosh et al. (2023) where we studied, in detail, the bar properties for these models.The maximum b/p strengths for the models are obtained from Eq. ( 1).We mention that all stellar particles (thin+thick) are used to calculate the maximum bar and b/p strengths for all models.The resulting distribution of the thin+thick models in maximum bar-maximum b/p strengths are shown in Fig. 9.As is seen clearly from Fig. 9, a stronger bar in a model produces a stronger b/p structure.This correlation holds true for all three geometric configurations considered here.Therefore, we also find a strong correlation between the (maximum) bar strength and the resulting (maximum) b/p strength in our thin+thick models, in agreement with past findings.In addition, we investigated the correlation (if any) between the lengths of the bar and the b/p in our thin+thick models.In Fig. 10 (top panel), we show the temporal evolution of the ratio of the b/p length (R b/p ) and the bar length (R bar ) for the model rthickE0.5.The bar length, R bar , is defined as the radial location where the amplitude of the m = 2 Fourier moment (A 2 /A 0 ) drops to 70% of its peak value (for a detailed discussion, see recent work by Ghosh & Di Matteo 2024).As is clearly seen, the ratio increases shortly after the formation of b/p, and the ratio almost saturates by the end of the simulation run (9 Gyr).Furthermore, we calculated the b/p length (R b/p ) and the bar length (R bar ), at the end of the simulation (9 Gyr) for all thin+thick models considered here.This is shown in Fig. 10 (bottom panel).For the rthickE and rthickG models, the ratio increases progressively with increasing f thick values.However, for the rthickS model, the ratio increases monotonically until f thick = 0.7, and then it starts to decrease.Lastly, we investigated the time delay between the bar formation and the onset of buckling instability for all thin+thin models considered here, and we studied if and how it varies with thick-disc mass fraction ( f thick ).In Appendix D, we show how the bar-formation epoch, τ bar , varies with different f thick values and with different disc configurations.Similarly, in Sect.3.2, we define the epoch of buckling instability when the peak in A buck occurs.The resulting variation of the time delay, τ buck − τ bar with f thick is shown in Fig. 11.For a fixed geometric configuration (rthickE, rthickS, or rthickG), the time interval between the bar formation and the onset of buckling instability becomes progressively shorter with increasing f thick values.This happens due to the fact that with increasing f thick , the bar forms progressively at a later stage (see Appendix D).In addition, for a fixed f thick value, the rthickS models almost always show shorter time delay (τ buck − τ bar ) when compared to other two geometric configurations considered here (see Fig. 11).

Kinematic signatures of buckling and its connection with b/p formation
Understanding the temporal evolution of the diagonal and offdiagonal components of the stellar velocity tensor was shown to be instrumental for investigating the formation and growth of the b/p structure (see Saha et al. 2013, and references therein).
Here, we systematically studied some key diagonal and offdiagonal components of the stellar velocity dispersion tensor and their associated temporal evolution for all the thin+thick models considered.
At a given radius R, the stellar velocity dispersion tensor is defined as (Binney & Tremaine 2008) where quantities within the brackets denote the average quantities, and the averaging is done for a group of stars.Here, i, j = R, φ, z.The corresponding stress tensor of the stellar fluid is defined as where ρ(R) denotes the local volume density of stars at a radial location R. τ n and τ s denote the normal stress (acting along the normal to a small differential imaginary surface dS ) and the shear stress (acting along the perpendicular to the normal to dS ), respectively (for further details, see Binney & Tremaine 2008;Saha et al. 2013).The components of τ n are determined by the diagonal elements of the velocity dispersion tensor, while the shear stress is determined by the off-diagonal elements of the velocity dispersion tensor (for details, see Binney & Tremaine 2008).Furthermore, the diagonal elements of the velocity dispersion tensor determines the axial ratios of the stellar velocity ellipsoid with respect to the galactocentric axes (ê R , êφ , êz ), whereas the orientations of the velocity ellipsoid are determined by the off-diagonal elements of the velocity dispersion tensor (for details, see Binney & Tremaine 2008;Saha et al. 2013).One such quantity of interest is the meridional tilt angle, which is defined as The tilt angle, Θ tilt , denotes the orientation or the deformation of the stellar velocity ellipsoid in the meridional plane (R−z-plane).In the past, it was shown for an N-body model that when the bar grows, it causes much radial heating (or equivalently, increasing the radial velocity dispersion, σ RR ) without causing a similar degree of heating the vertical direction (or equivalently, no appreciable increase in the vertical velocity dispersion, σ zz ).Consequently, the model goes through a vertical buckling instability causing the thickening of the inner part, which in turn also increases σ zz (e.g.Debattista et al. 2004;Martinez-Valpuesta et al. 2006  Fragkoudi et al. 2017b;Di Matteo et al. 2019).Therefore, it is of great interest to investigate the temporal evolution of the verticalto-radial velocity dispersion (σ zz /σ RR ) in order to fully grasp the formation and growth of b/p structure in our thin+thick models.In addition, using N-body simulations, Saha et al. (2013) demonstrated that during the onset of the buckling phase, the model shows a characteristic increase in the meridional tilt angle, Θ tilt , which in turn could be used as an excellent diagnostic to identify an ongoing buckling phase in real observed galaxies.In this work, we studied the temporal evolution of these dynamical quantities in detail for all the thin+thick models considered.
In Appendix E, we show the radial profiles of vertical-toradial velocity dispersion as a function of time for the model rthickE0.5.In order to quantify the temporal evolution of σ zz /σ RR , we computed them using Eq. ( 5 of the vertical-to-radial velocity dispersion (σ zz /σ RR ), calculated separately for thin, thick, and thin+thick particles, is shown in Fig. 12 (left panel) for the model rthickE0.5.As seen from Fig. 12, the temporal evolution of σ zz /σ RR displays a characteristic U shape (of different amplitudes) during the course of the evolution, arising from the radial heating of the bar (increase in σ RR ) and the subsequent vertical thickening (increase in σ zz ) due to the buckling instability.In addition, the temporal profiles of σ zz /σ RR for the thin disc show a larger and more prominent U-shaped feature when compared with that for the thickdisc stars.This is consistent with the notion that thin-disc b/p are, in general, stronger than the thick-disc b/p.The epoch corresponding to the maximum increase in the quantity σ zz /σ RR coincides with the peak in the A buck , denoting strong vertical buckling instability (see the location of vertical magenta line in Fig. 12).This further demonstrates that the b/p structures in our thin+thick models are indeed formed through vertical buckling instability.In Appendix E, we show the temporal evolution of σ zz /σ RR for all thin+thick models considered here.We checked the trends mentioned above hold true for all the models that formed a b/p structure during the evolutionary trajectory.Lastly, we investigated the temporal evolution of the meridional tilt angle, Θ tilt for the model rthickE0.5. Figure 12 (right panel) shows the corresponding temporal evolution of the tilt angle, Θ tilt , calculated separately for thin, thick, and thin+thick particles using Eq. ( 7) for the model rthickE0.5.The temporal evolution of Θ tilt shows a characteristic increase during the course of evolution.The epoch of the maximum value of the tilt angle coincides with the epoch of strong buckling instability (see the location of vertical magenta line in Fig. 12).This is in agreement with the findings of Saha et al. (2013) and is consistent with a "b/p formed through buckling" scenario.Furthermore, the temporal profiles of Θ tilt for the thin disc shows a larger and more prominent peak when compared with that for the thickdisc stars.This is expected as the thin-disc b/p is stronger than the thick-disc b/p.We checked that the trends, mentioned above, hold true for all the models that formed a b/p structure during the evolutionary trajectory.For the sake of brevity, they are not shown here.

X-shape of the b/p and relative contribution of thin disc
A visual inspection of Fig. 1 already revealed that for a fixed geometric configuration, the appearance of the b/p structure changes from more X-shaped to more boxy as the thick-disc mass fraction steadily increases.This trend is more prominent for the rthickE and rthickG models (see middle and right panels of Fig. 1).Here, we investigate this in further detail.In addition, we also investigate how the thin-and thick-disc stars contribute to the formation of the b/p structure.
To carry out the detailed analysis, we first chose two models, namely, rthickE0.1 and rthickE0.9.In the rthickE0.1 model, the thin-disc stars dominate the underlying mass distribution, whereas in the rthickE0.9model, the thick-disc stars dominate the mass distribution, thereby providing an ideal scenario for the aforementioned investigation.In Fig. 13 (top left panels), we show the density contours of the edge-on stellar (thin+thick) distribution (with bar placed along the x-axis) in the central region encompassing the b/p structure.This clearly brings out the stark differences in the morphology of density contours.In the rthickE0.1 model, the contours have a more prominent X-shaped appearance, whereas in the rthickE0.9model, the contours have a more prominent box-shaped appearance.Figure 13 (bottom left panels) shows the density profiles along the bar major axis, calculated at different heights (from the mid-plane) for these two thin+thick models.At larger heights, a bimodality in the density profiles along the bar major axis reconfirms the strong X-shaped feature for the rthickE0.1 model.For the rthickE0.9model, no such bimodality in the density profiles along the bar major axis is seen, thereby confirming that the b/p structure is more box-shaped in the rthickE0.9model.Furthermore, we calculated the vertical stellar density distribution at a radial location around the peak of the b/p structure (see vertical blue lines in top left panels of Fig. 13) for these two models.A careful inspection reveals that the vertical stellar density distribution for the rthickE0.1 model is more centrally peaked (with well-defined tails), whereas vertical stellar density distribution for the rthickE0.9model is broader, especially at larger heights (see right panel of Fig. 13).
To further investigate how the thin-and thick-disc stars contribute to the formation of the b/p structure, we calculated the thin-disc mass fraction, f thin (=1 − f thick ) at the end of the simulation run (t = 9 Gyr).Figure 14 shows the corresponding distribution of the thin-disc mass fraction in the edge-on projection (x−z-plane) for all thin+thick disc models.In each case, the bar is placed in side-on configuration (along the x-axis).As is seen clearly from Fig. 14, the thin-disc stars dominate in central regions close to the mid-plane (z = 0) and are responsible for giving rise to a strong X-shape of the b/p structure, in agreement with the findings of past studies (see e.g.Di Matteo 2016; Athanassoula et al. 2017;Debattista et al. 2017;Fragkoudi et al. 2017bFragkoudi et al. , 2020)).As one moves farther away from the mid-plane, the thick-disc stars start to dominate progressively.Furthermore, the appearance of the b/p structure changes from more X-shaped to more box-shaped as the thick-disc mass fraction steadily increases.These trends remain generic for all three geometric configurations (different thin-to-thick disc scale length ratios) considered here.

Discussion
In what follows, we discuss some of the implications and limitations of this work.First, our findings demonstrate that the b/p  7), for thin (in blue), thick (in red), and total (thin+thick) disc (in black) particles, for model rthickE0.5.The vertical maroon dotted line denotes the onset of bar formation (τ bar ), while the vertical magenta dotted line denotes the onset of buckling instability (τ buck ; for details, see text).
Fig. 13.Variation of the b/p morphology with varying thick-disc mass faction ( f thick ).Top left panels: density contours of edge-on stellar (thin+thick) distribution (with bar placed along the x-axis) in central region at t = 9 Gyr for models rthickE0.1 and rthickE0.9.For the rthickE0.1 model, the contours display more prominent X-shaped feature, whereas for the rthickE0.9model, the contours display more prominent box-shaped feature.Bottom left panels: density profiles (normalised by the peak density value, and in log-scale) along bar major axis, calculated at different heights (from |z| = 0 to 6 kpc, with a step-size of 1 kpc) from mid-plane, for models rthickE0.1 and rthickE0.9.The density profiles have been artificially shifted along the y-axis to show the trends as height changes, and they do not overlap.Right panel: vertical stellar density distribution (normalised by the peak density value, and in log-scale), at a radial location around peak of the b/p structure (marked by blue vertical lines in top left panels) for rthickE0.1 and rthickE0.9models.Fig. 14.Fraction of thin-disc stars, f thin (=1− f thick ), in the edge-on projection (with the bar placed along the x-axis) compared to the total (thin+thick) disc, at the end of the simulation run (t = 9 Gyr) for all thin+thick disc models with varying f thick values.Left panels show the distribution for the rthickS models, whereas middle panels and right panels show the distribution for the rthickE and rthickG models, respectively.The values of ( f thick ) vary from 0.1 to 0.9 (top to bottom), as indicated in the left-most panel of each row.For each model, the fraction of thin-disc stars decreases with height from the mid-plane (z = 0).In addition, the appearance of the b/p structure changes from more X-shaped to more boxy-shaped as the thick-disc mass fraction steadily increases.
structure can form even in the presence of a massive thick-disc component.This provides a natural explanation to the presence of b/p in high-redshift (z = 1) disc galaxies in the hypothesis that these high-z discs have a significant fraction of their mass in a thick disc (see e.g.Hamilton-Campos et al. 2023).A recent work by Kruk et al. (2019) estimated that at z ∼ 1, about 10% of barred galaxies would harbour a b/p.Therefore, our results are in agreement with the recent observational trends.In addition, bars forming in the presence of a massive thick disc (as shown in a recent work of Ghosh et al. 2023) and the present work showing that b/ps also form in the presence of a massive thick disc suggest that bar and b/p bulge formation may have appeared at earlier redshifts than what has been considered so far, as well as in galaxies dominated by a thick-disc component.The findings that the b/p morphology and length depend on the thick-disc fraction (the higher the thick-to-thin disc mass ratio, the more boxy the corresponding b/p and the smaller its extent) may be taken as a prediction that can be tested in current and future observations (the JWST, for example).
Secondly, the occurrence of b/p in disc galaxies is observationally found to be strongly dependent on the stellar mass of the galaxy in the local Universe (Yoshino & Yamauchi 2015;Erwin & Debattista 2017;Marchuk et al. 2022) as well as at higher redshifts (Kruk et al. 2019).This implies that a galaxy's stellar mass is likely to play an important role in forming b/ps via vertical instabilities.In our suite of thin+thick models, although we systematically varied the thick-disc mass fraction, the total stellar mass remained fixed (∼1×10 11 M ).Investigating the role of the stellar mass on b/p formation via vertical buckling instability would be quite interesting; however, it is beyond the scope of the present work.
Furthermore, in our thin+thick models, the stars are separated into two well-defined and distinct populations, namely, thin-and thick-disc stars.While this might be well-suited for external galaxies, this scheme is a simplification for the Milky Way.Bovy et al. (2012) showed that the disc properties vary continuously with the scale height in the Milky Way.Nevertheless, our adapted scheme of discretising stars with a varying fraction A196, page 13 of 19 of thick-disc stars provides valuable insight into the trends, as it has been shown for the MW that a two-component disc can already capture the main trends found in more complex, multicomponent discs (e.g.see Di Matteo 2016; Debattista et al. 2017;Fragkoudi et al. 2017bFragkoudi et al. , 2018aFragkoudi et al. ,b, 2020)).
Finally, if these b/p structures also formed at high redshift, two additional questions arise.These concern the role of interstellar gas, which is particularly critical in high-z discs that are gas-rich, and the role of mergers and accretions, which highz galaxies may have experienced at high rates -in maintaining/perturbing/destroying bars and b/p bulges.The role of the interstellar gas in the context of the generation/destruction of disc instabilities, such as bars (Bournaud et al. 2005) and spiral arms (Sellwood & Carlberg 1984;Ghosh & Jog 2015, 2016, 2022), has been investigated in past literature.In addition, the b/p bulges play a key role in evolution of disc galaxies by regulating the bar-driven gas inflow (e.g.see Fragkoudi et al. 2015Fragkoudi et al. , 2016)).Furthermore, bars can be weakened substantially (or even destroyed in some cases) as a result of minor mergers (Ghosh et al. 2021).Therefore, it would be worth investigating the b/p formation and evolution scenario in the presence of the thick disc and the interstellar gas and how likely they are to be affected by the merger events.

Summary
In summary, we investigated the dynamical role of a geometrically thick disc on the b/p formation and their subsequent evolution scenario.We made use of a suite of N-body models of thin+thick discs and systematically varied mass fractions of the thick disc and the different thin-to-thick disc scale length ratios.Our main findings are listed below.
-B/ps form in almost all thin+thick disc models with varying thick-disc mass fractions and for all three geometric configurations with different thin-to-thick disc scale length ratios.The thick-disc b/p always remains weaker than the thin-disc b/p, and this remains valid for all three geometric configurations considered here.
-The final b/p strength shows an overall decreasing trend with increasing thick-disc mass fraction ( f thick ).In addition, the b/ps in simulated galaxies with shorter thick-disc scale lengths form at earlier times and show a rapid initial growth phase when compared to the other two geometric configurations.Furthermore, we found a strong (positive) correlation between the maximum bar and b/p strengths in our thin+thick models.-For a fixed geometric configuration, the time interval between the bar formation and the onset of vertical buckling instability becomes progressively shorter with an increasing thick-disc mass fraction.In addition, for a fixed thick-disc mass fraction, models with shorter thick-disc scale length display shorter time delay between the bar formation and the onset of a buckling event when compared to the other two geometric configurations.-The final b/p length shows an overall increasing trend with increasing thick-disc mass fraction ( f thick ), and this remains valid for all three geometric configurations considered here.
In addition, for a fixed f thick value, the models with larger thick-disc scale lengths form a larger b/p structure when compared to the other two geometric configurations.Furthermore, the weaker b/ps are more extended structures (i.e.larger R b/p ).-The b/p structure changes appearance from being more X-shaped to being more box-shaped as the f thick values increase monotonically.This trend holds true for all three geometric configurations.Furthermore, the thin-disc stars are predominantly responsible for giving rise to a strong Xshape of the b/p structure.-Our thin+thick models go through a vertical buckling instability phase to form the b/p structure.The thin-disc stars display a higher degree of vertical asymmetry and buckling when compared to the thick-disc stars.Furthermore, the vertical asymmetry persists long after the buckling phase is over; the vertical symmetry in the inner region is restored relatively quickly, while the vertical symmetry in the outer region (close to the ansae or handle of the bar) is restored long after the buckling event is over.-The thin+thick models demonstrate characteristic signatures in the temporal evolution of different diagonal (σ zz /σ RR ratio) and off-diagonal (meridional tilt angle, Θ tilt ) components of the stellar velocity dispersion tensor, as one would expect if the b/p structure is formed via the vertical buckling instability.These kinematic signatures are more prominent when computed using only the thin-disc stars as compared to using only the thick-disc stars.
To conclude, even in the presence of a massive (kinematically hot) thick-disc component, the models are to able to harbour a prominent b/p structure formed via vertical buckling instability.This clearly implies that b/ps can form in thick discs at high redshifts and is in agreement with the observational evidence of the presence of b/ps at high redshifts (Kruk et al. 2019).Our results presented here also predict that at higher redshifts, the b/p will have a more boxy appearance than a more X-shaped one, which remains to be tested in future observations at higher redshifts (z = 1 and beyond).

Fig. 1 .
Fig. 1.Edge-on density distribution of all disc particles (thin+thick) at the end of the simulation run (t = 9 Gyr) for all thin+thick disc models with varying f thick values.Black dotted lines denote the contours of constant density.Left panels show the density distribution for the rthickS models, whereas middle panels and right panels show the density distribution for the rthickE and rthickG models, respectively.The thick disc fraction ( f thick ) varies from 0.1 to 1 (top to bottom), as indicated in the left-most panel of each row.The bar is placed along the x-axis (side-on configuration) for each model.The vertical black dashed lines denote the extent of the b/p structure in each case (for details, see the text in Sect.3.1.2).

Fig. 2 .Fig. 3 .
Fig. 2. Radial profiles of median of absolute value of distribution of particles in vertical (z) direction, |z|, (normalised by the initial value, z0 ), for thin-disc (left panels), thick-disc stars (middle panels), and total (thin+thick) disc stars (right panels), as a function of time (shown in colour bar) for the model rthickE0.5.Here, z0i denotes the initial value (used for the normalisation; for details, see text) where i = thin, thick, thin+thick, respectively.The thin disc b/p remains much stronger as compared to the thick disc b/p.

Fig. 5 .
Fig. 5. Strength of b/p, S b/p (top panel), and extent of b/p, R b/p(bottom  panel), calculated using thin+thick disc particles, at t = 9 Gyr shown as a function of thick-to-total mass fraction ( f thick ) and for different geometric configurations.With increasing f thick value, the b/ps progressively become weaker and larger in extent, and this trend remains true for all three geometric configurations considered here.The errors are calculated by constructing a total of 5000 realisations by resampling the entire population via a bootstrapping technique (for details, see text).

Fig. 6 .
Fig. 6.Temporal variation of b/p extent, R b/p , calculated using thin-, thick-disc, and thin+thick disc particles for model rthickE0.5.The b/p extent increases by a factor of ∼2 over the total simulation runtime.At t = 9 Gyr, the thick-disc b/p remains a bit larger than the thin disc b/p.The errors on R b/p are estimated by constructing a total of 5000 realisations by resampling the entire population via a bootstrapping technique (for details, see text).

Fig. 8 .
Fig.8.Distribution of mid-plane asymmetry, A Σ (x, z), in the edge-on projection (x−z-plane), computed separately for the thin (left columns) and thick (middle columns) discs, as well as thin+thick (right columns) disc particles using Eq.(4) at six different times (before and after the buckling event) for the model rthickE0.5.The bar is placed along the x-axis (side-on configuration) for each time-step.Black lines denote contours of constant density.A mid-plane asymmetry persists, even long after the model has gone through the buckling phase.

Fig. 9 .
Fig. 9. Bar-b/p strength correlation: distribution of all thin+thick models in maximum bar strength -maximum b/p strength plane.The maximum bar strengths are taken from Ghosh et al. (2023), whereas the maximum b/p strengths are determined from Eq. (1).The colour bar denotes the thick-disc mass fraction ( f thick ).Different symbols represent models with different geometric configurations (see the legend).The maximum strength of the bar correlates overall with the maximum b/p strength, and this remains true for all three geometric configurations.

Fig. 10 .
Fig. 10.Variation of b/p extent with varying thick-disc mass fraction.Top panel: temporal evolution of ratio of b/p length (R b/p ) and bar length (R bar ) for model rthickE0.5.Bottom panel: variation of ratio of b/p and bar length, calculated at the end of the simulation run (t = 9 Gyr), with thick-disc mass fraction ( f thick ).

Fig. 11 .
Fig. 11.Variation of time delay between the formation and onset of buckling instability, τ buck − τ bar , with thick-disc mass fraction ( f thick ), for all thin+thick models considered here.For a fixed geometric configuration, τ buck − τ bar becomes progressively shorter with increasing f thick values (for details, see text).
Fig. C.1.B/p extent, at the end of the simulation run (t = 9 Gyr), computed for the thin and thick discs, as well as for the thin+thick case, for all the thin+thick models considered here.Left panels show the distribution for the rthickS models, whereas middle panels and right panels show the distribution for the rthickE and rthickG models, respectively.Thin disc b/ps always remain shorter than the thick disc b/ps.The errors on R b/p are estimated by constructing a total of 5,000 realisations by resampling the entire population via a bootstrapping technique (for details, see Sect.3.1.2).

Table 1 .
Key structural parameters for the equilibrium models.
Quantifying the kinematic signature of teh vertical buckling instability and the subsequent b/p formation.Left panel: temporal evolution of vertical-to-radial velocity dispersion (σ zz /σ RR ), calculated within the b/p extent (R b/p ), (σ zz /σ RR ) 2 (t; R ≤ R b/p ), for thin (in blue), thick (in red), and total (thin+thick) disc (in black) particles, for model rthickE0.5.Right panel: temporal evolution of meridional tilt angle (Θ tilt ), calculated within R b/p using Eq. (