2.1 Governing equations

A calorically perfect fluid is injected into an unstratified quiescent ambient of the same fluid at lower temperature, from a circular orifice in an adiabatic wall. In order to model a flow with strong density variations but negligible compressibility, a low-Mach-number approximation of the compressible Navier–Stokes equation is used. This approximation, which retains all the effects of variable density in the convective terms, but discards the compressible dependency of density on pressure, was introduced by McMurtry et al. (Reference McMurtry, Jou, Riley and Metcalfe1986) for a study of non-buoyant jets in the limit of zero Mach number. It was then extended to include a buoyancy term by Nichols et al. (Reference Nichols, Schmid and Riley2007) and Chandler (Reference Chandler2010), and their formulation is used in the present investigation. The dimensional governing equations in this approximation are given by

where $\tilde{\unicode[STIX]{x1D70C}}$ , $\tilde{\boldsymbol{u}}$ , $\tilde{p}$ , $\tilde{T}$ denote the dimensional density, velocity, pressure and temperature, $\tilde{\unicode[STIX]{x1D70C}}_{\infty }$ is the ambient density, $g$ is the acceleration due to gravity, $\unicode[STIX]{x1D6FC}$ is the thermal conductivity, $C_{p}$ is the specific heat, $\unicode[STIX]{x1D707}$ is the dynamic viscosity, and $\mathscr{R}$ is the specific gas constant. Note that the pressure in this formulation is split into a thermodynamic component $p_{0}$ , which is constant throughout the flow, and a fluctuating hydrodynamic component $\tilde{p}$ . While the continuity and momentum equations (2.1a ) and (2.1b ) are of the same form as in the fully compressible case, the energy equation (2.1c ) simplifies to a simple advection–diffusion equation for temperature.

In dimensionless form, scaled with the nozzle radius $R$ , the inlet centreline velocity $\tilde{u} _{j}$ , the temperature difference $\tilde{T}_{j}-\tilde{T}_{\infty }$ between inflowing and ambient fluid, and the ambient density $\tilde{\unicode[STIX]{x1D70C}}_{\infty }$ , equations (2.1) become

A reduced temperature $T=(\tilde{T}-\tilde{T}_{\infty })/(\tilde{T}_{j}-\tilde{T}_{\infty })$ is used, and in all further calculations, this variable is expressed in terms of $\unicode[STIX]{x1D70C}$ by use of (2.2d ). The flow is characterised by the following parameters:

In the homogeneous limit $S\rightarrow 1$ , the state equation (2.2d ) prescribes $\unicode[STIX]{x1D70C}\rightarrow 1$ , and the system (2.2) is then identical to the Boussinesq equations used in Chakravarthy et al. (Reference Chakravarthy, Lesshafft and Huerre2015), provided $\unicode[STIX]{x1D70C}-1$ is taken to be of order $S-1$ ; the present formulation is therefore consistent with our earlier study. A rigorous derivation of the Boussinesq equations and a discussion on the underlying assumptions may be found in Tritton (Reference Tritton2012).

2.2 Base flow

In a cylindrical coordinate system $(r,\unicode[STIX]{x1D703},z)$ , the flow variables $q=(\unicode[STIX]{x1D70C},\boldsymbol{u},p,T)^{\text{T}}$ are split into steady and unsteady components as

The numerical tools used for the construction of the base flow, as well as for the linear perturbation analysis described in § 2.3, are very similar to those employed by Coenen et al. (Reference Coenen, Lesshafft, Garnaud and Sevilla2017), except that density variations are modelled as an effect of heating, rather than species mixing. Equations (2.2) are discretised with finite elements in FreeFEM $++$ , and a steady axisymmetric solution $\overline{\boldsymbol{q}}$ is obtained from Newton–Raphson iterations (Garnaud et al. Reference Garnaud, Lesshafft, Schmid and Huerre2013). The numerical domain is 80 radii long in the streamwise direction and 30 radii in the transverse direction. Iterations are performed until all flow quantities are converged to within $10^{-9}$ times their maximum values. At the inlet, $z=0$ , boundary conditions

are prescribed for the axial velocity $\overline{u}_{z}$ , the radial velocity $\overline{u}_{r}$ and the density $\overline{\unicode[STIX]{x1D70C}}$ . The initial shear layer momentum thickness is 10 % of the nozzle radius. All other boundary conditions are specified as

