The secular growth of bars revealed by flat (peak + shoulders) density profiles

In many real galaxies, the bar’s major-axis surface-density profile has a multi-part structure. The innermost part of the profile is steep and (often) quasi-exponential. Beyond a certain radius, it flattens to a shallower profile (an ‘up-bending’ transition) before turning to a steeper slope (a ‘down-bending’ transition) further out. Finally, the steep outer profile sometimes becomes shallower again at or just beyond the end of the bar as it transitions to the disc profile (another ‘up-bending’ transition). We define the shoulder as the combination of the middle two regions: the shallow part of the profile plus the steep falloff, which together are the outermost part of the bar. Such a profile therefore exhibits a central peak, followed by the shoulder (‘peak + shoulders’, see Erwin et al., in preparation). We argue that profiles of this type are essentially the same as the ‘flat’ profiles identified by Elmegreen & Elmegreen (1985); we remind the reader that the slope of the shallow region need not actually be close to zero.

To illustrate this, Fig. 1 shows the major axis profile of three galaxies. Spitzer IRAC1 (3.6µm) profiles of NGC 1387 and NGC 4340 were retrieved from the Spitzer archive (PI K. Sheth, Program ID 10043); the image of UGC 9661 came from the Spitzer Survey of Stellar Structure in Galaxies (S⁴G; Sheth et al., 2010). Only NGC 4340 has the shoulder profile described above, NGC 1387 has a convex profile, and UGC 9661’s profile is rather noisy. The red vertical solid lines represent the inner and outer boundaries of the shoulder as found by SRA we have developed, which is described below. The red vertical dot dashed lines represent the centre of what we term the clavicle¹¹1In human anatomy, the clavicle is the collarbone, connecting the shoulder blade and the breastbone. We use it here to denote the flattest part of the shoulder structure., the centre of the flat portion of the shoulder, again as found by the SRA.

In the next subsection we describe the SRA in detail. Readers not interested in these details can skip to Section 3 which describes the models.

We develop a shoulder recognition algorithm that is automatic, responsive to the signal-to-noise ratio of the underlying image (or possibly variations caused by dust and/or star formation), has tunable parameters that enable it to match by-eye detections, and provides a natural quantification of shoulder parameters. Our method is equally suited to simulations and observations but in the description which follows, we focus on its application to simulations. A list of symbols we use in the quantitative portions of this work and their meanings are given in Table 1.

Rather than make any a priori assumptions about parametric profile components (for example by fitting multi-parameter exponentials or Sérsic profiles), we use a non-parametric approach. This has the advantages of being less subjective, and requiring no visual inspection to determine fitting ranges. It can also be relatively easily automated.

Since we are investigating shoulders within the bar, we measure the bar radius beforehand (see Section 3.3 for details). We then compute the logarithmic stellar surface mass density profile, $\log\Sigma_{\bigstar}(x)$, along the bar’s major axis, which we take as $|y|\leq 1$ kpc (see Section 3.4). In order to compare different profiles in a uniform manner, we normalise each profile to have values between 0 and 1 and denote the result as $\log\Sigma_{\bigstar N}$. We also normalise the $x$-axis to the bar radius $R_{\mathrm{bar}}$, so both axes are dimensionless.

Our method depends on the derivatives of this profile, so smoothing is essential. Having investigated a number of smoothing techniques, we settled on using a Butterworth lowpass filter of order 2 (Butterworth, 1930). The algorithm then obtains derivatives of the smoothed profile by the method of central differences; a successful implementation depends on sufficient noise reduction to ensure smoothly varying derivatives, but not so much that the smoothed profile fails to faithfully follow the major profile features.

We achieved this balance for the simulations by examining several profiles, with a particular focus on model 2 of Debattista et al. (2020) (see Section 3.1) at $t=5$ Gyr. We calculated the smoothed profiles with many combinations of smoothing extent and filter order. For each combination, we calculated the root-mean-square residuals between $\log\Sigma_{\bigstar N}$ and its smoothed version, as well as the dispersion in the differences between adjacent values of the derivative $\mathrm{d}\log\Sigma_{\bigstar N}/\mathrm{d}x$, and the number of extrema in that derivative (too many extrema implying insufficient smoothing). We selected the combination which gave a reasonable balance of noise reduction (with particular attention to the number of extrema in the derivative), and a faithful representation of the overall profile shape, and validated the results by visual inspection of the smoothed versus original profiles. While this part of the analysis is subjective, it is set once and held constant throughout.

We use the first derivative of the smoothed profile, $\mathrm{d}\log\Sigma_{\bigstar N}/\mathrm{d}x$, to test for the presence of shoulders (see also Lütticke et al., 2000). The process is illustrated in Fig. 2. For clarity, we show the profile for $x>0$, but the procedure is similar for $x<0$. The bar radius and its error are calculated using Fourier analysis of the stellar surface density projected onto the $(x,y)$-plane, described in detail in Section 3.3.

