Eigenvectors, Circulation and Linear Instabilities for Planetary Science in 3 Dimensions (ECLIPS3D)

Context. The study of linear waves and instabilities is necessary to understand the physical evolution of an atmosphere, and can provide physical interpretation of the complex flows found in simulations performed using Global Circulation Models (GCM). In particular, the acceleration of superrotating flow at the equator of hot Jupiters has mostly been studied under several simplifying assumptions, the relaxing of which may impact final results. Aims. We develop and benchmark a publicly available algorithm to identify the eigenmodes of an atmosphere around any initial steady state. We also solve for linear steady states. Methods. We linearise the hydrodynamical equations of a planetary atmosphere in a steady state with arbitrary velocities and thermal profile. We then discretise the linearised equations on an appropriate staggered grid, and solve for eigenvectors and linear steady solutions with the use of a parallel library for linear algebra: ScaLAPACK. We also implement a posteriori calculation of an energy equation in order to obtain more information on the underlying physics of the mode. Results. Our code is benchmarked against classical wave and instability test cases in multiple geometries. The steady linear circulation calculations also reproduce expected results for the atmosphere of hot Jupiters. We finally show the robustness of our energy equation, and its power to obtain physical insight into the modes. Conclusions. We have developed and benchmarked a code for the study of linear processes in planetary atmospheres, with an arbitrary steady state. The calculation of an a posteriori energy equation provides both increased robustness and physical meaning to the obtained eigenmodes. This code can be applied to various problems, and notably to further study the initial spin up of superrotation of GCM simulations of hot Jupiter.


Introduction
The study of the influence and propagation of waves in planetary atmospheres is often performed using models with several simplifications, most notably the assumption of a zero or zonallysymmetric and constant initial zonal flow (e.g. Kasahara & Qian 2000), or restriction to a beta-plane solution (e.g. Lindzen 1967). Such simplifications allow analytical prediction of the key wave mechanisms, and in some cases a complete understanding of their structure (even in the mathematical sense; see Matsuno 1966).
Despite the simplifications, such studies have provided significant insight into atmospheric dynamics. Wheeler & Kiladis (1999) for example demonstrate that the propagation of waves can be linked to convective motions in the atmosphere of the Earth, and the resulting description matches analytical theories (e.g. Vallis 2006;Holton 1992) remarkably well. Baroclinic and barotropic instabilities are also known to have an impact on the circulations of planetary atmospheres (see Williams 2003).
However, for more complex mean flows or situations where individual terms in the hydrodynamical equations are not clearly dominant, analytical treatments rapidly become impractical. Additionally, some wave structures, or modes, are only supported by the full equations, being effectively "filtered out" by the simplifications. Notably, Wang & Mitchell (2014) numerically identify a Rossby-Kelvin wave mode for a planetary atmosphere that cannot be recovered in the quasi-geostrophic equations (see Gill 1980or Vallis 2006 for further details).
The detection and characterisation of a specific class of exoplanets, hot Jupiters, has provided impetus to the study of nonaxisymmetrically forced planetary atmospheres. Hot Jupiters are Jovian planets in short-period orbits close to their host star, and likely have synchronised rotational and orbital speeds such that a single hemisphere or day side faces the host star at all times (see Baraffe et al. 2010). The slow rotation (periods of ∼4 days) and Jovian radii suggest such atmospheres exist in a regime where the Rossby number is of order unity R o = U L f , where U is the characteristic flow velocity, L a characteristic length scale, and f the Coriolis parameter, meaning rotation is neither dominant nor negligible. Observational evidence has indicated the presence of fast zonal "jets" (zonally coherent flows) of a few km s −1 (Louden & Wheatley 2015). The mechanism for accelerating these zonal flows has been explored by Showman & Polvani (2011), building on the linear studies of Matsuno (1966) and Gill (1980) and using a two-layer equivalent depth approach.
To further study the linear waves and instabilities present in these atmospheres, we developed a public code, ECLIPS3D 1 (Eigenvectors, Circulation and Linear Instabilities for Planetary Science in 3 Dimensions), which we benchmark in this work. More globally, this code can be used in the study of linear stable or unstable modes within any planetary atmosphere from an arbitrary initial steady state.
We expand upon Thuburn et al. (2002), who studied propagating wave modes in an axisymmetric atmosphere at rest to include linear modes in an atmosphere with a steady background flow in three-dimensional spherical coordinates. We detail the structure of ECLIPS3D including the equations solved and the process of obtaining a solution. We detail the different sets of equations implemented (axisymmetric 2D, 3D, two-layer equivalent depth) as well as the time-dependent or independent solutions (waves, instabilities, and standing circulation). Finally, we implemented a posteriori calculation of an energy equation for a given solution. These semi-analytical results allow the verification of the frequency of the modes (and the growth and damping rate for instabilities), as well as isolation of the dominant mechanism providing insight into the physical phenomena driving the instability.
In Sect. 2, we outline the equations implemented in ECLIPS3D (with the full equations detailed in Appendix A), and the procedure for solving them, alongside the method of calculating the energy equation. We then benchmark ECLIPS3D against a range of classical calculations of waves, instabilities, and circulations in Sect. 3, including a setup similar to Showman & Polvani (2011). Finally, we draw conclusions and comment on future developments and applications for ECLIPS3D in Sect. 4.

