Stability of barotropic vortex strip on a rotating sphere

We study the stability of a barotropic vortex strip on a rotating sphere, as a simple model of jet streams. The flow is approximated by a piecewise-continuous vorticity distribution by zonal bands of uniform vorticity. The linear stability analysis shows that the vortex strip becomes stable as the strip widens or the rotation speed increases. When the vorticity constants in the upper and the lower regions of the vortex strip have the same positive value, the inner flow region of the vortex strip becomes the most unstable. However, when the upper and the lower vorticity constants in the polar regions have different signs, a complex pattern of instability is found, depending on the wavenumber of perturbations, and interestingly, a boundary far away from the vortex strip can be unstable. We also compute the nonlinear evolution of the vortex strip on the rotating sphere and compare with the linear stability analysis. When the width of the vortex strip is small, we observe a good agreement in the growth rate of perturbation at an early time, and the eigenvector corresponding to the unstable eigenvalue coincides with the most unstable part of the flow. We demonstrate that a large structure of rolling-up vortex cores appears in the vortex strip after a long-time evolution. Furthermore, the geophysical relevance of the model to jet streams of Jupiter, Saturn and Earth is examined.

We study the stability of a barotropic vortex strip on a rotating sphere, as a simple model of jet streams. The flow is approximated by a piecewise-continuous vorticity distribution by zonal bands of uniform vorticity. The linear stability analysis shows that the vortex strip becomes stable as the strip widens or the rotation speed increases. When the vorticity constants in the upper and the lower regions of the vortex strip have the same positive value, the inner flow region of the vortex strip becomes the most unstable. However, when the upper and the lower vorticity constants in the polar regions have different signs, a complex pattern of instability is found, depending on the wavenumber of perturbations, and interestingly, a boundary far away from the vortex strip can be unstable. We also compute the nonlinear evolution of the vortex strip on the rotating sphere and compare with the linear stability analysis. When the width of the vortex strip is small, we observe a good agreement in the growth rate of perturbation at an early time, and the eigenvector corresponding to the unstable eigenvalue coincides with the most unstable part of the flow. We demonstrate that a large structure of rolling-up vortex cores appears in the vortex strip after a long-time evolution. Furthermore, the geophysical relevance of the model to jet streams of Jupiter, Saturn and Earth is examined.
2018 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/ by/4.0/, which permits unrestricted use, provided the original author and source are credited.