We find all extrema of the slope $\mathrm{d}\log\Sigma_{\bigstar N}/\mathrm{d}x$ which are within the bar (i.e., at $|x|<R_{\mathrm{bar}}$) and at distances from the centre $>0.2R_{\mathrm{bar}}$ (we impose the latter condition as we do not wish to mix shoulders with density features associated with a BP bulge, see Erwin & Debattista, 2016). We then identify the flattest such extremum – that is, the one with slope closest to zero (this marks the flattest part of the profile). Finally, we classify the profile as having a shoulder if that extremum is sufficiently flat: we define this as having a slope $<T$ (for $x<0$) or $>-T$ (for $x>0$). After visual inspection of many profiles, we settled on $T=0.4$ as a reasonable threshold. Hence we do not find shoulders if (i) there is no bar or (ii) there is no region of the profile meeting the flatness criteria with respect to $T$. Note that we consider a profile which turns beyond zero slope with an up-bend to be a shoulder profile (an example can be seen in the top left panel of Fig. 3).

If we do identify a shoulder, we name this minimum point in $|\mathrm{d}\log\Sigma_{\bigstar N}/\mathrm{d}x|$ as the centre of the clavicle (point [1] in Fig. 2). We denote the slope at the clavicle centre by $a$ – this represents the ‘flatness’ of the shoulder.

This analysis is repeated for $x<0$. By definition, both sides’ clavicle centres must be located at $|x|<R_{\mathrm{bar}}$. Hence, ring-like structures with overdensities outside the bar are not considered shoulders. If we do not detect a clavicle on both sides of $x=0$, the SRA deems the profile to have no shoulders. So we only recognize shoulder pairs. Moreover, since we recognize a clavicle as the point of the smallest absolute value of the slope, the algorithm does not detect more than one pair of shoulders.

If a profile has shoulders, we determine their extent; for this we use the radius of curvature of the smoothed profile, which is defined as:

and is shown as the orange dot-dashed line in Fig. 2. The only instance in the literature we have found which uses this parameter in a similar context is Lucatelli & Ferrari (2019), who use the curvature ($\mathcal{R}_{c}^{-1}$) to identify different components in the radial density profiles of observed galaxies. Foyle et al. (2008) analysed density profiles, identifying break radii, and discussed the use of profile derivatives to determine them (their Appendix B). They concluded that numerical methods were ‘clearly promising’, although they used visual estimates for simplicity.

We set the inner boundary of the clavicle to be the location of the minimum in $\mathcal{R}_{c}$ nearest to the centre of the clavicle, on the side closest to $x=0$, and its outer boundary to the corresponding minimum further out. We set the outer edge of the entire shoulder to be at the second outward minimum in $\mathcal{R}_{c}$.

In Fig. 2 we annotate key points in the quantification process in the natural order in which they are calculated: [1] the clavicle centre (thick red dot-dashed vertical line), at the minimum of $\mathrm{d}\log\Sigma_{\bigstar N}/\mathrm{d}x$ which is greater than $-T$, between $x=0$ and the bar radius; the slope at this point ($a$) is greater than $-T$ and so the SRA recognizes a shoulder; the closer $a$ is to zero, the flatter the shoulder; [2] the inner boundary of the clavicle, where $\mathcal{R}_{c}$ reaches its first inward minimum from the clavicle centre (thick purple dot-dashed vertical line); we define the inner boundary of the entire shoulder to be located here also; [3] the outer boundary of the clavicle, where $\mathcal{R}_{c}$ reaches its first outwards minimum from the clavicle centre (thick purple dot-dashed vertical line); [4] the outer edge of the entire shoulder, where $\mathcal{R}_{c}$ reaches its second outward minimum from the clavicle centre (thin red vertical line); we denote this as $R_{\mathrm{sh}}$; [5] a simple linear extension of the inner profile connecting inner and outer boundaries of the shoulder, constructed to estimate the excess mass contained within the shoulder (see below); we consider [4] to be the point at which the profile has reached the value where it would have been, had the shoulder been absent.

Based on these definitions, the extent of the clavicle is shown by the thick red double-headed horizontal arrow, and that of the shoulder by its green counterpart. Very thin shoulders (where the clavicle width $\leq 0.05L_{\mathrm{{bar}}}$) are rejected on the basis that these are likely to be local transient density perturbations rather than the extended morphological feature we are seeking.

We quantify shoulder ‘strength’ via the ‘excess mass’ in the shoulder. We define this quantity by extending a line from the shoulder’s inner to its outer boundary (line [5] in Fig. 2) and treating this as a notional ‘original’ profile; the difference between the actual shoulder profile and this exponential is the ‘excess’. Denoting the mass along the original profile between these points as $m_{o}$, and the mass along the line as $m_{l}$, we calculate the excess mass as $m_{e}=m_{o}-m_{l}$, which we can then normalise to the original mass between these points. Hence our measure of shoulder strength is the dimensionless $\mathcal{S}=m_{e}/m_{o}$. Stronger shoulders have higher fractional excess mass.