These boundary conditions are obtained from the kinematic constraints on the axis (Khorrami, Malik & Ash Reference Khorrami, Malik and Ash1989), and by imposing zero normal stresses at $r_{max}$ and $z_{max}$ , together with Dirichlet and Neumann conditions for the density.

As the objective is to characterise the role of buoyancy in the instability dynamics, the main parameters to be varied are the Richardson number and the density ratio. The ranges of parameters $10^{-4}\leqslant Ri\leqslant 10^{3}$ and $1.05\leqslant S\leqslant 7$ will be investigated. The effect of the Reynolds number above a value of 100 is found to be weak, and a standard value of $Re=200$ (in some cases $Re=500$ ) will be used, while $Pr=0.7$ is maintained throughout.

Figure 1.

Axisymmetric steady base flows obtained for $Ri=10^{-4}$ (a,b) and for $Ri=10^{3}$ (c,d). Axial velocity (a,c) and density (b,d) are shown in the ( $r,z$ ) plane.

Figure 2.

Variations of the centreline velocity along $z$ . (a) $Ri=10^{-4}$ ; (b) $Ri=10^{3}$ .

For a strong density ratio $S=7$ , and the two extreme values $Ri=10^{-4}$ and $Ri=10^{3}$ , steady base flow profiles of axial velocity and density are shown in figure 1. The flow at low $Ri$ is dominated by the momentum of the injected fluid, which diffuses radially with axial distance under the effect of viscosity. This is clearly the case of a jet, in a configuration where buoyancy has no noticeable impact on the base flow dynamics, despite the strong density variations. The flow at high $Ri$ , in contrast, is principally driven by the buoyancy force, as the injected momentum is negligibly weak in comparison. The fluid in this case is not pushed out of the orifice, but rather pulled out by buoyancy, forming a slender rising column around the axis (note the different radial scales in figure 1). This flow is characteristic of a plume, and it is often called a ‘lazy’ plume, as the momentum at its base is much lower than that of a self-similar profile, where momentum and buoyancy are in balance.

Between the two extremes shown in figure 1, the aspect of the base flow at different Richardson numbers changes gradually. Cases with $Ri<1$ will be denoted ‘buoyant jets’ in the following, as opposed to ‘plumes’ with $Ri>1$ . Vertical variations of the centreline velocity are shown in figure 2 for the same two configurations as in figure 1. The jet, at $Ri=10^{-4}$ , exhibits a short potential core region, where the centreline velocity is constant, followed by hyperbolic decay. The plume flow, at $Ri=10^{3}$ , accelerates progressively with vertical distance from the inlet, and the centreline velocity approaches asymptotically a limit value in the self-similar regime (Yih Reference Yih1988). Buoyant jets at low but finite Richardson number behave as plumes at large distances from the nozzle, when their excess momentum is sufficiently diffused. While a jet entrains ambient fluid only through momentum diffusion, entrainment into a plume tends to be much stronger due to its continuous production of axial momentum. The plume base flow presented in figure 1(c,d) is particularly marked by radial entrainment close to $z=0$ .

2.3 Linear stability problem

Infinitesimal perturbations of a steady base flow are sought with a global ansatz

The integer $m$ denotes the azimuthal wavenumber and $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{r}+\text{i}\unicode[STIX]{x1D714}_{i}$ is a complex frequency. Upon linearising the governing equations (2.2), and substitution of (2.10), the linear perturbation equations are found as

Differential operators in the above equations are written with a subscript $m$ in order to indicate that azimuthal derivatives are replaced with a factor $\text{i}m$ ; these operators are documented in appendix A. While the base flow is taken to be swirl free, $\overline{u}_{\unicode[STIX]{x1D703}}=0$ , the azimuthal perturbation velocity $u_{\unicode[STIX]{x1D703}}^{\prime }$ may in general be non-zero. Homogeneous Dirichlet conditions are imposed on $\hat{\boldsymbol{u}}$ and $\hat{\unicode[STIX]{x1D70C}}$ at the inlet, $z=0$ , and a homogeneous Neumann condition is prescribed for $\hat{p}$ . On the axis, depending on the azimuthal mode considered, appropriate boundary conditions as detailed in Khorrami et al. (Reference Khorrami, Malik and Ash1989) and Chakravarthy et al. (Reference Chakravarthy, Lesshafft and Huerre2015) are enforced. On the outer boundaries at $r_{max}$ and $z_{max}$ , no-stress conditions consistent with (2.9) are used.