Introduction
Jet streams are the prominent flow structures observed in atmospheric flows on Earth and gas planets such as Jupiter & Saturn [1,2]. It is of fundamental significance to understand how these jet structures appear, persist for a long time and become unstable in the longitudinal direction. The polar jet stream on Earth often intrudes into mid-latitudes and affects cyclonic storms in the atmosphere [3], thus predicting their course has become an important part of weather forecasting [4]. In this regard, the stability of jets is highly important.
The first step to understanding the complex phenomena of atmospheric flows, theoretically, is by introducing a simple mathematical model using a limited number of physical quantities. One such model is the barotropic model which considers an incompressible fluid of uniform density on the surface of a sphere. A remarkable contribution of this model to the understanding of zonal jets was provided by Rhines [5], where he found that, under a certain circumstance, the rotational effect and the energy cascade give rise to anisotropic east-west zonal band flows with a characteristic width known as the Rhines scale. The first numerical attempt to reproduce the two-dimensional turbulent flows with zonal jets in a plane with a rotational effect was made by Williams [6] using the barotropic model subject to random forcing. Subsequent numerical investigations have verified his numerical result in the case of a rotating sphere. Yoden & Yamada [7] have shown that the dynamics of unforced barotropic turbulence leads to a flow associated with a big circumpolar vorticity region, not zonal jets, as its quasi-final state. By using a homogeneous random forcing, Nozawa & Yoden [8] found that the barotropic flow generates, as a quasi-static two-dimensional turbulent state, strong circumpolar vortices accompanied by weak zonal jets when the rotation speed is moderate. Huang & Robinson [9] also constructed persistent latitudinal jet structures. Later, further numerical computation of the two-dimensional forced turbulence [10] showed that those zonal jets merge after a long-time evolution, and only a small number of zonal jets can survive as its asymptotic state. Note that other mathematical models such as shallow water models, general circulation models, quasi-geostrophic models and deep convection models have been used for further studies of jet streams on planets [1,11,12].
In view of the geophysical importance of zonal jets, we study the evolution of a latitudinal band region of a constant vorticity, called a vortex strip and investigate its stability. The linear stability of zonal jets on a β-plane has been studied by Obuse et al. [13]. They showed that these zonal jets are linearly unstable, and their instability triggers the merger of jets after a long-time evolution. However, the stability of zonal jets on a rotating sphere has not yet been investigated in detail so far. Our model for a vortex strip is based on the barotropic flow of an incompressible inviscid fluid on a rotating sphere. The linear stability and nonlinear evolutions of a vortex strip and a polar vortex cap on a non-rotating sphere were presented by Dritschel & Polvani [14,15]. In this sense, the purpose of this study is to extend their results to the case of a rotating sphere.
Here, we apply the vortex contour dynamics [16,17] to study the linear stability and to compute the nonlinear evolution of a vortex strip. A continuous vorticity distribution by the effect of the solid body rotation is approximated by a set of zonal band regions with piecewise constant vorticity. Embedding the vortex structure into these background vortex bands, we consider their evolution by tracking the boundaries of the vortex strip and the vortex bands. This approach has been used successfully in previous studies. For instance, the motion of a pair of vortex points [18] and a ring structure of N point vortices [19] on a rotating sphere were investigated through this approach.
The vortex strip can be regarded as a simple model of the shear layer, and several authors have presented numerical simulations of it. The evolution of a vortex strip in the plane, in connection with the Kelvin-Helmholtz instability, was studied by Pozrikidis & Higdon [20] and Roberts & Christiansen [21]. They found that the strip rolls up into a chain of vortices connected by so-called, braids. Recently, Bosler et al. [22] and Xiao et al. [23] reported numerical results for a thin vortex strip and zonal jet on a rotating sphere by using a Lagrangian particle/panel method and the RBF-vortex method, respectively. Most of these studies considered only thin vortex strips and did not examine thoroughly the dependence of the evolution of the strip on the physical parameters. Our main concern in this paper is the investigation of the stability of a vortex strip, i.e. unstable modes, growth rates and dependence on the rotation speed and the width of the strip, etc. Our aim is not the reproduction of the numerical results mentioned above. The numerical computation for vortex strips is conducted only up to the moderately nonlinear stage, but for a wide range of parameters, to compare with the linear stability. We thus rely on the framework of the vortex contour dynamics for the nonlinear evolution. We will therefore not be using advanced numerical methods such as the spectral method [24] and the semi-Lagrangian or Lagrangian methods [22,23,25,26]. The vortex contour dynamics is efficient to accomplish our purpose, and provides an advantage over the spectral method because the stability of the vortex strip is described not as the growth of spectra of its global spherical harmonic representation, but rather, in terms of the local deformation of contour curves. Owing to this property, we can examine not only whether the zonal jets are unstable, but also which wavenumber of local perturbations becomes the most unstable. Additionally, it allows us to provide a detailed description of how the instability of local zonal jets affects the flow of the whole sphere.
This paper is organized as follows. In §2, we describe the flow geometry of the vortex strip on a rotating sphere and explain how to approximate the flow in the framework of vortex contour dynamics. In §3, we present the linear stability analysis for vortex strips. Numerical computations of the nonlinear evolution of vortex strips are given in §4. The geophysical relevance of the vortex strip to the jet streams of Jupiter, Saturn and Earth is also discussed in §5. In the final section, we present our conclusion.