We define the errors on all calculated parameters as half the difference between the values for $x<0$ and $x>0$. As a proof of concept, we have run the SRA against the three bar major axis profiles shown in Fig. 1. The SRA correctly identifies the shoulders in NGC 4340, and does not recognize shoulders in the other two (NGC 1387 does not pass the first derivative slope threshold, and the shallowest first derivative slope for $x<0$ in NGC 9661 is outside the bar).

Our first set of models are pure $N$-body models with no gas or star formation. Several of these have been published already. From Debattista et al. (2020), we use models 2, 3, 4 and 5, which are baryon-dominated pure disc Milky Way-like models. From the same paper we also use the thin$+$thick disc model T1, as well as the dark matter-dominated model HD2. From Debattista et al. (2017) we use the thin$+$thick disc model T4.

We include an unpublished thin$+$thick disc model, T6. This is similar to T1 in that it has a thin and a thick disc of equal mass, both having scale length $R_{d}=2.4\mbox{$\>{\rm kpc}$}$. The main differences are in the geometric and kinematic parameters of the two discs. T6 has an initially thicker thick disc, with scale height $h_{z}=900\mbox{$\>{\rm pc}$}$ (versus $h_{z}=400\mbox{$\>{\rm pc}$}$ in T1), a central velocity dispersion $\sigma_{0}=140\mbox{$\>{\rm km\,s^{-1}}$}$ (versus $\sigma_{0}=90\mbox{$\>{\rm km\,s^{-1}}$}$ in T1) which declines exponentially with a scale length $R_{\sigma}=3.5\mbox{$\>{\rm kpc}$}$ (versus $R_{\sigma}=2.5\mbox{$\>{\rm kpc}$}$ in T1). The thin disc differs from that in T1 by being thicker, with $h_{z}=300\mbox{$\>{\rm pc}$}$ (versus $h_{z}=100\mbox{$\>{\rm pc}$}$ for the thin disc in T1).

We also include two unpublished bulge$+$disc models. These have been constructed using GalactICs (Kuijken & Dubinski, 1995; Widrow & Dubinski, 2005). They are based on the Milky Way-like models in Widrow et al. (2008), which are described also in Hartmann et al. (2014). The profiles of the bulges are given by

where $b$ is always set such that $R_{b}$ is the half-mass (effective) radius. In projection, this results in a Sérsic profile.

Model CB1 is based on the Toomre-$Q=1.75$, $X=3.5$ model of Widrow et al. (2008), with a more compact and more massive bulge. We use $n=1$, $p=0.44$, and $R_{b}=300\mbox{$\>{\rm pc}$}$. The density scale $\rho_{0}$ is not set directly but is set via the scale velocity (see Hartmann et al., 2014), which we set to $\sigma_{b}=350\mbox{$\>{\rm km\,s^{-1}}$}$. The comparable model in Widrow et al. (2008) has $n=0.85$, $p=0.36$, and $R_{b}=670\mbox{$\>{\rm pc}$}$, with scale velocity $\sigma_{b}=215\mbox{$\>{\rm km\,s^{-1}}$}$. Model CB1 is comprised of 2 million disc particles, 0.4 million bulge particles and 1 million dark matter particles.

Model CB2 is based on the Toomre-$Q=1.25$, $X=2.5$ model of Widrow et al. (2008), with another compact bulge, having $n=1$, $p=0.44$, $R_{b}=300\mbox{$\>{\rm pc}$}$ and $\sigma_{b}=400\mbox{$\>{\rm km\,s^{-1}}$}$. In contrast the comparable model in Widrow et al. (2008) has $n=0.85$, $p=0.44$, and $R_{b}=580\mbox{$\>{\rm pc}$}$, with scale velocity $\sigma_{b}=240\mbox{$\>{\rm km\,s^{-1}}$}$. Model CB2 is comprised of 0.9 million disc particles, 0.2 million bulge particles and 1 million dark matter particles.

Model PB1 is also a new model, which is comprised of a pseudobulge and an exponential disc built using the version of GalactICS presented in Deg et al. (2019). The halo of model PB1 is a Hernquist model, defined using the GalactICS parameters $\sigma_{h}=550\mbox{$\>{\rm km\,s^{-1}}$}$, $R_{h}=30\mbox{$\>{\rm kpc}$}$, $\alpha=1$, $\beta=4$. The main disc has $M_{d1}=4.31\times 10^{10}\>{\rm M_{\odot}}$, $R_{d1}=2.67\mbox{$\>{\rm kpc}$}$, and $z_{d1}=0.35\mbox{$\>{\rm kpc}$}$ while the pseudobulge is another exponential disc with $M_{d2}=4.77\times 10^{9}\>{\rm M_{\odot}}$, $R_{d2}=0.26\mbox{$\>{\rm kpc}$}$, and $z_{d2}=0.23\mbox{$\>{\rm kpc}$}$. Model PB1 is comprised of 2 million disc particles, 200,000 particles in the pseudo-bulge disc and 15 million dark matter particles.