The system (2.11) is cast into the form of an eigenvalue problem

According to the ansatz (2.10), the real part of the eigenvalue, $\unicode[STIX]{x1D714}_{r}$ , represents the oscillation frequency while the imaginary part, $\unicode[STIX]{x1D714}_{i}$ , is the growth rate of the perturbation. As before, the variable $\hat{T}$ is eliminated through the equation of state (2.11d ). So-called global eigenmodes are computed by resolving $\hat{\unicode[STIX]{x1D70C}}$ , $\hat{\boldsymbol{u}}$ and $\hat{p}$ in both $r$ and $z$ , such that spatial variations of the base flow and the perturbation quantities are accounted for without further limiting assumptions (Theofilis Reference Theofilis2003). The eigenvalue problem (2.12) is then solved with an iterative shift-invert Arnoldi algorithm, in the same way as in Garnaud et al. (Reference Garnaud, Lesshafft, Schmid and Huerre2013), with an accuracy close to machine precision.

In addition, a local analysis is performed in § 5, in order to identify the absolute mode in a parallel base flow (Huerre & Monkewitz Reference Huerre and Monkewitz1990). Perturbations (2.10) are then Fourier-transformed in $z$ , leading to the standard ansatz

3.1 Axisymmetric perturbations

Eigenvalues $\unicode[STIX]{x1D714}$ obtained for two configurations, $Ri=10^{-4}$ and $10^{3}$ , with otherwise identical parameters $S=7$ and $Re=200$ , are presented in figure 3 as black symbols. These are the two extreme jet and plume cases discussed in § 2.2.

Figure 3.

Global spectra of axisymmetric perturbations ( $m=0$ ), for the two configurations shown in figure 1. (a) $Ri=10^{-4}$ , jet case; (b) $Ri=10^{3}$ , plume case. $Re=200$ and $S=7$ are set identically for both cases. True eigenvalues (○) are compared to their counterparts (+) that are obtained when the buoyancy term in (2.11b ) is removed.

The jet, at $Ri=10^{-4}$ (figure 3 a), exhibits one unstable mode with $\unicode[STIX]{x1D714}=0.558+0.025\text{i}$ . When buoyancy effects are eliminated, by setting $Ri=0$ in the perturbation equations but still using the same base flow, eigenvalues shown as red crosses are obtained; visibly, the buoyancy term in the perturbation equations has no significant impact on the instability dynamics. This observation, as well as the overall appearance of the spectrum, is fully consistent with the non-buoyant and slightly buoyant results of Coenen et al. (Reference Coenen, Lesshafft, Garnaud and Sevilla2017), and many more details on the unstable eigenmode of pure jets are provided in their study.

The plume, at $Ri=10^{3}$ (figure 3 b), possesses five unstable modes, one being strongly dominant with $\unicode[STIX]{x1D714}=29.4+11.3\text{i}$ . When the perturbation Richardson number is set to zero in this configuration, all unstable modes vanish from the spectrum. It can be concluded that buoyancy plays a determining role for the instability of this plume, not only by setting up the base flow, but also by the coupling of density and velocity perturbations. It is also noted that the high- $Ri$ plume base flow is not subject to the non-buoyant instability that affects the low- $Ri$ jet.

Unstable eigenvalues of the plume take on significantly higher values, both in their real and in their imaginary parts, than those of the jet. This is clearly a result of the scaling with the inflow centreline velocity, which appropriately characterises a jet, but is less pertinent for ‘lazy’ plumes. A common scaling is employed here for the sake of consistency across the full range of Richardson number values; if only high- $Ri$ plumes were considered, a buoyancy-based velocity scaling would be more suitable.

The distinct nature of the instability modes of the jet and the plume flows is also apparent in the shape of their eigenfunctions. Figure 4 shows the spatial distribution of the axial velocity amplitude in the dominant modes for the two respective cases. At low Richardson number, the maximum perturbation amplitude is found on the centreline, 12 radii downstream of the nozzle. Perturbations are confined inside the jet column, as documented in more detail by Coenen et al. (Reference Coenen, Lesshafft, Garnaud and Sevilla2017). In the high- $Ri$ plume, the spatial eigenmode structure is very different: the maximum amplitude is located inside the mixing layer very close to the inflow, in the region where the density gradient in the base flow is maximal.