Flow descriptions (a) Distributions of vorticity
We consider a simple geophysical fluid model, the so-called barotropic model, in which the fluid is inviscid and incompressible, and of constant density on the surface of the unit sphere. The vorticity ω and the stream function ψ satisfy (2.1) Here, ∇ 2 Σ denotes the Laplace-Beltrami operator on the sphere in spherical coordinates (θ, ϕ) ∈ [−π/2, π/2] × [0, 2π ). The zonal and the meridional velocities, u and v, are derived from the stream function via the following geostrophic relation: Let us describe the basic flow of a vortex strip on a rotating sphere considered in this paper. Suppose that the boundaries of a vortex strip are located between two latitudinal lines θ = θ b 1 , θ b 2 on the sphere and that the vorticity constants are given by ω N for θ b 1 < θ < π/2, ω b for θ b 2 < θ < θ b 1 and ω S for −π/2 < θ < θ b 2 . As the vorticity distribution induced by the solid body rotation of the sphere becomes f (θ , ϕ) = 2Ω sin θ , where Ω is its angular velocity, the absolute vorticity of a barotropic vortex strip is then represented by A schematic of the vorticity distribution is shown in figure 1a. Owing to Gauss's theorem, the total vorticity over the sphere S 2 should be zero. As the integration of the vorticity induced by the solid body rotation is zero, i.e. S 2 2Ω sin θ dA = 0, we may consider only the integration of  3) between z b 1 = 2/3 and z b 2 = 1/2 on the rotating sphere with the angular speed Ω = 0.5. The thin (blue) curves and the dotted (red) curves represent the velocity fields induced only by the zonal strip, i.e. (2.3) with Ω = 0, and by the solid body rotation, respectively. the constant vorticities. Hence, the vorticity constants, ω N , ω b and ω S satisfy The constraint (2.4) is the same as that for the sphere without rotation considered by Dritschel & Polvani [14]. We examine relation (2.4) among the vorticity constants and the profile of the velocity induced by the barotropic vortex strip. If we assume that the vorticity constants of the northern and the southern polar regions are equal, say ω N = ω S = 1, then equation (2.4) yields which is always negative owing to −1 < z b 2 < z b 1  to the presence of the vortex strip, which results in the change of the zonal velocity from the eastward direction, above the strip, to the westward direction, below the strip. The dotted curve demonstrates that the solid body rotation gives rise to an eastward zonal flow. When those two flows are combined, as shown by the thick curve, we still observe the strong shear flow inverting the flow direction through the vortex strip. Suppose next that ω N = 1 and ω S = −1. Then from equation (2.4), it follows When the vortex strip is located in the northern hemisphere, i.e. 0 < z b 2 < z b 1 < 1, ω b is well defined and positive for any z b 1 and z b 2 owing to z b 1 + z b 2 = 0. Hence, it is uniquely expressed In other words, ω b is always positive and diverges as the width of the vortex strip tends to zero. Similarly, for the vortex strip in the southern hemisphere, −1 < z b 2 < z b 1 < 0, we have −∞ < ω b (r) < 0 for r > 1, which means that ω b is always negative and diverges as the width vanishes. For −1 < z b 2 < 0 < z b 1 < 1, ω b has a value between −1 and 1. Special care should be taken in the limiting cases where z b 1 → 0+ and can take any value in (−1, 1) in the zero limit of z b 2 = rz b 1 . Additionally, when z b 1 is zero, ω b = ω S = −1 holds regardless of the location of the lower boundary z b 2 , which results in no difference in the vorticity constant across the lower boundary. The same situation occurs for ω b = ω N = 1 when z b 2 = 0. Therefore, we exclude a consideration of the vortex strip in the case of z b 1 = 0 or z b 2 = 0 when ω N ω S = −1. Figure 1c shows the velocity profile of the vortex strip between z b 1 = 2/3 and z b 2 = 1/2 for ω N = 1 and ω S = −1 with Ω = 0.5. All the velocity fields are positive and the solid body rotation enhances the eastward velocity of the vortex strip. Finally, when ω N = −1 and ω S = 1, ω b is given by , which means that ω b is positive and diverges as the width goes to zero. A similar situation arises in the limiting case when z b 1 → 0+ and z b 2 → 0− as that for ω N = 1 and ω S = −1. Figure 1d shows the velocity profile of the vortex strip for ω N = −1 and ω S = 1 with Ω = 0.5. The velocity profile corresponding to the piecewise constant vorticity field has the opposite direction to the velocity induced by the solid body rotation.

(b) Contour dynamics formulation
We approximate the continuous vorticity distribution induced by the solid body rotation of the sphere with the M + 1 zonal bands with piecewise-constant vorticity separated by M latitudinal boundaries. To track the evolution of those boundaries, the flow on the surface of a sphere is embedded into three-dimensional space, using the three-dimensional Cartesian representation of those boundaries as described in [18,19]. We here comment on a different approach by Surana & Crowdy [27] which maps the boundaries to the complex plane through the stereographic projection and considers the analytic representation of the flows. This approach may be extended to the flow of a rotating sphere, which will be discussed later.
Let k (α, t), k = 1, . . . , M denote the boundary curves where α is the Lagrangian parameter along the curve and t is the time. Then the evolution equation of the curve k (α, t) is described by the following boundary integrals with respect to the M boundaries: where ω m represents the vorticity jump between the boundaries m−1 (α, t) and m (α, t) for m = 1, . . . , M. Note that 0 and M+1 are regarded as the north and south poles, respectively. The initial configuration is given by  with Let us note that each band with (2.9) has the same area because of the same height. The vorticity jump constant ω k of the kth band is chosen to be the average of vorticity at the band boundaries in the piecewise constant vorticity. For numerical convenience, we assume here that the two boundaries of a zonal vortex strip coincide with some of the M latitudinal band boundaries approximating the solid body rotation at the initial moment. Therefore, for the vortex strip, the vorticity constant ω k for each vortex band 1 ≤ k ≤ M + 1 is given by . . , M} are the indices of the band boundaries of the vortex strip. The vorticity jump ω k is then specified by in which the vortex strip is embedded in the region between the two band boundaries z = z k b 1 and We check how many boundaries are required to approximate the continuous vorticity distribution corresponding to the solid body rotation. Figure