The new model SD1 is built using a modified version of the GalactICS code that generates collisionless discs using a Sérsic surface density profile:

where $M_{d}$ is the disc mass, $R_{d}$ is the scale length, $n$ is the Sérsic index, and $\Gamma$ is the gamma function. Model SD1 has $5\times 10^{6}$ halo particles and $4.4\times 10^{6}$ disc particles. The halo uses the same parameters as PB1, while the Sérsic disc has $M_{d}=5.74\times 10^{10}\>{\rm M_{\odot}}$, $R_{d}=0.265\mbox{$\>{\rm kpc}$}$, $z_{d}=0.25\mbox{$\>{\rm kpc}$}$ and $n=2.05$.

We also include two models in which we suppress most secular bar growth by setting the halo in full prograde rotation (Debattista & Sellwood, 2000; Long et al., 2014; Collier et al., 2018). This results in large spin parameters of the haloes, $\lambda$, which are rare in cosmological simulations (Bullock et al., 2001). Since our goal is to contrast these models with their evolving and growing versions with the more typical, unrotating haloes, we are unconcerned by the fact that fully rotating haloes are very improbable in nature. The two models for which we have done this are model 2 (which becomes 2S, with $\lambda\simeq 0.084$) and model SD1 (which becomes SD1S, with $\lambda\simeq 0.091$).

All numerical parameters in the collisionless runs are the same as those in models 2, 3, 4 and 5 which are described in Debattista et al. (2020).

Models PBG1 and PBG2 are built using the version of GalactICS released in Deg et al. (2019), which builds equilibrium models with initially isothermal gas discs. These two models consist of 2 stellar discs, a gas disc, and a dark matter halo. For both models the main stellar disc has $M_{d1}=4.31\times 10^{10}\>{\rm M_{\odot}}$, $R_{d1}=2.67\mbox{$\>{\rm kpc}$}$, $z_{d1}=0.32\mbox{$\>{\rm kpc}$}$, and an exponential disc pseudo-bulge with $M_{d2}=5.67\times 10^{9}\>{\rm M_{\odot}}$, $R_{d2}=0.35\mbox{$\>{\rm kpc}$}$, $z_{d2}=0.3\mbox{$\>{\rm kpc}$}$. They both have halos with $\sigma_{h}=550\mbox{$\>{\rm km\,s^{-1}}$}$, $R_{h}=30\mbox{$\>{\rm kpc}$}$, $\alpha=1$, and $\beta=4$. Both models have a gas disc of mass $M_{g}=4.79\times 10^{9}\>{\rm M_{\odot}}$. The difference between the two models is that model PBG1 has a gas disc with scale-length $R_{g}=2.67\mbox{$\>{\rm kpc}$}$, while model PBG2 has a more extended disc, with $R_{g}=6.7\mbox{$\>{\rm kpc}$}$.

We also use the star-forming simulation described in Cole et al. (2014), Gardner et al. (2014) and, most extensively, in Debattista et al. (2017). Uniquely, this model starts out with a gas corona but no stars at all, so its evolution is free of any potential biases that might arise from inserting a disc of stars ab initio. We adopt Gardner et al. (2014)’s designation for this model: HG1.

Models PBG1 and PBG2 were evolved with changa (Jetley et al., 2008, 2010; Menon et al., 2015). Stars form out of gas with a $5\%$ efficiency, once the cool gas density exceeds 0.1 cm^-3 in a converging flow. Thermal and chemical turbulent diffusion uses the prescription of Shen et al. (2010), with a mixing parameter $D=0.03$. Gas and star particles have softening of $\epsilon=50~{}\mbox{$\>{\rm pc}$}$, including newly formed stars, while the dark matter particles have $\epsilon=100~{}\mbox{$\>{\rm pc}$}$. Supernova feedback couples 15% of the $10^{51}$ erg per supernova to the interstellar medium. We use a base time step of $\Delta t=2.5$ Myr with timesteps refined such that $\delta t=\Delta t/2^{n}<\eta\sqrt{\epsilon/a_{g}}$, where we set the refinement parameter $\eta=0.175$ and $a_{g}$ is the gravitational acceleration at a particle’s position. We set the opening angle of the tree-code gravity calculation to $\theta=0.7$. Gas particle timesteps also satisfy the condition $\delta t_{gas}=\eta_{courant}h/[(1+\alpha)c+\beta\mu_{max}]$, where $\eta_{courant}=0.4$, $h$ is the SPH smoothing length set over the nearest 32 particles, $\alpha$ and $\beta$ are the linear and quadratic viscosity coefficients and $\mu_{max}$ is described in Wadsley et al. (2004).