Figure 4.

Spatial distributions of axial velocity eigenfunctions $\hat{u} _{z}$ , corresponding to the most unstable modes of the two respective cases shown in figure 3. (a) Jet at $Ri=10^{-4}$ , $S=7$ and $Re=200$ ; (b) plume at $Ri=10^{3}$ , $S=7$ and $Re=200$ .

In the reference experiments by Cetegen & Kasper (Reference Cetegen and Kasper1996) and Bharadwaj & Das (Reference Bharadwaj and Das2017), a slightly different definition is chosen for the Richardson number, which corresponds to $2Ri/S$ in our nomenclature. Variations of the dominant eigenvalue in the present analysis are therefore presented as functions of $Ri/S$ in figure 5 in order to facilitate a comparison. The baseline case, with $Re=200$ , $S=7$ and inflow conditions (2.8), is represented by solid circles. Eigenvalues of this configuration display continuous variations both in the Strouhal number $St=\unicode[STIX]{x1D714}_{r}/(2\unicode[STIX]{x03C0})$ and in the growth rate $\unicode[STIX]{x1D714}_{i}$ . Strouhal number values are asymptotically constant in the low- $Ri/S$ regime, whereas they follow a power law at values $Ri/S>0.1$ . A regression fit yields the dependence $St=0.55(Ri/S)^{0.43}$ , which is in good agreement with experimental results in the range $1\leqslant Ri/S\leqslant 250$ : after conversion to the present definition of the Richardson number, the power law determined by Cetegen & Kasper (Reference Cetegen and Kasper1996) in this regime is given by $St=0.52(Ri/S)^{0.38}$ (shown as a line in figure 5 a), and the corresponding measurements of Bharadwaj & Das (Reference Bharadwaj and Das2017) for plumes from an orifice convert to $St=0.51(Ri/S)^{0.39}$ . A different scaling, measured as $St\propto (Ri/S)^{0.28}$ for $Ri/S>250$ by Cetegen & Kasper (Reference Cetegen and Kasper1996), is outside the parameter range considered here. The growth rate in the baseline configuration, shown in figure 5(b), increases monotonically with $Ri/S$ , and it is positive throughout.

Eigenvalues from three other flow configurations are included in figure 5 in order to assess the sensitivity of the instability with respect to the Reynolds number, to the density ratio and to the inlet velocity profile. With the standard profile (2.8), parameter combinations $Re=500$ , $S=7$ (open circles) and $Re=200$ , $S=4.5$ (squares) are chosen. The Strouhal number values in figure 5(a) are barely affected by these changes, and the growth rates in figure 5(b) follow a similar trend as in the baseline case. The less heated configuration (squares) is stable at $Ri<0.1$ , consistent with the analysis by Coenen et al. (Reference Coenen, Lesshafft, Garnaud and Sevilla2017).

Triangle symbols indicate results for a special case where $Re=200$ and $S=7$ are maintained, but the velocity inlet profile is changed to a parabolic pipe flow, while the density profile is still given by (2.8). This flow case is introduced in order to better approach the experimental conditions of Subbarao & Cantwell (Reference Subbarao and Cantwell1992) and Cetegen & Kasper (Reference Cetegen and Kasper1996), where the fluid exits from a nozzle as a developed laminar pipe flow. This change in the velocity profile has barely any effect on the Strouhal number across all $Ri/S$ values, but it does inhibit the global instability in the low- $Ri/S$ regime. The latter is consistent with the absence of self-excited behaviour at low $Ri$ in the experiments of Subbarao & Cantwell (Reference Subbarao and Cantwell1992) and Cetegen & Kasper (Reference Cetegen and Kasper1996), and in the simulations of Satti & Agrawal (Reference Satti and Agrawal2006).

Figure 5.

(a) Strouhal number and (b) growth rate of the global eigenvalue $\unicode[STIX]{x1D714}$ , as a function of $Ri/S$ . Legend: (●) $Re=200$ , $S=7$ ; (○) $Re=500$ , $S=7$ ; (▫) $Re=200$ , $S=4.5$ ; (▵) $Re=200$ , $S=7$ , with parabolic inlet velocity profile. All other configurations take the inlet velocity profile (2.8). (——) Power law from the Cetegen & Kasper (Reference Cetegen and Kasper1996) experiments, rescaled to match the present definition of $Ri$ .