Linear stability analysis
Let us explain the meaning of the stability of flows in our study. The zonal band configuration z k0 in (2.9) represents a steady state of equation (2.7). We consider the linear stability of the boundaries k (α, 0) for k = 1, . . . , M, giving an infinitesimal perturbation proportional to exp[i(mα − σ t)] in the latitudinal direction. The temporal growth rate σ depends not only on Ω and m, but also on (θ b 1 , θ b 2 ) for the vortex strip. We mainly pay attention to how the growth rate for a given mode m, i.e. the imaginary part of σ , is affected by the angular speed of rotation Ω and the width of the vortex strip h : Let us denote the kth boundary as where ζ k represents the perturbation added to the kth boundary. Substituting it into the z-component of (2.7) and linearizing it, we obtain where is the corresponding angular velocity. For each k, the value of Ω k is determined from the direct calculation of the steady configuration of the band boundaries (2.9). For example, as equation (2.1) is written as we have Ω k b 1 +1 for the zonal flow: The third equation of (3.3) can be found by using the constraint (2.4). Let us set the perturbation vector ζ = (ζ 1 , . . . , ζ M ) as whereζ = (ζ 1 , . . . ,ζ M ) denotes the Fourier amplitude vector of the perturbation for a given longitudinal wavenumber m, and c.c. represents its complex conjugate. Then, we obtain the eigenvalue problem for the amplitude vector, Here, A is an M × M matrix with the complex components where z > = max(z a , z b ) and z < = min(z a , z b ). The expression of A is analogous to that given by Dritschel & Polvani [18]. Solving the eigenvalue problem, we obtain M eigenvalues, which are denoted as σ (m) k for k = 1, . . . , M. We take the vorticity constants ω N = ω S = 1 and use M + 1 vortex band boundaries with M = 95 to approximate the solid body rotation. We examine the cases when the vortex strip is located between two band boundaries with even indices, namely k b 1  k ] = 0 for all k = 1, . . . , M and 2 ≤ m < 100 up to a numerical tolerance, then the flow is regarded to be neutrally stable. However, the flow is linearly unstable if there exists an eigenvalue with a positive imaginary part. As a matter of fact, in all the numerical results shown below, there exist only two eigenvalues with non-zero imaginary parts for each m when the flow is not linearly stable; one has a positive imaginary part and the other has a negative imaginary part with the same magnitude. We here note that the existence of the unstable eigenvalue with a positive growth rate means that not only the boundaries of the vortex strip but also all band boundaries become unstable as a whole. This will be discussed in more details shortly.
For a given pair of even (k b 1 , k b 2 ) with k b 1 < k b 2 , we calculate the eigenvalues for modes 2 ≤ m < 100. As we mentioned above, if the flow is linearly unstable, there exists only one eigenvalue with its imaginary part being positive. Accordingly, for notational convenience, we assume that σ  When the imaginary parts of the eigenvalues are zero for all m < 100, i.e. the flow is neutrally stable, no number is displayed. For Ω ≤ 5, we find in figure 3a-c that the flow is unstable for almost all pairs of (k b 1 , k b 2 ). For Ω = 60, the unstable region shrinks greatly as we see in figure 3d. Therefore, the rotation of the sphere gives a stabilizing effect on the flow with the vortex strip, which will be examined later. When we observe the stability diagram for Ω = 0.5 more closely, the most unstable mode gets larger as z k b 2 approaches the diagonal line along a fixed z k b 1 . This behaviour is similar to the result found by Dritschel & Polvani [14] for the non-rotating case. As the rotation speed Ω increases, the monotone behaviour of the most unstable mode along a fixed z k b 1 is broken and higher modes appear in the stability diagram. For Ω = 60, the most unstable modes are observed in a small band region below the diagonal line z k b 1 = z k b 2 . The growth rate corresponding to the most unstable mode m max in figure 3 is plotted in figure 4. We observe that, for fixed z k b 1 , the growth rate below the diagonal line becomes the largest and it decreases monotonically as z k b 2 goes southwards. This indicates that the flow tends to be less unstable as the width h increases. When we pay attention to the growth rates along the diagonal line, the peak of the growth rate is attained at the equator for all the cases of Ω. Note also that the small growth rates in the right low region in figure 4b for Ω = 0.5 increase slightly in figure 4c for Ω = 5. Figures 3a and 4a show the stability diagram and the associated growth rates for a small rotation speed Ω = 0.01, which is taken to see the behaviour in the vanishing limit of Ω. The growth rates in the right lower region are very small, which implies a convergence to 0 and the vortex strip thus becomes stable as Ω → 0. In the stability diagram, the modes in the lower right region correspond to small growth rates and are coloured with blue, and the modes in the remaining region are coloured with red. The shape of the red (unstable) region and the behaviour of the mode increasing along a fixed z k b 1 are similar to those of the non-rotating sphere [14]. Furthermore, from figure 4a, we find that, for Ω > 0, not too large, the strip of any width is unstable, whereas the strip with a large width is stable in the non-rotating case. Therefore, we conclude that the rotation produces an instability on the wide vortex strip when the rotation speed is small.
The eigenvector for the linearized equation gives rise to the unstable flow profile. Suppose that X 1 = 2.008 and σ 1 = 1.585, respectively, and the growth rates for m > 6 are mostly zero, we plot the eigenvectors for m = 2, 3, 4, 5 in figure 5a. All the cases in figure 5 show that the distribution for each m has a peak in the middle of the lower and upper boundaries of the vortex strip.
The unstable eigenvector tells us how each boundary evolves under the initial perturbation. To describe it, let us recall a basic theory of differential equations. Suppose now that A is a matrix whose distinct eigenvalues and eigenvectors are denoted by γ k with γ k = γ j (k = j) and u k for In order to discuss the limiting behaviour of the eigenvector in terms of the angular speed Ω, we choose two values, a large one Ω = 20 and a very small one Ω = 0.01. Figure 6 shows the distributions of the magnitude of the eigenvectors |X for 2 ≤ m ≤ 6. The vortex strip is located between z 16 = 2/3 and z 24 = 1/2. In figure 6a for Ω = 20, the eigenvectors for m = 2 and 3 are not shown, because their corresponding growth rates are small. The peak slightly moves to the lower boundary of the vortex strip at z 24 = 1/2 for all the modes m. For Ω = 0.01 in figure 6b, the peak of each distribution is located in the middle of the vortex strip, similar to the case of Ω = 0.5. This figure shows that the distribution of the eigenvector in the limit as Ω → 0 does not converge to that for the non-rotating sphere. Note that, in the non-rotating sphere, the eigenvector has only two components of the same magnitude at the two boundaries of the vortex strip, which means the most unstable boundaries. Consequently, these facts reflect a fundamental difference of the vortex strips in the rotating sphere from the non-rotating sphere. 1 | associated with the most unstable eigenvalue σ (4) 1 .
We check the convergence of the eigenvalues and their corresponding eigenvectors with respect to M. Figure 7a is the plot of the growth rates of the unstable eigenvalue σ where z (m) ∞ is the continuous limit of (3.1): Let us comment again on the choice of the number of discretizations M used in the stability analysis and following numerical computations. Although the result of M = 23 fairly resolves the distribution of the components of the eigenvector with respect to the profile with the peak, it does not provide good accuracy for other properties of the flow such as the stability diagram. The main reason for this is because of the inaccuracy in the velocity field as we see in figure 2.
We now observe the effect of the angular speed of the solid body rotation Ω on the eigenvalues and the eigenvectors. Figure 8 shows the growth rates of the unstable eigenvalues σ   for all m decreases with Ω, exhibiting an intermittent behaviour, and becomes zero when Ω exceeds 290. This indicates that the flow tends to be stable for large Ω, although such a fast rotation would be unphysical, in the geophysical sense. On the other hand, as Ω → 0, the growth rates converge to certain finite values. We have checked that these values agree with the non-rotating case [14].  Finally, we take the vorticity constants ω N and ω S as (ω N , ω S ) = (1, −1) and (−1, 1). Figure 9 shows the stability diagram of the most unstable modes for the two cases. The angular speed of the rotating sphere is set to Ω = 0.5 and the number of band boundaries is fixed at M = 95. According to the argument in §2a, for ω N = 1 and ω S = −1, the unstable region is split into three regions separated by the line z 48 = 0, and there is no unstable region in the middle. We find the monotone behaviour of the most unstable mode along a fixed z b 1 . For ω N = −1 and ω S = 1, all the cases considered for the regions are unstable and the monotone behaviour holds for the upper and lower left regions. The growth rates corresponding to the most unstable modes are also plotted next to the stability diagrams in figure 9. The maximum growth rates are attained when the vortex strip is located in the vicinity of the north and south poles. Figure 10 shows the plots of the distribution of the magnitude of the eigenvectors |X (m) 1 | associated with the unstable eigenvalues σ (m) 1 for 2 ≤ m ≤ 6 corresponding to the two cases of ω N and ω S in figure 9, in which the vortex strip is embedded between z 16 = 2/3 and z 24 = 1/2, and the most unstable mode is m max = 4. In figure 10, the eigenvalue profiles for m max = 4 have a peak in the middle of the upper and lower boundaries of the vortex strip for both cases. Let us take a closer look at the eigenvector of the most unstable mode. The value at the lower strip boundary z 24 is smaller than the value at the upper boundary z 16 . The components of the eigenvector between z 24 and z 36 have similar magnitudes and, thus they form an unstable band region. An interesting behaviour can, however, be found in the profile of the eigenvectors for ω N = 1 and ω S = −1. The distributions for the modes m = 2 and 3 have strong peaks near the equator. This suggests a possibility that the flow above the equator, not in the middle of the vortex strip, becomes the most unstable, which is not observed for ω N = ω S = 1. The distributions of the eigenvector for ω N = −1 and ω S = 1 also exhibit a similar behaviour.