The simulations listed above do not include all the models we have studied. We include here models exhibiting interesting behaviours which help us understand the shoulder phenomenon. We emphasise that our simulations are isolated, so each model’s evolution is entirely secular.

We found no obvious differences between the collisionless models, the two models with gas and the pure star-forming model HG1, so in this study we do not compare collisionless models and models with gas. We examine model HG1 in Section 4.3.6.

As we are investigating profiles within the bar, we require measurements of the bar strength and radius. Since the bar is a bisymmetric deviation from axisymmetry, we define the bar strength, $A_{\mathrm{bar}}$, as the amplitude of the $m=2$ Fourier component of the stellar particle surface density distribution, projected onto the $(x,y)$-plane. Having rotated the models so the bar major axis is oriented along the $x$-axis, we calculate the azimuthal angle (with respect to the $x$ axis) $\phi_{k}$ of each stellar particle, of mass $m_{k}$, and then calculate:

Several methods have been developed for measuring bar sizes (e.g. Aguerri et al., 2000; Athanassoula & Misiriotis, 2002; Erwin, 2005; Michel-Dansac & Wozniak, 2006). We adopt the bar radius as the average of two calculations: the first is the (cylindrical) radius at which the amplitude of the $m=2$ Fourier component reaches half its maximum value after the peak. The second is the (cylindrical) radius at which the phase of the $m=2$ component deviates from a constant by more than $10^{\circ}$, with the uncertainty as half the difference. We refer to the resulting bar radius as $R_{\mathrm{bar}}$; the resulting average uncertainty is $10\%$. We calculate the time of bar formation, $t_{\mathrm{bar}}$, as the time when $\log(A_{\mathrm{bar}})$ changes from a positively sloped line to flat – i.e. when the instability has saturated and formed a steady bar.

We restrict our analysis to the bar by taking a cut along its major axis, $|y|\leq 1$ kpc. To show that our choice of cut in $|y|$ does not materially alter our results, we show in Fig. 4 the profile for various cuts in $|y|$, ranging from 500 pc to 2 kpc for model 2. Little additional profile information is gained by increasing the width of the cut while using too thin a cut only increases the noise. We therefore adopt $|y|\leq 1$ kpc for all models in this study. We call this the ‘major axis cut’ for ease of reference, and it always implies $|y|\leq 1$ kpc.

The bars in some of the models suffer the buckling instability (Raha et al., 1991; Martinez-Valpuesta & Shlosman, 2004; Smirnov & Sotnikova, 2019). Buckling manifests as a deviation from vertical symmetry, followed by a rapid increase in thickness, and generally the development of a box/peanut (BP) bulge (Bureau & Freeman, 1999; Debattista et al., 2004; Athanassoula, 2005; Bureau et al., 2006) in the inner part of the bar. We define the buckling amplitude as:

where $z_{k}$ is the $z$ coordinate of the $k^{th}$ particle. $A_{\mathrm{buck}}$ has dimensions of length.

We quantify the strength of a BP bulge following the method of Fragkoudi et al. (2017): we measure the maximum of the median of absolute heights ($\lvert{z}\rvert$) for the particles in radial bins of $0.1\mbox{$\>{\rm kpc}$}$. Normalising by the global median of absolute heights at $t=0$ ($|z|_{0}$) for each model makes comparison between models possible. So we define the BP strength, with tilde representing the median, as:

We measure the radial extent of the BP bulge, $R_{\mathrm{BP}}$, as the radius at which $|z|$ (and therefore $\mathcal{B}$) is a maximum. As a measure of vertical heating, we denote the global median $|z|$ for a model scaled to its value at $t=0$ as $\mathcal{Z}_{\mathrm{global}}$.

Fig. 3 (upper panels) shows examples of the original and smoothed profiles and their residuals for a selection of models. The residuals outside the bar are generally larger than those within, but this is not a concern as we are interested in the profiles within the bar. We also find that the central regions of highly concentrated models can have larger residuals; however these inner regions are not of concern either because the shoulders do not reside there.

Judging whether a profile has shoulders is somewhat subjective; as discussed earlier, the SRA uses the clavicle slope as a threshold. Consider the three profiles shown in the lower panels of Fig. 3. In all cases the profile slope changes at $|x/R_{\mathrm{bar}}|\sim 0.75$. At this point, the profile on the left has barely perceptible shoulders, the middle panel exhibits clear shoulders, while the right panel has strong shoulders. The slopes’ absolute values decrease as the shoulders strengthen. Thus the slope lends itself to shoulder identification and quantification – there is an absolute slope value above which we consider shoulders to be absent.