As the effect of heating enters the problem both through the density ratio $S$ and through the Richardson number $Ri$ , the onset of instability is examined for independent variations of these two parameters. The main results of the local (see § 5) and global instability analyses are summarised in the state diagram in the ( $Ri,S$ ) plane shown in figure 6. The thin line delineates the neutral boundary separating locally convectively unstable inlet base flows (in white below the curve) from locally absolutely unstable ones (in blue and red above the curve). The thick neutral line separates the globally stable states (in white and blue below the curve) from the globally unstable states (in red above the curve). In the white area, convective instability prevails, and in the blue and red areas, absolute instability prevails. Along a diagonal line in the state diagram, the flow changes from convectively unstable to absolutely unstable to globally unstable, once the absolutely unstable region is large enough. For buoyant jets (low $Ri$ ) and plumes (large $Ri$ ) the same sequence takes place as $S$ increases at a given $Ri$ . Note that the globally unstable domain is reached much ‘sooner’ for plumes than for buoyant jets. The dip in the global stability neutral curve for $Ri$ of order unity signals a gradual shift from a shear-dominated instability to a buoyancy-dominated instability, as further discussed in § 4.

Figure 6.

Instability regimes in the $S$ – $Ri$ plane for the axisymmetric mode at $Re=200$ . The thick line denotes the global stability boundary, and the thin line denotes the boundary between local convective and absolute instability of the inlet profile.

3.2 Helical perturbations

The local analysis of plumes in the Boussinesq limit $S\rightarrow 1$ by Chakravarthy et al. (Reference Chakravarthy, Lesshafft and Huerre2015) concluded that absolute instability only occurs for helical perturbations, $m=1$ , whereas axisymmetric perturbations in that setting were found to be only convectively unstable. Although the global instability of axisymmetric eigenmodes in non-Boussinesq situations, as documented above, appears to be fully consistent with experimental and numerical observations of self-excited behaviour, the possibility of helical global instabilities remains to be explored.

Figure 7.

Global spectrum of helical perturbations ( $m=1$ ), for the configuration $S=7$ , $Ri=1$ and $Re=200$ . Results from two calculations are shown: (○) without absorbing layer; (●) with absorbing layer at $z>30$ .

Eigenvalues pertaining to helical instability modes are displayed in figure 7 for two different calculations, both for the same physical parameter setting $S=7$ , $Ri=1$ and $Re=200$ . One case, represented by blue unfilled circles, was computed with the same boundary conditions as all previous results. A branch of regularly spaced modes is seen to be unstable over the interval $0.8\leqslant \unicode[STIX]{x1D714}_{r}\leqslant 4$ . The features of this branch are very similar to several jet cases discussed by Coenen et al. (Reference Coenen, Lesshafft, Garnaud and Sevilla2017), as well as observations made in many different flow cases, especially in the boundary layer calculations by Åkervik et al. (Reference Åkervik, Ehrenstein, Gallaire and Henningson2008). In a recent study (Lesshafft Reference Lesshafft2017), the occurrence of such ‘arc branches’ is ascribed to the presence of spurious pressure feedback from the downstream end of a truncated flow domain, and it is predicted that artificial damping near the outflow should be effective in reducing the growth rate of such unphysical modes.

Indeed, if artificial damping is applied in an ‘absorbing layer’ (Colonius Reference Colonius2004) at $z>30$ , the growth rates of the arc branch modes decrease, and all modes recede to the stable half-plane if the damping is sufficiently strong. Such a case is represented by red filled circles in figure 7, where a damping term $-\unicode[STIX]{x1D706}(z)\hat{\boldsymbol{q}}$ has been added to the right-hand side of (2.12). The damping coefficient $\unicode[STIX]{x1D706}(z)$ ramps up from zero to 16, over the interval $30<z<50$ , according to equation (2.4) of Chomaz (Reference Chomaz2003). No unstable helical modes are found with this boundary treatment. In contrast, the same artificial damping has no discernible effect on the unstable eigenvalues for axisymmetric perturbations shown in figure 3. This behaviour is in full agreement with the arguments of Lesshafft (Reference Lesshafft2017), as $m=0$ perturbations are locally stable in the downstream flow region, whereas $m=1$ perturbations remain convectively unstable, as will be shown in § 5.