Nonlinear evolutions of the vortex strip
We conduct numerical computations for the nonlinear evolutions of the vortex strip by solving the contour dynamics equation (2.7). We then compare the nonlinear evolutions with the linear stability analysis presented in the previous section. In all results of this section, we approximate the contour vorticity distribution corresponding to the solid body rotation with M = 95 band boundaries. We set the angular speed at Ω = 0.5. The initial configuration of the band boundaries is given by for a given mode m. Note that it corresponds to ζ k (α, 0) = 0.01 cos mα in (3.1), which does not coincide with the unstable eigenvectors discussed in the stability analysis. Each boundary is discretized by N = 512 points and the temporal integration is carried out by using the fourthorder Runge-Kutta method with the time step t = 0.005 or 0.01. The derivative of the boundary with respect to α in the integral of (2.7) is calculated through its Fourier series representation in three-dimensional Cartesian coordinates. To compare the nonlinear evolution with the linear stability analysis, let us consider the pth Fourier coefficientζ   where ζ k (α, t) denotes the perturbation added to the z-component of the kth boundary z k (α, t) in equation (3.1). We focus on the Fourier amplitudes of mode p, which is given by We here assume that the vorticity constants are fixed as ω N = 1 and ω S = 1 and then we change the location of the vortex strip and the mode m of the initial perturbation. We first take the boundaries of the vortex strip to be z 16 = 2/3 and z 24 = 1/2. Figure 11a shows the early time evolution of the Fourier amplitudes |ζ (±4) k (t)| for the most unstable mode m max = 4. Most of the Fourier amplitudes grow exponentially with the same rate predicted by the linear stability analysis. Figure 11b gives the snapshots of |ζ   figure 5a. Therefore, the nonlinear evolution at an early stage is in good agreement with the linear stability analysis. Figure 12 shows the evolution of the band boundaries from t = 1 to 2.2, where the two boundaries of the zonal strip are shown in blue, while the other band boundaries are shown in red. In figure 12, each panel shows every two contours viewed from the north pole (right) and the south pole (left). After t = 1.4, the boundaries in the vortex strip deform largely and evolve into a structure with rolling-up vortex cores. The number of the vortex cores is the same as the mode of the initial perturbation m max = 4. When we choose the other initial modes m = 2 and 3, we observe the similar results as shown in figure 13. The   exponential growth rates of the Fourier amplitude for both the cases are similar to those given by the linear stability analysis, although they are not shown here. The band boundaries in the middle of the vortex strip become unstable and develop into large rolling-up vortex cores in the strip. The number of the vortex cores again coincides with that of the initial mode. Now we consider the case when the vortex strip is located in the southern hemisphere between z 64 = −1/3 and z 72 = −1/2. In this case, the width of the vortex strip is unchanged. The most unstable mode for this case is m max = 5 as shown in figure 5a. Figure 14a shows that the band boundaries in the middle of the vortex strip become unstable and they evolve into a structure with five rolling-up vortex cores. The evolution and the snapshots of the Fourier amplitudes at the early stage are also plotted in figure 14b,c, respectively. The exponential growth rates of the spectra are almost the same as the linear growth rate and the snapshots |ζ  the middle of the vortex strip, whose profiles are also similar to the eigenvector |X Let us consider vortex strips with larger widths when the exponential growth rate is small as in figure 4a. When we consider the vortex strip between z 40 = 1/6 and z 56 = −1/6 for the initial perturbation of the most unstable mode m max = 3 as shown in figure 15, the exponential growth rate of the Fourier amplitudes still agrees with the linear stability analysis. Their snapshots |ζ (±3) k (t 0 )| at t 0 = 1, 2 and 3 in figure 15c show that the peak location is very near the equator, which is consistent with the eigenvector |X 1 | in figure 5c. After the long-time evolution of the zonal jet, the band boundaries in the middle of the vortex strip form three large rolling-up vortex cores at t = 5. We can however find an inconsistency between the linear stability analysis and the nonlinear evolution when we consider the vortex strip with a larger width. The nonlinear evolution of the vortex strip between z 16 = 2/3 and the equator z 48 = 0, and the growth of the Fourier amplitudes for the initial perturbation of the most unstable mode m max = 2 are shown in figure 16. We observe a large deformation of the band boundaries in the middle of the vortex strip as in figure 16a, but no rolling-up vortex cores appear by t = 8. Figure 16b demonstrates a deviation of the growth rates of the Fourier amplitudes |ζ   We now change the signs of the vorticity constants ω N and ω S , fixing the two boundaries of the vortex strip at z 16 = 2 3 and z 24 = 1 2 . We find more complex behaviours in this case. Figure 17 shows the snapshots of the Fourier amplitudes and the nonlinear evolutions of the vortex strip for the initial perturbation of the most unstable mode m max = 4 when the vortex constants are given by (ω N , ω S ) = (1, −1) and (ω N , ω S ) = (−1, 1). In both cases, the growth rates of the Fourier amplitudes agree with the exponential growth rate given by the linear stability analysis, and the snapshots |ζ  change of the band boundaries to elliptic curves, which cannot be identified clearly when plotted on the sphere. The result of figure 10 demonstrates that the band boundaries, above the equator, away from the vortex strip may deform largely.
Let us make a technical remark regarding the rotation speed. In the numerical computations, we have considered only one rotation speed Ω = 0.5. In fact, if we take a sufficiently large rotation speed, the computation is not sustainable for a long time. We have found that when the rotation speed is high, the boundaries do not stay on the sphere and tend to deviate from the surface of the sphere, because they are represented in three-dimensional Cartesian coordinates in the vortex contour dynamics model. This causes the breakdown of the computation at a late time. To overcome this problem, a different mathematical formulation for the equations in the spherical coordinate system is required [28]. Another possible method is to map the zonal band boundaries, approximating the solid body rotation of the sphere, to circular concentric boundaries in a plane by using a conformal stereographic projection as in [27]. We leave these works for future studies.

