On the spontaneous symmetry breaking of eastward propagating dipoles

The spontaneous symmetry breaking of weak eastward propagating dipoles is revealed from high-resolution numerical simulations in an equivalent-barotropic quasigeostrophic beta-plane model. The evolution of initialized Larichev–Reznik dipoles is found to depend on the initial intensity, characterized by the dipole drift speed divided by the maximum Rossby wave phase speed. We found that weak dipoles lose symmetry and eventually disintegrate due to growing perturbations. This instability is attributed to a critical, linear growing mode that is extracted from the potential vorticity anomaly. The growth rate of perturbations is found to decrease with the increasing dipole intensity

Dynamics of coherent vortices in rotating and stratified fluids is a classical research topic that remains active and well-motivated by the need to understand ubiquitous vortex phenomena in geophysical fluids such as oceans and atmospheres.It is generally recognized that vortices are crucially affected by the large-scale gradient of the ambient potential vorticity (PV) (the beta-effect).The dynamics of both monopolar and multipolar vortices have been investigated in the literature through the application of various observational, laboratory, analytical, and numerical methods.In this Letter, we describe a new type of beta-related instability of steadily propagating dipoles consisting of a closely packed pair of counter-rotating vortices.
Wide families of steady-state dipoles originate from the Lamb-Chaplygin solution to the two-dimensional (2D) Euler equations. 1,2This solution is stable in the presence of 2D perturbations; 3 however, 2D unstable modes have been found in the presence of viscosity. 4As an analog to the Lamb-Chaplygin dipole, a stationary dipole solution was obtained on the barotropic b-plane by Stern, 5 and even more generally, a zonally drifting equivalent-barotropic solution was obtained by Larichev and Reznik. 6In the barotropic limit when the Rossby radius is infinite, the Larichev-Reznik dipole (hereafter, LRD) is only capable of eastward drift; otherwise, it can drift in either zonal direction.A review of results surrounding LRDs can be found in Ref. 7. Numerical experiments suggest that the eastward propagating LRDs remain stable during limited times of integration, for example, when exposed to short-wave perturbations, 8 weak friction, 9 and topographic perturbations. 10][13] On the b-plane, it was shown that the fate of unsteady f-plane dipoles depends on their initial direction of propagation and intensity. 14In particular, initially strong eastward propagating dipoles decelerate and adjust into a steadily propagating state, whereas initially weak eastward propagating dipoles breakup into two well-separated monopolar vortices drifting in the opposite, westward direction.Another study with a focus on LRDs tilted away from the zonal direction found that their ability to adjust to steady eastward propagating states depends on the initial tilt and intensity. 15These studies motivate systematic research on stability of various dipole configurations, and our study is within this category.
Here, we demonstrate that LRDs can experience spontaneous symmetry breaking and become unstable.We obtained our results by utilizing superior numerical techniques and higher spatial resolution, surpassing those used in previous studies.
To begin our investigation, we consider the upper-ocean equivalent-barotropic quasigeostrophic (QG) model on the beta-plane 16,17 as The above QG approximation assumes smallness of the Rossby number: e ¼ L Û = f 0 ( 1, where L and Û are natural length and velocity scales, respectively, and f 0 is the leading-order Coriolis parameter.We nondimensionalized the model such that where variables without hats are dimensionless.Substituting Eq. (3) into Eq.( 1), the leading-order balance yields the flow evolution governance, where w ¼ wðx; yÞ is the velocity streamfunction, q ¼ r 2 w À Sw is the corresponding PV anomaly (PVA), S ¼ ð L= Rd Þ 2 ; b ¼ bL 2 = Û ; ¼ = L Û , and JðÁ; ÁÞ denotes the Jacobian operator.The quadratic invariants of energy and enstrophy can be obtained by integrating over the employed doubly periodic domain and they are conserved for the inviscid case ð ¼ 0Þ.Moreover, the kinetic moment, is a linear invariant even in the presence of viscosity.The dipole steady states of interest are exact solutions in the coordinate system translating uniformly with the dipole center ðx c ðtÞ; y c ðtÞÞ.Therefore, they satisfy (4) under the transformation ðx; yÞ !ðx À x c ; y À y c Þ, which reads as where c x and c y are constant zonal and meridional propagation speeds of the dipole, respectively.Spatially localized PV solutions to (7), with ¼ c y ¼ M ¼ 0, are often referred to as equivalent-barotropic "modons"; these objects maintain their shape by avoiding radiating linear (Rossby) waves.The well-known Larichev-Reznik solution 6 is incapable of meridional drift, but can propagate in both eastward and westward directions, with the only restriction being that the westward drift speed must exceed the maximum zonal speed of the linear waves.These steady states are derived from the relative vorticity, where ðr; #Þ are polar coordinates, c ¼ c x , p 2 ¼ b=c þ S, and k is a positive constant satisfying the following nonlinear equation: where J l and K l denote order-l Bessel function and modified Bessel function of the first kind, respectively.This equation has an infinite number of solutions for k, but for the purposes of our analysis, we consider the dipole solutions with the smallest k.The corresponding streamfunction can be expressed as where We carried out numerical simulations with doubly periodic boundary conditions, on the uniform 4N Â N ¼ 8192 Â 2048 rectangular computational grid in a ð4L y ; L y Þ ¼ ð60; 15Þ domain.The numerical model is described in Ref. 18, and its numerical convergence properties can be found in Refs.19 and 20.The advantage of the numerical algorithm is that the advanced CABARET advection scheme allowed us to use coarser grids for the same numerical accuracy and minimized numerical dissipation.We evaluated the numerical solutions by comparing how well the invariants are conserved for different resolutions and small values of viscosity .Since there were no significant differences for a range of considered (see the supplementary material), we present solutions without explicit viscosity (but with small and unavoidable numerical viscosity).
Eastward propagating LRDs with speed c ¼ 0.1 were imposed in the model at initialization.We focus on scales L $ R d that are reflective of oceanic vortices; hence, in nondimensional form S ¼ 1 and v R ¼ b.In Ref. 14, strong and weak dipoles were defined by keeping b fixed and varying c such that the corresponding dipole intensities were c=b ¼ 1 and c=b ¼ 0:25, respectively (since the nondimensional maximum Rossby wave speed is v R ¼ b).In our case, we fix c and consider weak dipoles with b ¼ 0:4; b ¼ 0:3, and b ¼ 0:2 as well as a strong dipole with b ¼ 0:1.
During the evolution of a weak dipole, we observe that for a significant amount of time, it appears to steadily propagate eastward.However, beyond a certain point, it begins to oscillate and decelerate, up until t ¼ t s , which is the time of zonal drift decaying to zero and the eddy pair starting to separate.Using parabolic interpolation, we found t s $ 215 for b ¼ 0:4 and larger t s $ 327 for smaller b ¼ 0:3 [changes in domain size do not significantly change the value of t s , Fig. 1 (Multimedia view)].For times larger than t s , the dipole disintegrates into two westward propagating monopolar vortices, similar to those described in Ref. 14.Eventually, these vortices disintegrate into the background, but their further analyses go beyond the scope of our study.
To analyze the symmetry breaking of the dipole flow, we implement boundary conditions for a no-flow-through zonal wall along the dipole axis and draw comparisons between solutions.With the wall, the LRD remains stable (see Fig. 2), which suggests that one potential cause for the spontaneous symmetry breaking is the exchange of PV between the vortex pair.We checked the hypothesis of PV exchange, resulting from the spontaneous symmetry breaking, by considering the PVA integrals in the upper and lower domain halves.Since total PVA is conserved (also in our numerical solutions), we expect these quantities to evolve in opposite directions, if PV exchange occurs.Indeed, we found that the upper-/lower-half integral experiences small but significant decrease/increase over time, but verifying whether this exchange causes disintegration of the dipole remains for future studies.
We compared the time series for several integral invariants ( 5) and ( 6) obtained using N ¼ 1024 and N ¼ 2048 (the latter ones are in Table I) to evaluate the accuracy of our numerical model.If we consider integral energy and enstrophy conservation, as well as (material) extremum PV conservation, for solutions obtained with N ¼ 1024, then, over the interval t ¼ 50 À 150, the following changes in values are noted: DE ¼ 0:609%; DZ ¼ À0:545%, and DP max ¼ À1:40%, respectively.In comparison with solutions obtained using N ¼ 2048, we find that DE ¼ 0:381%; DZ ¼ À0:0300%, and DP max ¼ À0:575%.This means that the energy conservation is roughly 37% better on a 8192 Â 2048 uniform grid, while the conservation of enstrophy and extremum PV is approximately 95% and 59% better, respectively.When approaching t ¼ t s , we see a faster decrease in energy and other invariants, though this is more controlled for higher resolution.Such energy and enstrophy losses are consistent with increased dissipation of the advection scheme on steep gradients emerging on violent growing perturbations.To better understand the development of oscillatory motion for weak dipoles, we extracted the coordinates of the positive and negative dipole extrema, (x 1 , y 1 ) and (x 2 , y 2 ), respectively, using local 2D-parabolic interpolation.The dipole center is (x c , y c ), where x c ¼ ðx 1 þ x 2 Þ=2 and y c ¼ ðy 1 þ y 2 Þ=2.Evolution of the extrema and center are described by the corresponding trajectories.We get physical insight from time series of the following quantities: c ¼ dx c =dt (zonal drift speed), Dx ¼ x 1 À x 2 (the characteristic of the dipole tilt), and Dy ¼ ðy 1 À y 2 Þ=ðy 1 ð1Þ À y 2 ð1ÞÞ (the separation of the dipole partners, normalized to start at unity) (see Fig. 3).
For the 0.25 intensity dipole, we see small and slowly growing oscillations up until t $ t s À 50, beyond which c decreases down to zero at t ¼ t s (t s $ 215 for b ¼ 0:4).In comparison with the 0.33 intensity dipole, we see similar behavior, but with some time delay: t $ t s À 70 (t s $ 327 for b ¼ 0:3).The dipole wobbling is apparent from the time series, and we see that the symmetry breaking begins before the partner separation, as illustrated by the dipole center oscillations.As t approaches t s , oscillations in Dx become more pronounced, indicating further development of the asymmetry, and Dy clearly indicates the increasing separation of partners.For both weak-dipole intensities considered, at t ¼ t s the meridional separations between the dipole extrema are similar, suggesting some critical distance reached before the dipole breakup.
Due to the base flow symmetry, the PV field in the domain can be uniquely divided into two components, following (Ref.4): where q A is A-component, which has even PV relative to the zonal axis, and q S is S-component, which has odd PV relative to the zonal axis.The interpretation of these components is the following: S-component perturbation leads to symmetric deformation of the vortices around the zonal axis, whereas A-component perturbation leads to the antisymmetric deformation of the vortices.This is equivalent to the reversed notation in Ref. 21.
The energy E A and enstrophy Z A of the A-component begin to grow from zero due to the spontaneous breaking of symmetry [Fig.4(a)].After the period of small initial oscillations and during an intermediate period of linear growth (t 1 < t < t 2 ), Z A grows at a nearly constant rate r Z .For the weakest dipole (b ¼ 0:4), we estimated r Z $ 0:07; t 1 $ 100; t 2 $ 180.As the value of b decreases, we found that r Z decreases.
In Fig. 4(b), the characteristic length scale derive from the Acomponent, , indicates that initially antisymmetric perturbations appear at grid scale, grow toward the scale of the most unstable mode until t $ t 1 , and then remain around L A $ 0:2.Nonlinear effects that develop after that will be investigated in a separate paper.The dipole evolution for other b values look qualitatively similar, but oscillations before the linear stage are larger, periods of oscillations in the linear stage are longer, and the follow-up decrease in L A is slower for stronger dipoles.These differences are most notable for the 0.5 intensity dipole (b ¼ 0:2), which maintains themselves for significantly longer and decelerate slower, before the partners separate (t s $ 562).
Figure 5 presents q A snapshots in the linear growth stage for b ¼ 0:4, when small oscillatory motion becomes pronounced in Fig. 3. From the panels in Fig. 5, one can see a nearly time-periodic spatial structure for the growing A-perturbations, analogous to the growing A-mode shown in Fig. 12 of Ref. 4. The corresponding period is T $ 16, which is estimated from phase evolution, and the growth rate is estimated from the increase in A-perturbation amplitude over the period: r ¼ ðln ðq m ð159Þ=q m ð143ÞÞ=T $ 0:035, where q m ¼ maxðq A Þ. Therefore, r Z $ 2r, which is consistent with quadratic invariants having a growth rate of 2r.In addition, similar estimates  We can associate the A-perturbation pattern with a linear instability mode of the dipole by seeking q A of the following form: q A ¼ Aðx; yÞ cos ðxtÞ À Bðx; yÞ sin ðxtÞ ½ e rt ;

Physics of Fluids
where A and B are the real and imaginary parts of the mode, respectively, x ¼ 2p=T is the oscillation frequency, and r is the exponential growth rate.Fitting this equation with only two snapshots from Fig. 5 allowed us to obtain A and B (see Fig. 6), such that Eq. ( 13) produces all the other snapshots consistent with Fig. 5.This makes us to conclude that the spatiotemporal structure of the growing A-component can be interpreted as a critical linear-instability mode with a period of T $ 16.The spatial structure of A and B (Fig. 6) is substantially more complicated than the 2D viscous linear instability A-mode of the Lamb-Chaplygin dipole (Ref.4, Fig. 8).The real part of the latter Amode corresponds to meridional displacement of the dipole, and the imaginary part corresponds to tilting of the dipole axis.The wake formed behind the drifting dipole is seen both in the real and imaginary parts of the A-mode, and it is due to viscosity.In contrast, the double quadrupolar structure of the PV perturbation of the LRD A-mode (Fig. 6) describes elongation and contraction of counterrotating PV cores in the dipole partners (see supplementary animation in the supplementary material).Note that emergence of quadrupolar instability in circular f-plane equivalent-barotropic vortices has been reported in Ref. 22, but here the novel aspect is the A-mode synchronization of two such instabilities in the dipole vortex partners.
In the LRD A-mode, there is also a symmetric (about the zonalaxis) wake.Its formation may be a consequence of Rossby wave radiation mechanism in contrast with the viscous mechanism forming the wake in the Lamb-Chaplygin dipole.Moreover, the wake associated with the LRD A-mode fans out (perhaps, due to radiating meridional Rossby waves), whereas the Lamb-Chaplygin viscous wake thins out, perhaps, due to the absence of Rossby waves.In the Lamb-Chaplygin dipole, the instability is attributed to interaction between the wake and dipole. 4The formation of the wake affects the dipole motion that amplifies the wake, so that this positive feedback mechanism supports the instability growth.Our weak dipoles also produce wakes (Fig. 2); however, analysis of these wakes and their feedback roles is left for future studies.
Next, we considered a strong initial LRD corresponding to c=b ¼ 1 (Fig. 7, Multimedia view), and the solution propagates nearly steadily, like a weak dipole in the presence of a wall (Fig. 2).With our current methodology, it is difficult to assess whether this solution is actually stable or just has a very small growth rate, because significantly longer time integrations need longer domains (to avoid artificial interaction with the wake in periodic domain) and finer grid resolution (to reduce long-term effects of numerical viscosity).Regardless, with this in mind and considering the characteristic differences for different values of b [Fig.4(b)], there could exist a critical dipole intensity threshold in the range 0:5 < c=b < 1, where dipoles with intensities less than this critical value experience the symmetry breaking phenomena, while those with greater intensities remain steadily propagating.
Finally, our LRD instability results are new and disagree with previous studies; therefore, we additionally checked the methodology by replicating the results for an f-plane dipole in the presence of b-effect 14 (see the supplementary material).We also imposed a no-flow-through zonal wall in the middle and found that it has no effect on the f-plane dipole dynamics, which is consistent with the absence of growing Acomponent perturbations.
In summary, we have shown that weak eastward propagating dipoles on the b-plane in an equivalent-barotropic QG model experience the phenomena of spontaneous symmetry breaking.We initialized our simulations with symmetric steady-state LRD flow, and despite such dipoles appearing to maintain a steady translation in the short term, they begin to decelerate and exhibit increasing oscillatory motion.At t ¼ t s , the vortex cores reach a critical separation and then quickly split and start propagating in the opposite, westward direction; thus, the dipole breaks apart.Comparative simulations, in the presence of a no-flow-through zonal wall along the LRD axis, demonstrate the physical significance of the partner separation, while a resolution comparison of invariant conservation laws proves the convergence of our numerical solutions.In Fig. 4, we clearly see the intermediate stages of linear growth, with growth rates ranging from r $ 0:035 for 0.25 intensity dipoles and r $ 0:02 for 0.5 intensity dipoles.
Furthermore, our experiments for stronger dipoles suggest looking for some critical dipole intensity threshold, in the interval of 0:5 < c=b < 1 that may divide eastward propagating dipoles into unstable and stable ones.However, it may be the case that very strong dipoles have such slow growth rates that instabilities are simply not picked up by the limited time span of our numerical solutions.Since it is still unclear from our analysis whether these stronger dipoles are stable or unstable, this motivates the study of the stronger dipole evolution in future works.
Given the limited scope of this short paper, we could only present a few examples of dipole instability.[25] See the supplementary material for animations of the dipole propagation and additional figures that support the results of this study.

FIG. 2 .
FIG. 2. Nondimensional contour plots of the dipole PV field at: t ¼ t s À 2T [panels (a) and (d), where T is the period of oscillation related to the growing critical mode discussed later], t ¼ t s [panels (b) and (e)], and t ¼ t s þ 2T [panels (c) and (f)].Top panels correspond to solutions in the absence of a zonal wall, and bottom panels correspond to solutions in the presence of a zonal wall.Solutions obtained using b ¼ 0:4 are shown in a fixed frame of reference.

FIG. 3 .
FIG. 3. Time series for dipoles with b ¼ 0:4 (red) and b ¼ 0:3 (blue), where t s is the time of dipole disintegration (t s $ 215 for b ¼ 0:4 and t s $ 327 for b ¼ 0:3): panel (a) shows the zonal propagation velocity, panel (b) shows the meridional center coordinate, panel (c) shows the zonal extrema difference, and panel (d) shows the meridional extrema difference (adjusted to start at unity).

FIG. 4 .
FIG. 4. Time series for dipoles with b ¼ 0:4 (red), b ¼ 0:3 (blue), and b ¼ 0:2 (black): panel (a) shows the logarithm of enstrophy for the A-component, and panel (b) shows the corresponding length scale L A ¼ ffiffiffiffiffiffiffiffiffiffiffiffi E A =Z A p .Time intervals differ since weaker dipoles disintegrate faster.
LETTERscitation.org/journal/phf for stronger dipoles result in approximately the same period but a slower growth rate: r $ 0:027 for b ¼ 0:3 and r $ 0:02 for b ¼ 0:2.As for the S-component, we see insignificant change in comparison with the initial pure symmetric state.

TABLE I .
Invariant estimates for various values of t, all obtained with N ¼ 2048 grid resolution.