The algorithm
2.1. Linearised equations Thuburn et al. (2002) show that even the simplest atmospheric waves exhibit behaviour that cannot be accurately expressed by separating variables (requiring a height-dependent shift in latitude; see Thuburn et al. 2002, for details). Therefore, in the general case, no assumption can be made on the mathematical expression of the wave regarding spatial coordinates. As our steady state is arbitrary, we linearise the full equations with no simplification. However, to more easily describe the main capabilities of ECLIPS3D we detail how ECLIPS3D solves the Euler equations, omitting diffusion or viscosity, although dissipative processes have been implemented (and discussed in the steady circulation case in Sect. 3.5 and Appendix A.5). The basic equation set is as follows.
where u, v, and w are the components of velocity in the longitudinal (λ), latitudinal (φ), and vertical (r) directions, ρ is the density, p the pressure, T the temperature, θ the potential temperature, p 0 a reference pressure, R the gas constant divided by mean molecular weight, c p the heat capacity, g the gravitational acceleration (and is a function of r, see Appendix A.1), Ω the rotation rate of the planet, r the radial distance from the centre of the planet, λ the longitude, φ the latitude, and finally Q is the heating rate (if present). Equations (1a)-(1c) represent momentum conservation, (1d) mass conservation, (1e) conservation of energy, (1f) is the equation of state (here an ideal gas), and (1g) defines potential temperature, closing the set. Thuburn et al. (2002) showed that potential temperature is more appropriate than normal temperature for studies of the linear modes. Finally, D/Dt is the Lagrangian or material derivative and t is time. Solving for waves or instabilities then requires these equations to be linearised. We follow the definitions of Thuburn et al. (2002) for the perturbed variables (this scaling comes from Daley 1988), which greatly simplify the equations when the steady state is axisymmetric and at rest. This choice was made for easier comparison and benchmarking with Thuburn et al. (2002), but a user of the code can change the implemented equations easily without affecting the method of solution. Namely, we write where a prime denotes a linearised variable and an i subscript the initial steady state. If the heating rate Q is non-zero, it has to be properly included in the linearised equations. When solving for waves and instabilities, we simply linearise Q and include it in the left-hand side of the equations. This is detailed in Appendix A.4. When looking for steady, linear circulation (and not free or forced waves) Q(r, φ, λ) is specified and considered small enough to only trigger a linear response. A dissipative mechanism must also be added in order to reach a linear steady state. This setup is similar to that of Showman & Polvani (2011), which is one of our benchmark cases, and is detailed further in Appendix A.5. Thuburn et al. (2002) considered the linearised equations for the case where the initial atmospheric state is axisymmetric, in hydrostatic balance and at rest. In the more general case, linearisation of each of the terms from Eqs. (1a)-(1g) must be completed as shown in Appendix A, alongside the resulting final equation set, Eqs. (A.2)-(A.5). These final, linearised equations are then discretised and solved within ECLIPS3D as detailed in Sect. 2.4. Additionally, we implemented a twolayer model following Showman & Polvani (2011), based on their Eqs. (9) and (10) (linearised versions of which are given in Appendix A.3). Other equation sets (e.g. shallow water, anelastic, and so on) and geometries could be implemented within ECLIPS3D with relative ease if required.

Boundary conditions
For inviscid flows, there must be no normal flow at the limits of the domain in order to obtain a well posed problem with complete boundary conditions. Namely, we impose that v cos(φ) tends to zero at the poles, and w is zero at the top of the atmosphere as a no escape condition. At the inner boundary, we impose a solid boundary with w equal to zero. For hot Jupiters, this requires the modelled domain to extend to high-enough pressures for the atmosphere to reach a quiescent region not involved in the acceleration of superrotation. This inner boundary condition can be easily changed if mass flows or energy transfer with the deep atmosphere need to be considered. If the density of the upper atmosphere is too low, unphysical velocities might arise. In the physical applications of ECLIPS3D so far, we have solved this problem by reducing the extent of the atmosphere but a smoothing of the higher atmosphere could be implemented (as done for example in GCMs, see Mayne et al. 2014a,b).
With the above choice of boundary conditions and implemented equations, the code will only recover standing waves in the vertical direction. However, Wu et al. (2001) expressed the importance of vertical wave propagation in the context of Matsuno-Gill structures (Matsuno 1966;Gill 1980), relevant for understanding the climate of Earth (see notably Sarachik & Cane 2010) but also regarding the spin-up of superrotation in hot Jupiters (Showman & Polvani 2011). Additionally, a numerical way to mimic evanescent waves is to impose a damping region at the top of the atmosphere, as is done for example for GCM studies of hot Jupiters (see Mayne et al. 2014a), which prevents the wave from reflecting but allows it to propagate.
In this paper, the boundary conditions have been chosen to enable comparison with previous studies, but any user of the code can easily apply different boundary conditions. Additionally, adapting the equations to include a damping layer in ECLIPS3D poses no theoretical or numerical issue. Vertically propagating waves can therefore be recovered with ECLIPS3D.