Geophysical relevance to planetary jet streams
We examine the geophysical relevance of the model for a vortex strip on a rotating sphere to jet streams in planets. The parameters are estimated using the Rossby number R o = U/2ΩL, which is a non-dimensional parameter with Ω, L and U being the angular speed, the radius of a planet and a representative latitudinal wind speed in the jet stream, respectively. The planetary parameters of Jupiter, Saturn and Earth are found in [2,4,12] and are summarized in table 1. As Ω = 0.5 and L = 1 in our model, it follows from for Jupiter is estimated as u J = 0.005. Similarly, we obtain the representative speed u S = 0.005 for Saturn and u E = 0.06 for Earth. Note that the stability analysis for Jupiter and Saturn yields the same results owing to u S = u J . In our model, the velocities of the boundaries of the vortex strip at z = z b 1 and z b 2 are given by Suppose ω N = ±ω S ; then we determine the vorticity constant ω N as follows: , (5.1) where u = u b 1 − u b 2 is equivalent to the representative wind speed defined above. Figure 19 is the stability diagrams for the most unstable modes and their corresponding growth rates for the parameters relevant to Jupiter/Saturn and Earth. The continuous vorticity distribution corresponding to the solid body rotation is approximated by M = 95 band boundaries. The rotational speed is Ω = 0.5 and we assume ω N = ω S . For the vortex strip between z b 1 and z b 2 with b 1 = 4p and b 2 = 4q (1 ≤ p < q < 23), the value of the vorticity constant ω N is obtained from (5.1) for u = u J = u S = 0.005 (Jupiter/Saturn) and u = u E = 0.06 (Earth).
In figure 19, we first pay attention to the stability of vortex strips with small width, i.e. for b 1 = 4p and b 2 = 4(p + 1), which are aligned just below the diagonal line of the stability diagrams, because those localized jet structures are often observed in many planetary atmospheric flows. In Jupiter/Saturn case, the most unstable modes for polar vortex strips tend to be low with large growth rate, whereas those for equatorial strips become high with small growth rate. On the other hand, in the case of Earth, the most unstable modes for polar vortex strips are lower than those for equatorial strips, which is similar to the case of the Jupiter/Saturn. However, their corresponding growth rates are different: the growth rates for vortex strips in the mid-latitude of the northern hemisphere are larger than those for polar vortex strips. Note that the stability diagrams for ω N = −ω S are similar to those in figure 19a,b for 0 < z b 2 < z b 1 < 1 and are not shown here. In figure 19a, we also find a stable region along the line z b 1 = z b 2 for Jupiter/Saturn case, where all perturbations up to m = 100 wavenumber become neutrally stable. Another stability region appears along the region of z b 2 = −0.5 as we see in figure 19b for the Earth case. It is not evident whether the presence of such a stable region has a certain geophysical significance, which should be investigated in the future.
We next examine the stability of specific vortex strips corresponding to strong polar jets observed at about 75-80 • N of Saturn [12]. As the width of vortex strips is so narrow that their boundaries do not coincide with those of the band boundaries for M = 95, we use a finer approximation M = 383. As the boundary locations 1 ≤ b 1 ≤ b 2 ≤ M are chosen arbitrarily, we have checked the three cases of (b 1 , b 2 ) = (3, 5), (4,6) and (6,8)    rapidly. In [12], a hexagonal circumpolar jet stream of Saturn was numerically demonstrated from a general circulation model. Here, we also confirm the hexagonal instability of the polar jet from the barotropic vortex contour dynamics. as the width gets smaller, regardless of the locations of the vortex strip. When the narrow vortex strip is located around the equator, it has the largest growth rate. For a small Ω ≈ 0, the rotation generates an instability on the wide vortex strips, unlike the non-rotating case. However, as Ω increases, the vortex strip tends to be less unstable and the most unstable mode gets higher. For a sufficiently large Ω, the vortex strip becomes neutrally stable, which reveals a stabilizing effect owing to the rotation of the sphere. We have performed the numerical computations of the nonlinear evolution of vortex strips on the rotating sphere. When the width of the vortex strip is small, the exponential growth rate of the Fourier amplitudes in the nonlinear evolutions at early times fits well with that predicted by the linear stability analysis. Furthermore, the profiles of the Fourier amplitudes are similar to the eigenvector corresponding to the unstable eigenvalue, and the band boundaries in the middle of the strip are the most unstable among all band boundaries. Hence, the unstable eigenvector suggests which boundary becomes the most unstable in the nonlinear evolution. After a longtime evolution, the vortex strip develops into a varicose large structure with rolling-up vortex cores, whose number is the same as the mode of the initial perturbation. However, a deviation of the growth rate of the nonlinear evolution from the linear stability analysis is apparent when the width is large. This may be due to small growth rates and excitation of some nonlinear effects.