As an example from our simulations, the left panels of Fig. 5 show the logarithmic surface density plot in the $(x,y)$-plane for model 2 at $t=5\mbox{$\>{\rm Gyr}$}$ (the same model and time as in Fig. 4). The model has a strong bar; the shoulder profile is visible at $|x|\sim 7\mbox{$\>{\rm kpc}$}$ within the bar radius (indicated by the grey shaded regions), and is roughly symmetric about $x=0$. This is accompanied by an increase in the contour spacing along the major axis. The downwards bend occurs just after the end of the bar, as was found in the simulations of Athanassoula & Beaton (2006).

Note that we consider an ‘up-bending’ profile within the bar to be a shoulder profile. An example can be seen in the top left panel of Fig. 3 (model 5 at $t=5$ Gyr). The SRA recognizes this as a shoulder.

Outer rings are not considered shoulders. The right hand panels of Fig. 5 show model HD2 at $t=7.6$ Gyr. The flat overdensities at $|x|\sim 6$ and $9$ kpc are not shoulders, since they lie outside the bar. The surface density plot in the upper panel shows that these are rings. Moreover, shoulders are not the same phenomena as Freeman Type II profiles (Freeman, 1970; Erwin et al., 2008), down-bending disc breaks in the azimuthally averaged profiles which occur much further out in the disc than the bar.

We start our study of the simulation suite by describing the evolution of the bars in the models. Fig. 6 shows the evolution of the global bar and buckling amplitudes, $A_{\mathrm{bar}}$ and $A_{\mathrm{buck}}$. We group the simulations into those that strongly buckle (group B) and those that have weak or no buckling (group WNB). Five of the group B models (models 2, 2S, 3, 5 and T4) undergo a second major buckling. In most group WNB models, bar strength grows smoothly (top row, right two columns), but there are some with flat or declining bar strengths. We have explored other ways of grouping the models, including the bar ‘speed’ (fast or slow as determined by $\mathcal{R}$, the ratio of corotation radius to bar radius, with a separation at $\mathcal{R}=1.4$), and heights within and outside the bar. However we were unable to find any more insightful grouping than those that buckle, and those which do not (or do so very weakly). We also checked central concentration via the $C_{28}$ parameter (Kent, 1987), defined as $5\log(R_{80}/R_{20})$, with $R_{80}$ and $R_{20}$ being the cyclindrical radii containing 80% and 20% of the stellar mass, respectively. A high central mass concentration suppresses buckling (Berentzen et al., 2007; Iannuzzi & Athanassoula, 2015; Seo et al., 2019), so this is degenerate with the B/WNB grouping we adopt, with $\mbox{$\left<{C_{28}}\right>$}=2.78(3.78)$ for B (WNB) models at $t=0$. We therefore rely on the B and WNB groupings throughout.

Fig. 6 also shows the evolution of $\mathcal{Z}_{\mathrm{global}}$ (bottom row, dashed line); for model HG1, we normalise this to $\widetilde{|z|}$ at $t=3\mbox{$\>{\rm Gyr}$}$ since at $t=0$ we have no stellar particles, and by $3\mbox{$\>{\rm Gyr}$}$ we have a stable bar. This ratio increases with time in both groups. Each major buckling episode in group B is accompanied by a large increase in $\mathcal{Z}_{\mathrm{global}}$, as the buckling heats the disc vertically. In contrast, in group WNB models the increase in $\mathcal{Z}_{\mathrm{global}}$ is gradual. Table 2 presents an overview of some key evolutionary parameters.

We have shown that persistent shoulders develop as a part of bar growth, and are accompanied by BP bulges. We now examine quantitatively the relationship between bar/BP growth and shoulder evolution.

We now explore whether projection affects the recognition and observed properties of shoulders. We examine model 2 at various inclinations $i<60^{\circ}$ and relative bar position angles $\Delta$PA. The SRA recognizes shoulders for every combination of $i$ and $\Delta$PA, implying that shoulders, if they exist, would be observed at all inclinations up to $\sim 60^{\circ}$. For example, Fig. 17 shows the model and shoulder detection for $i=60^{\circ}$ at $t=5\mbox{$\>{\rm Gyr}$}$ (when the BP bulge and face-on shoulders are well established), for various bar position angles. Erwin & Debattista (2013) showed that moderately inclined galaxies with a BP bulge have isophotes with a boxy inner part from the thickened portion of the bar, accompanied by spurs. They also showed that when the bar’s position angle is beyond $\sim 50^{\circ}{}$ from the major axis, such a morphology may not be apparent. Fig. 17 shows that in such cases, the shoulders are still recognized.