Energy equation
Following the method of Thuburn et al. (2002), we calculated an energy equation by combining the linearised Euler equations and integrating them over the whole atmospheric volume with appropriate boundary conditions. We derive this equation in the same context as Thuburn et al. (2002), with an initially axisymmetric, hydrostatically balanced atmosphere at rest. The general case is shown in Appendix B. We assume that the linearised variables X can be expressed as X (r, φ, λ, t) = X(r, φ)e −i(σt+mλ) , where the real part of σ is the mode frequency and its imaginary part the growth rate (if non-zero), X(r, φ) ∈ C, as explained in Sect. 2.4, and m ∈ Z: where the star symbol denotes the complex conjugate, with f = 2Ω cos(φ), F = 2Ω sin(φ), a the radius of the planet, H the height of the top of the atmosphere, and E = 1/2 (|u| 2 + |v| 2 + |w| 2 )/ρ i + |θ| 2 /ρ i N 2 i + |p| 2 /(ρ i c 2 i ) the sum of the kinetic, thermobaric, and elastic energies of the perturbation, with N 2 i and c 2 i the initial buoyancy frequency and sound speed, respectively. Using integration by parts we can express: with V the volume. This is possible only if V Er 2 cos(φ) drdφdλ 0 which is the case where N 2 i > 0. For this work, we assume a stably stratified atmosphere (N 2 i > 0). As stated by Thuburn et al. (2002), Eq. (4) shows that σ can only be real in this case and no instability can grow around an atmosphere at rest with no heating and the boundary conditions we have described.
Once the variables u , v , . . . are known Eq. (4) can then be integrated numerically in ECLIPS3D and be used to verify the obtained frequency and identify the dominant physical processes (e.g. in Sect. 3.1 we show that an acoustic wave is largely dominated by the terms involving the pressure, whereas a Rossby wave is dominated by the f and F terms).
When using the code, the use of this a posteriori energy equation is therefore a powerful tool both for diagnosis of the dominant physical mechanism, as well as a validation of the numerical results. It requires an interpolation of the output variables and their derivative on a common grid, which in the current version is performed with linear interpolation, but the consistency of the results with the energy equation confirm that a more sophisticated interpolation would not change the physical interpretation provided by this energy equation (see notably Sect. 3.2).

Method of solution
The derived linearised equations of motion can be expressed as: where D is a differential linear operator. If we introduce the operator A being: it is clear that A commutes with D as the latter is only timedependent through the ∂/∂t terms. Therefore, the vector subspaces of D remain stable upon application of A. As A is diagonal the kernel of D can be decomposed on the eigenvectors e σ of A. Such eigenvectors are well known: e σ ∝ e −iσt where σ ∈ C (the −i term is simply the convention we choose; the real part of σ is therefore the frequency and the imaginary part the growth rate). Therefore, coupled with appropriate boundary conditions, we can then solve Eq. (5) by decomposing them as: withû,v,ŵ,p,θ ∈ C and remembering that the evolution of the perturbed quantity is then the real part of the above expression. The actual solution is an infinite sum over all σ but the eigenmodes, and therefore the waves, are the individual projections to a single value. It is worth noting that, a priori, σ could take continuous values. For an atmosphere initially at rest, Matsuno (1966) show that only discrete values are solutions, but on the other hand baroclinic waves exhibit a continuous range of frequencies (see e.g. Charney 1947)

Time-dependent solution
From the discussion above, we can re-write our equations as a complex eigenvalue-eigenvector problem: Equation (8) can become difficult or even impossible to solve analytically. However, Thuburn et al. (2002) discretise this equation using a staggered grid of points, turning the analytical matrix B into a finite numerical matrix. This allows spatial derivatives to be calculated using finite differences. The staggered grid was selected carefully for precision and stability by Thuburn et al. (2002), where the 2D axisymmetric version is presented. We adopt a staggered grid in ECLIPS3D, shown in Fig. 1, which resembles the one used in Thuburn et al. (2002) but adapted to 3D and with different staggering of the u and v variables at the poles to simplify the boundary condition. We ) points in our grid, with N λ , N φ , and N r being the number of points in each coordinate, meaning the matrix B will be of size (5 * N tot ) 2 as there are five inter-dependent variables. However, each variable at a given point only depends on the values of all variables over the closest points in the grid. Therefore, B is an extremely sparse matrix.
Once the matrix B is filled with discretised values from Eq. (8), at the staggered grid points we must find the eigenvectors of this matrix. ECLIPS3D uses the ScaLAPACK 2 library for parallel linear algebra (Blackford et al. 1997). To express the eigenvectors of a complex matrix we first calculate the upper Hessenberg form of the matrix, then find the Schur decomposition before identifying the eigenvectors themselves 3 . Finally, the eigenvectors are returned to their original form via multiplication with the matrix of transformation.
Example cell of the 3D staggered grid adopted in ECLIPS3D, based on Thuburn et al. (2002). i discretises the longitudinal variable λ, j the latitude φ, and k the radial variable r. Here, u and v are staggered in latitude, with v running from the south to north pole and thereby having an additional latitude point. Also, p is staggered in height with w and θ with the latter two variables running from the bottom to the top of the atmosphere resulting in an additional height point.
However, this process of solving for the eigenvectors yields an eigenvector for each row in the matrix. From this set we must select those of interest, representative of physical modes in the atmosphere in question. To identify the interesting eigenvectors we employ two methods. Firstly, for instabilities with positive exponential growth rates, we assume that the modes that will lead the dynamical instability have the highest growth rate, and only select the fastest-growing modes. For the case of no instability we first filter out modes arising from numerical errors (e.g. extreme values at the poles), and then manually select modes from the solution set. Analytical expectations then determine the modes of interest. Globally, this selection process needs to be performed thoroughly and be based on analytical expectations of the modes to look for. In the current version of the code, the selection is performed independently of the matrix calculation, and therefore allows different eigenvectors to be isolated in a single ECLIPS3D run.