Conclusion
We have also studied the stability of the vortex strip when the polar vorticity constants have opposite signs, (ω N , ω S ) = (1, −1) and (ω N , ω S ) = (−1, 1). The linear stability analysis shows that the exponential growth rate corresponding to the most unstable eigenvalue becomes larger as the width is smaller; however, the most unstable location of the vortex strip is in the polar regions, and not around the equator. The eigenvector similarly provides information which band boundary becomes the most unstable, but it has different and complex profiles from those for ω N = ω S = 1, where there is a peak outside of the vortex strip and an unstable band region. In the nonlinear evolution, the vortex strip of a small width in the northern hemisphere evolves into a structure of rolling-up vortex cores in both cases; however, the rolling-up direction for (ω N , ω S ) = (1, −1) is opposite to that for (ω N , ω S ) = (−1, 1). Importantly, boundaries far away from the vortex strip can be unstable, for a different choice of mode of the initial perturbation. This phenomenon suggests that the polar jet stream yields an instability of atmospheric flows around the equator under certain circumstances. We conclude that the solid body rotation gives significant influences on the stability and the evolution of the flow of the vortex strip.
We have examined the geophysical relevance of the model for a vortex strip to the planetary jet streams such as Jupiter, Saturn and Earth. We find that the stability diagram and growth rates of Jupiter and Saturn are very different from those of Earth. For Jupiter and Saturn, the polar vortex strips have the largest growth rate, while for Earth, the strips in the mid-latitude of the northern hemisphere are the largest. This behaviour might be related to appearance of the jet stream near the north pole for Saturn and in the mid-latitude for Earth. Furthermore, the linear stability analysis for Saturn demonstrates the hexagonal instability for the vortex strip at 74-76 • N, which is in a reasonable agreement with the observation of the polar jet stream [1]. Therefore, our simple model based on vortex contour dynamics not only describes qualitative behaviours but also provides quantitative prediction for the stability of jet streams.
Data accessibility. This article has no additional data. Authors' contributions. S.S. and S.K. conducted the linear stability analysis and T.S. computed the nonlinear evolutions numerically. All the authors gave their final approve for publication.
Competing interests. We declare we have no competing interests.