In contrast, model 2S with a spinning halo does not manifest shoulders at any time step. Fig. 18 shows the model and shoulder detection for $i=60^{\circ}$ at $t=5\mbox{$\>{\rm Gyr}$}$, with the first column showing the face-on projection. Shoulders are not recognized at any combination of $i$ and $\Delta$PA.

Fig. 19 shows the ratios of inclined to face-on values for some of the observable shoulder parameters, when viewed in projection. Most parameters are not substantially ($\gtrsim 15\%$) affected. However, $a$ (i.e. flatness, top left panel) is strongly affected, and this parameter determines whether a shoulder is recognized. Values less than 1 indicate that a shoulder becomes flatter, and greater than 1 indicate that it becomes steeper, i.e. less shoulder-like. For higher $i$ and smaller $\Delta$PA ($\lesssim 30^{\circ}$), shoulder profiles become flatter. For higher $\Delta$PA, shoulder profiles appear less flat – potentially resulting in them disappearing altogether if they are intrinsically weak. Some caution is therefore needed when interpreting flat bar profiles in projection.

Erwin & Debattista (2016) presented the first direct detection of buckling bars, in NGC 4569 and NGC 3227 through examination of their isophotes. Using simulations, they noted that buckling galaxies are expected to exhibit spurs offset on the same side of the projected major axis, together with trapezoidal, rather than boxy, inner isophotes (their Fig. 2). Fig. 20 shows the surface brightness profile along the bar major axis (on the sky, not deprojected) of NGC 4569 ($i=69^{\circ}$, $\Delta\mathrm{PA}=26^{\circ}$). It is a system with shoulders within the bar (albeit with significant star formation contaminating the $x>0$ profile). Is it possible to have both the isophotal morphology of buckling and shoulders?

Fig. 21 demonstrates that it is indeed possible: we show model 4 shortly after peak buckling, at $t=4.2\mbox{$\>{\rm Gyr}$}$, rotated to the same orientation as NGC 4569.

The model presents the same isophotal morphology as NGC 4569, and the SRA recognizes shoulders in this projected profile. We find a similar signal in model 3, and in models 5 and T1 during their second bucklings (for a short time before the shoulders dissolve), but not their first. Thus, the models are consistent with those of NGC 4569, and the coexistence of shoulders and buckling does not necessarily imply a second buckling.

A key insight of this work is that shoulders are a manifestation of the bar’s secular growth. The simulations with prograde spinning haloes (Fig. 8) show that shoulders do not form if the bar is unable to grow. Since bars do not form with shoulders in place, shoulders are not part of the bar instability itself, but rather appear if the bar grows. Section 4.4.1 shows that shoulders end just outside the bar radius, with a relatively small variation in $R_{\mathrm{sh}}/R_{\mathrm{bar}}$ ($\sigma\sim 5\%$ of the mean), showing that the shoulder edge tracks the end of the bar rather closely.

Fig. 12 shows that shoulders typically become flatter as they evolve. This is not a general rule, however: models 3 and 4 undergo periods of weakening, and model PBG2 shows no strong growth trend. Furthermore, the data include episodes of steepening (i.e. becoming less shoulder-like) during secondary buckling and periods of spiral interference. Fig. 13 shows that for many models, in general the stronger and longer the bar, the stronger and flatter becomes the shoulder. This is consistent with Kim et al. (2015) who found that longer bars tend to show flatter profiles (defined by the Sérsic index of the bar profile) in their sample of 144 face-on barred galaxies.

However, Figs. 12 and 13 also show that it is not possible to use the observed profile flatness alone to determine the age of a bar, despite the general evolutionary trend towards flattening. The rate at which the slope flattens varies considerably between the models, and the decrease in slope is not monotonic in time (Fig. 12) or as the bar strengthens (Fig. 13). Furthermore, the dissolution of shoulders via secondary buckling and the profile’s subsequent return to exponential (grey markers in Figs. 12 and 13) bolsters the argument that exponential profiles are not necessarily indicators of young bars. We conclude that the flatness of the bar’s profile cannot be used as a chronometer, as suggested by Kim et al. (2015).

Buckling is preceded by a reduction in $\beta=\sigma_{z}/\sigma_{R}$, the ratio of the stellar vertical to radial velocity dispersions, to $\sim 0.5$ (Sellwood, 1996). We have confirmed that $\beta$ along the bar’s major axis declines before the first buckling. In our models, as was seen in simulations by Łokas (2019), the trigger value is closer to $0.6$. Buckling is followed immediately by a rise in $\beta$ in the buckled region. As the bar continues to grow post buckling, an increase in anisotropy (a renewed reduction in $\beta$) is seen in the shoulder region, as orbits of newly trapped particles become radially elongated (Figs. 10, 22). The increase is focused in the shoulder region since this is the location of highest $\sigma_{R}$. For the models which undergo a second buckling, $\beta$ eventually reaches $\sim 0.6$ in the shoulders once more, leading to buckling in this region, which destroys the shoulders (followed again by an immediate rise in $\beta$). So bar growth triggers both shoulder formation and eventually secondary buckling.