5.1 Axisymmetric perturbations, $m=0$

For the standard setting $S=7$ and $Re=200$ , and for Richardson numbers of $10^{-4}$ and $10^{2}$ , variations of the absolute frequency $\unicode[STIX]{x1D714}_{0}$ along the streamwise  $z$ direction are shown in figure 9. In the more extreme case $Ri=10^{3}$ , numerical difficulties led to unreliable results. In both configurations, the flow is found to be absolutely unstable ( $\unicode[STIX]{x1D714}_{0,i}>0$ ) over a streamwise interval of 6 or 7 radii downstream of the inlet. In the weakly non-parallel case of $Ri=10^{-4}$ , the real part $\unicode[STIX]{x1D714}_{0,r}$ shows moderate variations around a value of 0.5, in reasonable agreement with the global frequency $\unicode[STIX]{x1D714}_{r}=0.56$ as given in figure 5(a). In the strongly non-parallel setting $Ri=10^{2}$ , $\unicode[STIX]{x1D714}_{0,r}$ displays a variation between 1.3 and 80 over the absolutely unstable interval. This is not inconsistent with the global frequency $\unicode[STIX]{x1D714}_{r}=10.99$ , but it does not provide a means of predicting $\unicode[STIX]{x1D714}_{r}$ at leading order.

Nonetheless, the link between local absolute and global instability of axisymmetric perturbations is very plausibly confirmed by these results, both in the low- and in the high-Richardson-number regime. The neutral curve for the onset of absolute instability in the ( $Ri,S$ ) plane is reported in figure 6 for axisymmetric perturbations at $z=0$ . It is found that the transition from convective to absolute local instability at the inlet, with increasing $Ri$ and $S$ , occurs before global instability sets in. This is consistent with the common observation, both in model problems (Chomaz, Huerre & Redekopp Reference Chomaz, Huerre and Redekopp1991) and in slowly varying flow (e.g. Lesshafft, Huerre & Sagaut Reference Lesshafft, Huerre and Sagaut2007), that absolute instability must prevail over a sufficiently long streamwise region with sufficient growth rate in order to prompt a global instability.

5.2 Helical perturbations, $m=1$

The saddle point in the complex $k$ plane that is associated with helical absolute instability in the study of Chakravarthy et al. (Reference Chakravarthy, Lesshafft and Huerre2015) is also recovered in the analysis of the present inflow profiles. The absolute growth rate $\unicode[STIX]{x1D714}_{0,i}$ for $m=1$ perturbations at $z=0$ is displayed in figure 10(a) as a function of $S$ , for parameters $Ri=1$ and $Re=200$ . A unity Richardson number is chosen here for comparison with the analysis in § 3 of Chakravarthy et al. (Reference Chakravarthy, Lesshafft and Huerre2015), but higher values of $Ri$ lead to the same conclusions. The density ratio $S$ characterises the departure from the Boussinesq limit, and it is seen to have a very weak effect on the growth rate $\unicode[STIX]{x1D714}_{0,i}$ of the absolute helical mode.

The spatial variation of $\unicode[STIX]{x1D714}_{0,i}$ , over a short interval of $z$ adjacent to the inlet, is shown in figure 10(b) for the highly non-Boussinesq setting $S=7$ ; values beyond this streamwise region could not be obtained with sufficient confidence, due to numerical difficulties. The absolute helical growth rate at $S=7$ displays the same qualitative characteristics as the one described in Chakravarthy et al. (Reference Chakravarthy, Lesshafft and Huerre2015) for the Boussinesq limit: the growth rate declines quickly downstream of the inlet, perhaps asymptotically tending towards zero. It is not surprising that under these conditions the weak absolute instability of $m=1$ perturbations does not lead to a global flow destabilisation, as was found in § 3.2.

Figure 10.

Absolute growth rate $\unicode[STIX]{x1D714}_{0,i}$ of helical perturbations, for $Re=200$ and $Ri=1$ : (a) at the inlet, as a function of the density ratio $S$ ; (b) as a function of $z$ , for $S=7$ .