Time-independent solution
For the case of a steady circulation, without time-dependence and with constant heating rate, Eq. (5) can be expressed as:  Thuburn et al. (2002) and those identified in this work using ECLIPS3D.
where C is similar to B in Eq. (5) with the inclusion of a drag term if required (see Appendix A.5). Solving this problem is much easier than the time-dependent case, as we simply need to express C on the staggered grid and invert it to obtain the unique solution to these equations.

ECLIPS3D resolution
In order to achieve the highest possible resolution, we implemented two versions of ECLIPS3D: one that solves for the whole eigenvector spectrum and another one dedicated to solving a reduced number of selected eigenvectors.
For the case of numerically solving for all potential eigenvectors, the computational expense (both in computation time and memory) increases steeply with the number of points in the matrix, and can rapidly become inhibitive, limiting the resolution. Specifically, within our computational framework, ECLIPS3D can be used to solve 2D problems (axisymmetric, two layer, or barotropic equations; see following section) in cases with up to 100 × 100 points within a day of real time. The efficient parallelisation of eigenvector calculations is still an active area of research in the computer science community, and increasing the number of processors does not significantly accelerate calculations of this type. Higher-resolution problems therefore take much longer, with our benchmarking tests suggesting the time taken scales with the number of points as roughly N 2 tot or N 3 tot . Additionally, eventually with increasing numbers of processors, the communication between the processors becomes the primary overhead or limitation in the calculation, whereas using too few processors leads to saturation of the available memory (see the ScaLAPACK documentation for details). These computational limitations are amplified in 3D, where calculations at resolutions of 25×20×20 points on 64 processors require around four days. We are currently working on adapting ECLIPS3D to employ sparse matrix libraries (as discussed, the matrix we are solving is sparse) and are following the developments in computer science research regarding eigenvector calculations.
However, by solving for a reduce set of selected eigenvectors we can commensurately increase the resolution of the ECLIPS3D setup, whilst retaining a similar computational expense to the case where the target eigenvectors are not restricted. This approach can be taken when there are some existing or prior constraints on the frequency and/or growth rate of the eigenvectors, for example in the case of an instability where one requires the fastest growing mode. Solving for ten eigenvectors with a resolution of 25×20×20 points takes less than two hours to converge on 64 processors, although these calculations are still limited by the physical memory available on the processors. As mentioned, we are working on implementing sparse matrix-solving libraries in ECLIPS3D which will overcome this limitation.
Globally, the search for particular waves can be optimised by the combination of both a complete eigenvector, and specified or restricted eigenvector setup of ECLIPS3D. First, one would calculate a whole spectrum of low-resolution eigenvectors and identify the most interesting ones, before studying them in much higher resolution with the faster version of the code. However, for the time-independent solution, hence matrix inversion, the resolution can be much higher because inverting a matrix is a well parallelised and efficient process,. For example, for 25 × 25 × 25 points the matrix inversion takes less than an hour on 16 processors.
Finally, we stress that the symmetries of the problem can allow us to restrict our study to only one hemisphere, doubling the effective resolution with the same number of points. In this paper for example, the modes we present are always symmetric about the equator, but ECLIPS3D can solve for both symmetric and antisymmetric modes at the equator.

Benchmarking
We have applied ECLIPS3D to five well-studied cases from the literature to benchmark the code. These tests are explained in this section. First we reproduce the results of Thuburn et al. (2002) for a simple, hydrostatically balanced, and zonally symmetric atmosphere at rest (Sect. 3.1). This is followed by the case of an unstable jet providing a steady initial circulation as introduced by Wang & Mitchell (2014) (Sect. 3.2). We also present results for the baroclinic instability test of Jablonowski & Williamson (2006) and Ullrich et al. (2014) (Sect. 3.3). In order to implement longitudinal variation in the steady state, we study the stability of Rossby-Haurwitz waves (Haurwitz 1940), as done in Thuburn & Li (2000). Finally, ECLIPS3D is applied to the case of a linear steady-state circulation with atmospheric drag following Komacek & Showman (2016) (Sect. 3.5).