Persistent shoulders are often, but not always, accompanied by a BP. So a BP can be present without accompanying shoulders, perhaps for them only to emerge after a considerable time (e.g. model T1 which forms a BP $\sim 2$ Gyr before shoulders appear). In other models, particularly many of the WNB ones, shoulders appear before a BP. The spread in the ratio of shoulder edge to BP radius is twice that with the respect to the bar radius (Section 4.4.1); thus shoulders track the bar radius rather than the BP bulge radius. It appears that both shoulders and BPs – although both trace bar growth – form independently.

x₁ orbits are those primarily responsible for supporting the bar and have been the subject of many studies (see Contopoulos & Papayannopoulos, 1980; Sellwood & Wilkinson, 1993; Binney & Tremaine, 2008). Fig. 22 suggests that stars lingering in the looped region of x₁ orbits are responsible for the shoulders. Contopoulos (1988) found that x₁ orbits with loops appeared when the potential was strongly barred, and Athanassoula (1992) found the same in her models with high quadrupole moment. Valluri et al. (2016) discuss examples of x₁ orbits, one or two of which qualitatively resemble the shoulder-supporting morphology (their Fig. 4). They are not part of the backbone of the BP structure, being located near the ends of the bar (Gajda et al., 2016; Parul et al., 2020), which is near the outer edge of the shoulder structure (Fig. 2). Hence we expect BPs and shoulders to develop independently despite them both being impacted by bar growth.

As we have shown in Fig. 23, shoulder dissolution appears to transform these orbits, possibly into box orbits librating about the x₁ orbits. The buckling prevents the bar from growing (temporarily), and therefore renders it unable to capture additional particles into resonant looped x₁ orbits. Once the orbits already in the shoulder have evolved into more box-like orbits, the shoulder morphology is no longer supported and the bar profile becomes exponential once more, until fresh stars are trapped onto x₁ orbits.

We also note that the location of the inner limit of the shoulder in almost all models moves outwards as the bar evolves (Figs. 7, 8). This also supports the notion that the particles trapped into shoulder orbits over time evolve away from looped x₁ orbits, into more uniform box-like orbits (Łokas, 2019), thus eliminating the overdensity in that part of the bar as it continues to grow outwards.

We have used isolated galaxy simulations (16 collisionless $N$-body and three with gas and star formation) to study the outer regions of galactic bars, where the surface density profile along the bar major axis becomes shallow (or ‘flat’) and then breaks to a steep falloff, a pattern we term ‘shoulders’. Our main results are:

1. 1.

   Shoulders form as part of the bar’s secular evolution – they are a sign of a growing bar. They are not present when a bar first forms, and do not subsequently appear if a bar does not grow after formation (see Section 8.1).

2. 2.

   In many models, the strength and flatness of the shoulders increase as the bar evolves (although not monotonically). Most of our models are consistent with the observational findings of Kim et al. (2015) that stronger and longer bars have flatter profiles (see Section 4.4).

3. 3.

   Shoulders often – but not always – appear alongside a BP. Some models take a considerable time after BP formation before shoulders emerge and in some models, shoulders appear before a BP, so they are independent tracers of a growing bar. In non-buckling models, stronger and radially larger BP bulges are accompanied by stronger, flatter shoulders (see Section 4.4.3).

4. 4.

   Secondary buckling dissolves shoulders, either temporarily or longer term (depending on how effectively the bar captures additional material afterwards), returning the bar to an exponential profile without significantly weakening it. Strong spirals perturbing the bar can have the same effect. A thick disc can mask shoulders owing to significant mass being present at large heights. This destruction, the large variety of flattening rates of the bar as it evolves, and the fact that the shoulder growth is not monotonic, means that it is not possible to use the flatness of a bar’s profile in a simple way to determine its age. Notably, an exponential bar profile is not necessarily an indication of a young bar (see Sections 5, 8.1).

5. 5.

   For models with both BPs and shoulders, face-on shoulders are evident in projection, even though a ‘box$+$spurs’ morphology in the ($x,y$)-plane isophotes (a BP bulge indicator) may not be. The shoulder slope is strongly affected by projection, particularly at $i\gtrsim 45^{\circ}$. Caution is therefore urged when observing flat bar profiles in projection (see Section 6).

6. 6.

   We showed evidence hinting that shoulders are the manifestation of particles being trapped by the growing bar around looped x₁ orbits, where the time spent at apogalacticon results in the overdensity in the shoulder region. In time, these orbits transform, probably into librating boxes, so more material must be trapped for the shoulders to persist.

7. 7.

   We have verified that our conclusions are consistent with results from the fully self-consistent star forming model, and so are not peculiar to the collisionless models (see Section 4.3.6).