Initial atmospheric rest state
We first apply the 2D, axisymmetric version of ECLIPS3D to solve for the eigenmodes of an initially axisymmetric, isothermal, and hydrostatically balanced atmosphere at rest following Thuburn et al. (2002). Namely, the atmosphere is 80 km high with the bottom boundary at the Earth radius a = 6371 km, and the temperature is T = 250 K corresponding to N 2 = 3.83 × 10 −4 s −2 . The value of the other parameters are R = 287.05 J kg −1 K −1 , c p = 1005.0 J kg −1 K −1 , Ω = 7.292 × 10 −5 s −1 , and g = 9.8062 m s −2 . This version of ECLIPS3D needs to assume an integer wavenumber m in longitude, as in Thuburn et al. (2002). Table 1 shows the frequencies of the modes from both this study and that of Thuburn et al. (2002) revealing agreement better than 3%; the discrepancies are due to slightly different initialisations and grid staggering. When matching their setup exactly we return matching results to within machine precision. Additionally, our resulting eigenfunctions have the same shape in height and latitude and global values as those found in Thuburn et al. (2002). For example, we isolate and present both an acoustic and Rossby wave recovered by ECLIPS3D in Figs. 2 and 3, to be compared to Figs. 1 and 2 of Thuburn et al. (2002). The acoustic mode shows a vertical compression mode, with little energy in the horizontal velocities and an opposite phase between the pressure and vertical velocity perturbations. The tilt in the zonal velocities close to the pole confirms the impossibility of obtaining solutions with separate functions in the latitudinal and vertical directions in spherical geometry, as noted by Thuburn et al. (2002). The Rossby modes on the other hand are dominated by pressure and horizontal velocity perturbation, and if projected on a latitude-longitude plane we would recover the rotating winds around pressure maxima or minima in the mid-latitudes. Once again the shape of the waves clearly confirms the impossibility to separate variables.
The difference in the forcing mechanisms between acoustic and Rossby waves, namely the pressure gradient and Coriolis force, respectively, mean we expect their global features to differ. The calculations from our energy equation, shown in Table 1, are in excellent agreement with the obtained numerical frequencies. These calculations also allow one to investigate the restoring force. For example, for the Rossby waves the f and F f terms of Eq. (4) account for 90% of the value of σ, whereas they are negligible compared to the terms involving the pressure, and its derivative, for the acoustic waves.
Although this problem is axisymmetric, it can be used to test the 3D version of ECLIPS3D. In Fig. 4 we present a single mode from the 3D case, to be compared to Fig. 2. For this mode (and the other modes not presented explicitly here) we obtain the same height and latitude behaviour. For the additional dimension, longitude, we recover oscillatory modes with an arbitrary integer number, 2m, of zeros in longitude (corresponding to a wavenumber m). The frequency of the obtained modes again, as with the 2D version of ECLIPS3D, match those of Thuburn et al. (2002) to better than ∼3% (in this case errors are also introduced by the discretisation in longitude).

Unstable jet
The next benchmark case is an initial state which includes an initial velocity field. Here we follow Wang & Mitchell (2014) who identified an exponentially growing linear mode, bringing eastward momentum to the equator under axisymmetric forcing. This study essentially identifies unstable modes in an atmosphere similar to Thuburn et al. (2002) but including a mid-latitude unstable jet. Namely, the initial velocity is controlled by a given latitude φ 0 through: . Pressure (colour scale) and wind (vector arrows) for the most unstable mode obtained with ECLIPS3D from the setup of Wang & Mitchell (2014), which is to be compared with their Fig. 1a. where α controls the decay of the velocity field towards the pole. The value of α is not given in Wang & Mitchell (2014); here we chose α = 50 which mimics the shape of their initial velocity field in their Fig. 1. Wang & Mitchell (2014) identify two instabilities, firstly a well-known baroclinic instability (such as studied in Sect. 3.3), and secondly a new instability not captured by analytical treatments under the β-plane approximation. This new mode results in the convergence of eastward momentum at the equator, and is related to the Rossby and Froude numbers. Wang & Mitchell (2014) term this new instability the Rossby-Kelvin instability as it emerges from interaction between the mid-latitude Rossby waves and the Kelvin wave (an equatorially confined gravity wave with zero meridional velocity). In Fig. 5 we show the characteristics of the mode identified by our own study using the 2D axisymmetric ECLIPS3D, which is to be compared with Fig. 1a of Wang & Mitchell (2014). Figure 5 demonstrates the excellent agreement of the structure of the mode found using both ECLIPS3D and that of Wang & Mitchell (2014).
Following Wang & Mitchell (2014) we explore the effect on the most unstable mode of varying the planetary parameters. Figure 6 shows results for different φ 0 , a characteristic latitude of the initial flow (see Wang & Mitchell 2014, for definitions), and the Burger number Bu = ((N i H)/(2Ωa)) 2 where H is a characteristic height, revealing a change in the growth rate as H is altered (see Wang & Mitchell 2014, for more details). Figures 6a and b therefore represent the most unstable mode for a broad inital jet in latitude, up to 50 • , the growth rate of which is about Ω/5 and the characteristic height a third of the total height. On the other hand, Figs. 6c and d represent the most unstable mode for a narrower inital jet, with φ 0 = 35 • , the growth rate of which is about Ω/20 and characteristic height a fifth of the total height. The obtained growth rates are consistent with those presented in Fig. 2a of Wang & Mitchell (2014).

Baroclinic instability
Jablonowski & Williamson (2006) detail a baroclinic instability test for GCMs using pressure as a vertical coordinate. Ullrich et al. (2014) adapted this test for height-based GCMs. In this test a perturbation to a steady longitudinal wind at mid-latitudes leads to a dynamical instability growing in a few days (for Earth-like conditions). In full 3D GCM simulations many phenomena act simultaneously, meaning that reproducing an instability with the exact same behaviour and time evolution is unlikely. However, we can expect to reproduce the most unstable modes, which will drive the evolution of the atmosphere in the simulations.
We implemented the initial state prescribed in Ullrich et al. (2014) in both the axisymmetric 2D and 3D versions of ECLIPS3D (without the prescribed perturbation of Ullrich et al. 2014, as ECLIPS3D intrinsically perturbs steady states). We only show the 2D results as, similarly to the first test case, 2D and 3D are in excellent agreement. In this test, ECLIPS3D identifies the stable modes of Thuburn et al. (2002), slightly modified by the mean flow and the angular dependency of pressure and temperature. For the unstable modes, a continuum in frequency is returned (discretised by the numerical precision of the algorithm, controlled by the number of points in the matrix B), as expected from analytical treatment. ECLIPS3D identifies the most unstable mode at m = 5 with a growth rate of 6.4 × 10 −6 s −1 (∼2 days) which is presented in Fig. 7. Figures 4 and 5 of Ullrich et al. (2014) demonstrate that the instability dominates the flow after 8 days, consistent with our growth rate. Additionally, the shape of our instability has similar features to the thermodynamic state of the atmosphere in Ullrich et al. (2014) after 8 days. Our instability indeed exhibits a tilt in the plot of height versus latitude (Fig. 7b) as can be seen in Fig. 6 of Ullrich et al. (2014). The pressure also exhibits a similar sharp decrease just above the surface. One must note here that our results include an uncontrolled phase in longitude coming from the axisymmetry of our setup.
The only difference between 2D and 3D is the precision due to discretisation. In 3D, ECLIPS3D is limited as the size of the matrix to invert is much bigger than in 2D. Due to the shape of the staggered grid (Fig. 1), the first point in pressure is not the surface pressure and the sharp decrease is less obvious than in 2D (not shown).
Comparison between the semi-analytical calculation from our a posteriori energy equation and the numerical eigenvalue obtained for this mode reveals close agreement, on the real and imaginary part of σ, within a few percent. As explained in Appendix B, we can decompose this energy equation into three components (five in the case with meridional and vertical velocities) comprising the terms coming from the equations at rest with no angular dependency on the thermodynamic variables, the terms arising from the angular dependency of the steady state, or the terms coming from the initial zonal velocity. From analytical considerations, we expect the frequency to be dominated by the velocity terms, where the global mean flow excites the modes at specific phase velocity. However, the growth rate should be controlled by the angular-dependent terms as the baroclinic instability arises from horizontal gradients in the pressure and temperature (see e.g. Vallis 2006).
Our energy-based calculation gives a frequency of 1.09 × 10 −5 s −1 where ECLIPS3D finds 9.93 × 10 −6 s −1 . In the calculation, the velocity terms account for more than 80% of the frequency, confirming the analytical predictions. The calculated growth rate is 6.63 × 10 −6 s −1 , close to the ECLIPS3D value of 6.39 × 10 −6 s −1 , with the angular terms accounting for 96%. These results show that our a posteriori energy equation can be a powerful tool to obtain insight into the physics of numerical eigenvectors.

Rossby-Haurwitz waves
The previous steady states we studied are axisymmetric and include no background velocity (Sect. 3.1) or only a zonal A36, page 8 of 15 F. Debras et al.: ECLIPS3D   Fig. 6. Pressure perturbations (colour scale) and wind vectors (arrows, panels b and d only) for the case of a planet with the same radius and rotation rate as the Earth but an isothermal temperature pressure profile set at 500 K. Panels a and c: latitude against height at a longitude of 200 • ; panels b and d: longitude against latitude at a height of 25 000 m. We report the values of φ 0 and Bu as defined in Wang & Mitchell (2014) and the growth rate σ growth : (a) and (b): φ 0 = 50 • , Bu ∼ 0.2 and σ growth 2Ω = 0.12, and (c) and (d): φ 0 = 35 • , Bu ∼ 0.05 and σ growth 2Ω = 0.026. These results indicate that the growth rate of the most unstable mode is dependent on the characteristic height of the wave as found in Wang & Mitchell (2014).  .2 and 3.3). Ideally, we would also benchmark ECLIPS3D using a test including meridional velocities as well as a dependency on longitude. Unfortunately, there are no such non-linear steady states in 3D, in which the analytical theory can provide us with predictions to compare with. We therefore consider a 2D non-axisymmetric problem, with steady zonal and meridional winds: Rossby-Haurwitz waves. Here we identify the most unstable modes around two steady configurations of this setup which we detail below. Rossby-Haurwitz waves are analytical solutions of the nonlinear barotropic vorticity equations. They were discovered by Haurwitz (1940) by perturbing the non-divergent equations and solving them non-linearly. If the flow remains incompressible at all times, these waves remain analytical solutions of the full equations and propagate without changing their form at constant speed. With an appropriate choice of parameters, this speed can be zero, and these waves become a stationary, steady solution of the non-divergent barotropic vorticity equations, hence another test case for ECLIPS3D.
The non-divergent barotropic equations are simply Eqs. (1a) and (1b) with w = 0 with an imposed null divergence: This latter equation does not involve any time derivative, and therefore we have to slightly adapt the structure of the code for this setup. Instead of solving an eigenvector problem, we solve a generalised eigenvector problem where where x is an eigenvector, B the linearised matrix of equations, and C a diagonal matrix with some zeros in the diagonal. The linearisation is straightforward, as the equations are similar to the full set of equations and the divergence equation is already linear.
The interest of this test lies in the stability analysis of these waves. Hoskins (1973) showed that Rossby-Haurwitz waves are stable for longitudinal wave number R < 5 and unstable for higher wave numbers. However, in the analysis of Hoskins some triad interactions were missing, as simplifications were required to treat the problem analytically. Inspired by Baines (1976), Thuburn & Li (2000) resolved the issue by showing that a Rossby-Haurwitz wave of wavenumber 4 is unstable, because of an interaction with wave numbers 1, 3, and 5 (see also Lynch 2009).
ECLIPS3D does not make any assumption on the shape of the perturbation needed to trigger an instability, nor on the instability itself. We therefore expect to obtain unstable modes around a steady R = 4 wave. Our setup is similar to the classical benchmarking test of Williamson et al. (1992) and Thuburn & Li (2000), with a vorticity ψ being: where K and ω are constants and R is the longitudinal wavenumber. In order to obtain a steady, stationary wave we must also impose (see Haurwitz 1940): With R = 4 and Ω = 7.29 × 10 −5 s −1 , the Earth rotation rate, this leads to ω ≈ 5 × 10 −6 , close to the value chosen by Williamson et al. (1992) and Thuburn & Li (2000) ω ≈ 7.8 × 10 −6 . With R = 2, ω ≈ 1.5 × 10 −5 . In accordance to their setup, we impose arbitrarily K = ω. We present the initial steady states for R = 4 and R = 2 in Figs. 8a and c. As expected from Thuburn & Li (2000), we obtain a linear instability for the R = 4 setup, displayed in Fig. 8b. This instability oscillates with a period of ≈3 days, with an exponential growth timescale of 6 days. The timescale for instability is globally coherent with the results of Thuburn & Li (2000), which find that the flow becomes significantly altered after day 20. For R = 2, we also find an instability with striking resemblance to the R = 4 instability in shape, a period of just under a day, and growth timescale of 4 days. This mode is shown in Fig. 8d, and this result is in contradiction with Hoskins (1973) but in accordance with Baines (1976) who shows analytically that all the R ≥ 2 (n ≥ 3 in their study, where we have R = n − 1 here) Rossby-Haurwitz waves can be unstable. Interestingly, we also found an instability for a R = 3 Haurwitz wave but with a growth timescale of more than a hundred days. Such an instability would probably be smoothed out by any source of dissipation or diffusion in a GCM. Globally, our code agrees well with the theoretical study of Rossby-Haurwitz waves, and demonstrates the proper treatment of longitudinal-dependent steady state and meridional velocities in ECLIPS3D.

Linear steady circulation with drag
As discussed in the introduction, our development of ECLIPS3D was largely driven by studies of the acceleration of zonal flows in hot Jupiter atmospheres. For these planets, analytical studies have shown linear steady states to be of vital importance (see Showman & Polvani 2011). Therefore, we have also implemented the capability to calculate linear steady states in ECLIPS3D (see Sect. 2), which is a much simpler process compared to the identification of eigen modes. To benchmark this section of the code, we compare our results to those obtained in the study of Komacek & Showman (2016), in particular the case they present in their Fig. 5 where they solve the full Navier-Stokes equations, but with such a low heating rate that only the linear terms contribute. This requires the addition of a linear drag in the linearised equations following the depthdependent behaviour of that adopted by Komacek & Showman (2016). Additionally, the heating is performed via a Newtonian relaxation with a height-dependent radiative constant. The resulting equation set is shown in Appendix A.5. Figure 9 then presents the resulting linear circulations obtained using ECLIPS3D which show excellent qualitative agreement with the results of Komacek & Showman (2016;see their Fig. 5). Komacek & Showman (2016) do not present the vertical structure of their circulation.

Conclusion
In this paper, we introduce and benchmark ECLIPS3D: a parallel code for identifying linear instabilities, waves, and circulations around a steady state of the Navier-Stokes equations in planetary atmospheres. The linearised equations only omit viscosity and an a posteriori energy equation is used to identify contributions from each component. The time-dependent eigenvector solution or time-independent matrix inversion calculations are performed through discretisation onto a staggered grid and subsequently ScaLAPACK routines.
The benchmarks cover various well-studied wave and instability tests, namely a simple atmosphere at rest (Thuburn et al. 2002), a Rossby-Kelvin unstable jet (Wang & Mitchell 2014), Fig. 8. Pressure (colour scale, a and c in pascals and arbitrary units for b and d) and horizontal winds (arrows) for Rossby-Haurwitz waves and most unstable modes. Panels a and c: initial steady states for R = 4 and R = 2, respectively (see text). Panels b and d are the most unstable mode obtained with ECLIPS3D for the R = 4 and R = 2 setup, respectively. a baroclinically unstable jet (Ullrich et al. 2014), an unstable Rossby-Haurwitz wave (Thuburn & Li 2000), and a linear circulation with atmospheric drag (Komacek & Showman 2016). For all these setups, ECLIPS3D is able to produce excellent qualitative agreement with the previous works. We demonstrate that our a posteriori energy equation is a viable tool to verify the results and identify the dominant terms. We are currently preparing a follow-up study to explore the momentum transfer in hot Jupiter atmospheres and explore the stability of the initial conditions for GCMs using ECLIPS3D (Debras et al., in prep.). ECLIPS3D currently has several limitations, primarily its computational efficiency, leading to limitations on resolution, particularly for 3D cases. We are working on several methods to improve this issue for example using libraries adapted to sparse matrices, or splitting the eigenvector solution into several submatrices as opposed to a single large matrix (potentially useful as the time taken to solve this type of problem increases faster than linearly with matrix size). This splitting of the matrix may be particularly well suited to a spectral decomposition as we are searching for the most unstable mode, and are not necessarily trying to capture the entire "shape" of the mode. This could be done through spherical harmonics in the horizontal or Chebyshev's spectral decomposition in the radial direction.
Despite its limitations, ECLIPS3D in its current version still represents a powerful resource which can be used to study instabilities in 2D situations under axisymmetry or cases where a two-layer model is applicable, or for low-resolution 3D problems. The code itself can easily be adapted to different situations, in spherical coordinates, with additional physics or alternative boundary conditions. As the structure of the code is independent of the underlying equations, meaning alternatives can easily be implemented in terms of symmetries and coordinate systems.
Finally, ECLIPS3D could be applied to a wide range of astrophysical problems. The most obvious one, for which ECLIPS3D was designed, is the study of instabilities and linear circulations for planetary atmospheres, but the range of applicability is greater. Asteroseismology for example requires the need to linearise the equations of motion and identify the leading modes, sometimes with complicated circulation or thermodynamic state inside the star. Adapted to cylindrical geometry, ECLIPS3D could be a powerful tool to identify the possible instabilities in protoplanetary disks, where instabilities creating pressure traps are proposed to be an efficient way of making planets through core accretion. The addition of a magnetic field in the equations implemented in the code would not pose any theoretical A&A 631, A36 (2019) Fig. 9. Pressure (colour scale, Pascals) and horizontal winds (arrows) as a function of latitude and longitude at a height of 5 × 10 6 m, for a steady-state circulation with a forcing of ∆T = 100 K (see Komacek & Showman 2016, for definition). Results to be compared to Fig. 5 of Komacek & Showman (2016). (a) τ drag = 10 5 s and τ rad = 10 4 s. (b) τ drag = 10 3 s and τ rad = 10 3 s. (c) τ drag = 10 6 s and τ rad = 10 4 s. (d) τ drag = 10 4 s and τ rad = 10 6 s. challenge either, which could provide a significant amount of information on the linear behaviour of astrophysical fluids in more general cases.

Appendix A: Equations in ECLIPS3D
A.1. General 3D case In this appendix we detail the derivation of the full linearised equations for the longitudinal component of the momentum equation, and in the interests of brevity provide only the final expressions for the remaining components. These equations assume a dependence of g on r as g ∝ 1/r 2 , and every other quantity is dependent on the three spatial variables, longitude, latitude, and radial distance from the centre of the planet. Here, we consider that Q is equal to zero for simplicity, and relax this assumption in Appendix A.4.
The longitudinal equation of momentum is as follows.
Therefore, the final five perturbed equations are with ∂X /∂t = iσX .

A.2. Two-dimensional axisymmetric case
For the axisymmetric case, the equations are directly obtained from the 3D case by choosing a longitudinal wavenumber m such that X (t, r, φ, λ) = X (r, φ)e iσt e i(m/2π)λ .

A.3. Two-layer equivalent depth
The reference for this particular case can be found in Showman & Polvani (2011) or Vallis (2006). We consider a dynamic layer above a quiescent layer, reservoir of mass or energy, and study the horizontal winds u and v as well as