{'fraction_non_alphanumeric': 0.03915263321203915, 'fraction_numerical': 0.011903754478011903, 'mean_word_length': 4.201291079812207, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Finite system size effects on Drift Wave stability

Unstable electrostatic resistive modes, driven by density gradients, are identified in a bounded sheared slab. The boundary conditions play a crucial role and are shown to change the nature of the problem, which is related to so called"universal"mode. The dispersion relation and the structure of the eigenmodes of the instability are derived and are shown to depend on a limited set of dimensionless parameters. The occurrence and possible impact of these modes on numerical simulations and actual plasmas are discussed.


I. INTRODUCTION
Density gradients are a necessary feature of confined plasmas and are also common in natural systems. The inhomogeneity is a source of free energy, which can amplify naturally occurring plasma oscillations and lead to turbulence and transport of particles, energy and momentum. It is therefore no surprise that density driven instabilities, i.e. unstable driftwaves, were the subject of intense theoretical research in plasma physics since its very beginning.
The controversial story of the so called "universal" instability started with the pioneering work by Krall and Rosenbluth [1], which first identified it in a slab with a shearless magnetic field and in an electrostatic collisionless regime (the drive was electron wave-particle resonances coupled with the density gradient). The same authors introduced in a later paper [2] the effect of the magnetic shear, which turned out to be a very efficient mechanism to stabilize the mode in experimental conditions. In the following years, it was suggested that the instability could still manifest in a sheared geometry as a transient effect [3], i.e. not as a normal eigenmode. These so called convective modes were described as wave packets which initially grow in amplitude when they are in the region of local instability (close to the resonant surface), but eventually decay as they move in the stable region, coherently with the overall stability of the system.
Coming back to the normal eigenmodes, the mathematical framework of the problem was significantly complicated by the magnetic shear [2], since the derivation of the dispersion relation now required the solution of a differential equation (the magnetic shear prevents the use of Fourier transforms in the cross filed direction). This implies that proper boundary conditions must be imposed in order to find the correct result. In [2], and in all the following collisionless works, the ion Landau damping was not explicitly treated but it was assumed to stabilize the mode far away from its resonant surface. Indeed, the Landau damping is more effective for modes with a large parallel wave number, which is proportional to the cross field distance from the rational surface. The boundary conditions for the "universal" instability must therefore allow a smooth connection between the region where the Landau damping is neglected and the region where it is important. It is with this argument that Perlstein and Berk [4] challenged the results in [2], where the eigenmodes were assumed to decay to zero sufficiently far away from the resonant surface (i.e. at the entrance of the Landau damping region). The authors of [4] sustained that the correct boundary conditions should allow for "outgoing waves" which would transport their energy outwards where it would eventually be dissipated by wave-particle interactions. This simple correction led to a new dispersion relation predicting unstable normal eigenmodes (more unstable than the convective modes), thus resurrecting the "universal" instability.
Ten years later, two independent groups [5,6] came to the conclusion that the electron dynamics was not properly treated in [4] (the principal part of the plasma dispersion function was not included, an inheritance from [2]). With "outgoing waves" and the complete electron response, under no condition normal eigenmodes could be found in the sheared slab. These results led to the apparent conclusion that the "universal" instability was not so universal after all.
With the realization that the collisionless "universal" mode is stable in sheared magnetic fields, the only piece of physics missing was the effect of the collisions, which was studied in [7,8]. This was done on the basis that, in drift wave dynamics, dissipative effects in the collisional regime play a destabilizing role similar to resonant electron-wave interactions in the collisionless regime [9]. In these papers the approach is completely fluid (no wave-particle interactions), an "outgoing wave" boundary condition is used and the problem was solved with both perturbative techniques (separating the resistive layer from the outer solution) and WKB methods. The surprising result was that the resistivity increases even further the stability of the system. However, other effects such as toroidicity [8,10,11] or coupling to trapped electrons [12] can overcome the magnetic shear stabilization and produce linearly unstable modes. In addition, drift waves can be metastable and sustain turbulence through non linear coupling [13,14].
The work that we present is an ideal continuation of this line of research. In particular, we describe how wave reflection in a finite size system can destabilize the resistive "universal" instability. Differently from the standard calculations, the boundary conditions that we use could be interpreted as an overlapping of "outgoing" and "ingoing" waves which create an unstable standing wave in the plasma. Our analysis is performed in the framework of a linearised reduced two fluid model in which finite resistivity is included and wave particle interactions are neglected. Both our model and some of the analytical techniques that we use have similarities with those employed in [7,8], although the final conclusions are opposite (because of different boundary conditions).
The most straightforward application of our results is to numerical simulations, which are intrinsically bounded systems. Also periodicity in the cross field direction, a common boundary condition in codes, would entrap energy and produce the unstable standing wave that we describe below. In addition, our boundary conditions might also be compatible with the energy exchange that occurs between coupled resonant surfaces in toroidal geometry [10,11], so that the results of our simple slab analysis could have counterparts in realistic plasmas.

II. MODEL AND EQUATIONS
We investigate a reduced electrostatic two-fluid model, with isothermal electrons, cold ions and valid for small but finite β. Our equations were originally derived by Hazeltine et al. [15] and are consistent with those employed in previous studies of the "universal" instability [5,7,8]. The dimensionless equations evolve the plasma density, n, the plasma vorticity, U, which is related to the electrostatic potential, φ, and the parallel ion velocity, v. they are: Here J is the parallel current density (J eq is assumed to be constant) and Note that [f, ψ eq ] ≡ ∇ represents the parallel gradient since ψ is the poloidal magnetic flux. The model contains three dimensionless parameters, which are the normalized resistivity, η = 1/S, where S is the Lundquist number, the plasma β ≡ 4πp/B 2 (note the unusual coefficient 4 instead of 8) and ρ ≡ ρ s /L s , which is the ion Larmor radius calculated with the electron temperature (which we call ρ s ) normalized to an equilibrium length scale (here the magnetic shear length, L s ). The problem is solved in a slab domain which extends in the radial direction, x, along which the equilibrium density changes, and a transverse ("poloidal") direction y, which is assumed to be periodic. Both coordinates are normalized with respect to L s . The time is normalized with respect to the Alfven time, with v A ≡ B/ √ 4πn 0 m i (n 0 is a characteristic density and m i is the ion mass). More details on the normalizations and on the physics described by the different terms can be found in [17,18]. Note that our model is an extension of the well known Hasegawa-Wakatani equations [16], that are retrieved in the limit β = 0. In other words, with respect to [16], we treat the parallel compressibility in the density equation in a more complete way, allowing the coupling with the parallel ion velocity.
Assuming that the perturbations can be expressed in the form, φ(x, y, t) = φ(x)exp(γt + ik y y), and linearizing the equations we obtain: where: and the prime symbol represents derivation with respect to x. In the last equation, we have introduced the electron diamagnetic frequency (normalized to the Alfven frequency, , ω * ≡ (ρ s /L n )k y c s τ A , with c s ≡ T e /m i . This quantity describes the equilibrium density gradient (assumed to be locally constant) through the length scale L n ≡ L s n eq /n ′ eq (remember that ′ represents derivation with respect to the normalized radial coordinate).
The normalized diamagnetic frequency is another dimensionless parameter, which controls the dynamics of the system together with the other three introduced above. Note also that all the quantities with a hat symbol are further normalized with respect to k y . This notation will be dropped in the following but η, γ and ω * should always be assumed thus normalized unless otherwise stated.
The magnetic field is assumed to be sheared so that the equilibrium component of B y depends on x and it vanishes at x = 0 [B y (x) = ψ ′ eq ≈ x around the resonant surface, using a Taylor expansion], while B z is a constant. This also implies that the parallel wave number of the perturbation is given by: k = k y x since ik = B eq · ∇ and z is taken to be an ignorable direction.
The problem 6-7 is completely defined only when the boundary conditions are imposed, which, as discussed in the introduction, can strongly affect the final result. While the "outgoing wave" boundary conditions would require a solution that carries energy outwards, we impose that the wave is completely reflected at a finite distance L x from either side of the resonant surface, which is equivalent to setting φ(±L x ) = 0. Therefore, the energy does not leak out of the system, not differently from a configurations with periodic boundary conditions in x (which would lead to very similar conclusions as discussed in the Appendix).
We solved Eqs.6-7 in various limits in order to clarify the effect of the different terms.
In some cases, it is possible to obtain an exact solution of the problem, while in others, the solution can only be found by making use of small parameters and applying perturbative techniques.
For convenience, and similarly to [7,8], the potential Q is rewritten in the following form: with:

A. Dispersion relation in local approximation
To introduce the problem, we first study the dispersion relation in the limit of zero magnetic shear. This can be done by assuming a constant parallel wavelength, k = k ,0 .
This approximation removes all the radial dependencies from the equilibrium coefficients of the equations, thus allowing a Fourier transformation also in the x direction. The problem becomes scalar and trivial to solve once x is replaced by k ,0 /k y and φ ′′ by −k 2 x φ in Eqs.6 and 7. The dispersion relation is: where k 2 ⊥ ≡ k 2 x + k 2 y . In order to find a simple solution of this equation, we assume η ≪ 1, and we expand in this small parameter, so that γ = γ 0 + γ 1 , with γ 1 of order η. This gives: where the solution with the plus sign in Eq.15 describes a resistive instability driven by the density gradient. Note that in the case β = η = 0 we recover the standard drift waves, while for ω * = η = 0 we find sound waves. Equation 13 has a third solution, which is singular for η → 0, but it is stable and therefore not interesting. The growth rate in Eq.15 gives some insight on the nature of the instabilities we discuss in this paper. They are resistive (for η = 0, γ 1 = 0) and driven by the density gradients (through ω * ). The same features are observed also in the presence of magnetic shear, as the following Sections show.

III. EXACT RESULTS
A. Case with: ρ = β = 0 At low temperature ρ and β are small and they can be neglected in Eq.7 so that: where δ ≡ δ − Λ/x 2 R and ǫ 2 ≡ −Λ/x 4 R are evaluated at β = 0. In this limit, Eq.6 coupled with Eq.16 has an exact solution in terms of Whittaker functions [19]: where κ = 1 4 δ ǫ , C 1 and C 2 are complex integration constants. Modes with even and odd parity with respect to x are characterized by φ ′ (0) = 0 and φ(0) = 0, respectively. The first condition corresponds to: while the second gives: Imposing also the reflecting boundary condition, we obtain the dispersion relation M κ,±1/4 ( ǫL 2 x ) = 0, where the plus and minus signs are for odd and even modes, respectively. In the limit of small resistivity, ǫ ≫ δ, this simplifies to: J ±1/4 (−i ǫL 2 x /2) = 0. We use now approximate relations to find the zeros of the Bessel function [19], obtaining the dispersion relation: x , with n ∈ N and the usual parity convention for the plus and minus sign. Using the definition of ǫ, we find: where we have used the fact that ℜ( ǫ) = 0 implies that ℑ( ǫ 2 ) = 0 and therefore ℑ( ǫ) 2 = −ℜ( ǫ 2 ). In the previous expressions, γ r ≡ ℜ(γ) and γ i ≡ ℑ(γ) are the growth rate of the mode and its rotation frequency.
Equation 20 shows that, in the spectrum, the unstable eigenvalues are located on a curve that does not depend on the resistivity, but only on the diamagnetic frequency. Consequently, also the growth rate of the fastest growing mode depends only on ω * : γ r,max ≈ 0.3ω * .
Note also that both even and odd modes lie on the same curve, described by Eq.20. For η = 0 Eqs.8-12 yield a quadratic potential Q(x, γ) that is structurally similar to the one discussed in the previous section (i.e Eq.16). Hence, the procedure described above can be straightforwardly extended to the ideal case (η = 0), valid for high temperature collisionless plasmas. In this regime, the condition ℑ(ǫ 2 ) = 0 translates into γ i γ r = 0 so that the mode either rotates, but does not grow, or vice versa. The condition on ℜ(ǫ 2 ), analogous to Eq.21, would give: As the right-hand side of the previous equation is negative defined, the only acceptable solution is |γ r | < |γ i |, which leads to γ r = 0 and implies no growing modes. This result implies that in the fluid limit the resistivity is a key ingredient to destabilize the modes discussed in this paper. A proper investigation of the effect of finite system size in collisionless regimes would require a kinetic treatment which is beyond the scope of the present study.

IV. PERTURBATIVE APPROACH
As suggested by the previous Sections, a finite resistivity is essential to produce unstable modes in finite size fluid systems. This contrasts with earlier analytic results [7,8], which used outgoing wave boundary conditions and predicted stability and even a damping due to η. We now want to asses the properties of the unstable modes when all the physics is included in the problem. Unfortunately, in this case an exact solution of Eq.6-7 is not available. On the other hand, an asymptotic matching approach, similar to the one used in [7,8], provides an approximate solution.

A. Outer region
We start by studying the outer region, located at x ≫ x R . In this case, Q reduces to: the solution of which can be expressed as: (compare with Eq.17), with the following indexes: κ ≡ 1 4 δ ǫ , µ = − √ 1 + 4Λ/4, C 1 and C 2 are the constant of integration.

B. Inner Region
The inner region is characterized by x ∼ x R ≪ 1. In this limit, we have: and the solution of Eq.6 can be expressed in terms of Legendre functions: where C 3 and C 4 are constants of integration and ν ≡ 1/2( The Legendre functions are multivalued on on the real axis between −1 ≤ ℜ(x/ix R ) ≤ 1, so here they have a branch cut. This can be taken into account by taking ℑ[(ix R ) −1 ] < 0 (valid for unstable modes, ℜ(x R ) > 0) and x > 0, i.e. we go to zero from the third or fourth quadrant. At x = 0, this leads to: which implies that the odd modes [φ in (0) = 0] have: Similarly, for the even modes [φ ′ in (0) = 0]:

C. Matching and dispersion relation
Let's now match the inner and the outer solutions. To do this, we first take the small x limit of the outer solution: The large x limit of the inner solution is: Reminding that −ν = 1/2 + 2µ, we match the first term of φ in with the second of φ out and the second term of φ in with the first of φ out .
This gives: From the first equation, we find: If we assume a finite box size, we have to impose the condition: which leads to: Note that W κ,µ (ǫL 2 x )/M κ,µ (ǫL 2 x ) goes to zero when L x goes to infinity so that in this limit we correctly recover the dispersion relation in [7,8]. It is convenient to express W κ,µ in terms of M κ,±µ , which gives a more compact version of the new dispersion relation:

D. Dimensionless parameters and relevant regimes
Only a limited number of combinations of dimensionless parameters appear in the dispersion relation, Eq.38. They are: with γ = γ/(iω * ) the rescaled complex frequency. In an abstract form, we can therefore write Eq.38 as f (γ, √ ηω * /ρ 2 , ρk y , L x /ρ, ρ √ β/ω * ) = 0 where only four parameters control the problem. While some of them have a straightforward interpretation, it is useful to point out that ρ √ β/ω * = L n /L s (in the following we use only the latter definition). These new parameters determine the properties of the modes even when no perturbative expansion is involved, as can be seen by properly renormalizing Eqs.6 and 7. We note, for example, that the conclusions in Section III A are valid when L x /ρ ≪ L s /L n together with L x /ρ ≪ ( √ ηω * /ρ 2 )γ 1/2 . Similarly, the regime described in Section III B corresponds to the opposite limit: L x /ρ ≫ ( √ ηω * /ρ 2 )γ 1/2 . In addition, when L x /ρ ≫ γ 1/2 (L s /L n ) 1/2 , we are able to reproduce the results discussed in [7,8]. In the next Sections, we study other relevant regimes in which the complete dispersion relation, Eq.38, has a transparent interpretation.
From Eq.39 it is clear that in this limit also |ν| ≪ 1, which leads to several simplifications (we discuss solutions for which |γ| is of order unity or smaller) and that µ ≈ −1/4. We start by studying the even modes. In this case, Eq.38 becomes: where: h ≡ −i 3/2 (γ + 1) −1/2 (L n /L s ) 1/2 . This can also be written as: Note that while M κ,ǫ (z) refers to Whittaker functions, the notation M(a, b, z) is used for confluent hypergeometrics [19].
We discuss now the odd modes, which give the following dispersion relation: In this case, the right hand side is higher order and we need to treat explicitly the order ν term at the left hand side, which is dominant. This observation allows us to determine the correct dispersion relation directly from Eq.38, in which we can neglect the right hand side.
The solution is obtained by evaluating the zeros of M κ,−µ (ǫL 2 when ν is small. The zeros of a confluent hypergeometric are given by an approximated expression [19]: when n ∈ N is sufficiently large. Applying this formula to our problem, we find: which is independent from the parameter L n /L s . Note, however, that this solution is valid only if the order ν term on the left hand side of Eq.42 is larger than the right hand side. This is true when πν 1/2 h ≪ 1, which corresponds to π 2 √ ηω * /ρ 2 ≪ (L s /L n ) 1/2 . If we remain in this limit, Eq.44 is valid and it can be written as: with k 2 ⊥,n ≡ k 2 y + n 2 π 2 /L 2 x and n ∈ N. The maximum growth rate for the odd modes can be easily calculated from Eq.45 and, in the limit with corresponding rotation frequency ℑ(γ) max ≈ −ω * /2.
1. L n /L s → 0 limit for the even modes

Finite L n /L s corrections
Despite our best efforts, we could not identify a limit for Eqs.41 and 42 which contained L n /L s effects and that could be expressed with a transparent analytical formulation. As a consequence, this regime was studied numerically with a thorough characterisation of the complex frequency as the dimensionless parameters are varied. The results of this investigation are reported in Section V B 1.

V. NUMERICAL RESULTS
In order to verify our analytic results, we solved the linear version of Eqs.1-5 with a finite difference spectral code. Its output provided the full spectrum of the eigenvalues and eigenfunctions of the system, including sub-dominant instabilities and stable modes. In this Section we verify numerically the results obtained in Sec. III and IV and we also determine how the modes behave at finite L n /L s values, where no simple analytic limit was found.
A. small ρ and β Inspection of Eq.21 reveals that the growth rate of the modes depends only on the combi- Indeed, using the dimensionless complex frequency introduced in

Sec.IV D, we have that the unstable branch of the dispersion relation is
with a ≡ 4π 2 (ηω * /L 4 x )(n ± 1/8 − 1/4) 2 . This dispersion relation perfectly matches the numerical spectrum, given in Fig.1. In Fig.2 we give an example of the odd and even parity eigenfunctions associated with the n = 7 mode number. A similar oscillating structure of the perturbations is retrieved also in the more complicated regimes described in the following sections.

B. Small
√ ηω * /ρ 2 regime In order to compare the numerical spectrum with the dispersion relations Eqs.45 and 50, we extract from them the expressions relating ℜ(γ) to ℑ(γ): where Eq.51 describes the even modes, Eq.52 the odd modes and Θ ≡ π( √ ηω * /ρ 2 )(ρ/L x ), Fig.3 we compare the previous expressions with the numerical data for different values of √ ηω * /ρ 2 , finding excellent agreement. The small discrepancy around the maximum growth rate is due to the fact that the regimes investigated are not sufficiently asymptotic. We also tested the predictions of the dispersion relations Eqs.45 and 50 for the complex frequency of the modes (not just the shape of the spectrum), obtaining again a good match with the numerical results.

Finite L n /L s corrections
In this Subsection, we discuss the effect of a finite L n /L s on the stability of the system and we therefore complete the characterisation of the modes that was interrupted in Subsection IV E 2. In the limit of small √ ηω * /ρ 2 , we identified numerically that the growth rate can be written as: where γ 0 is the L n /L s = 0 growth rate given by Eqs.45 and 50. In other words, the effect of the resistivity and of L n /L s are independent from each other and F is the correction function to be numerically characterised. In addition, it is worth noticing that, in the fluid limit, k y ρ has only a weak effect on the correction function F . Indeed, for our equations to be valid, k y ρ must be much smaller than unity and this parameter appears in Eqs.41 and 42 only in κ as (1 + k 2 y ρ 2 ) (see Eq.39). For modes with k y ρ ∼ 1 a kinetic treatment would be needed, but this is outside the scope of the present work. Note also that F has an extra hidden parameter, the mode number n and it is different for even and odd eigenmodes.
In the numerical investigation we studied cases with L x /ρ = [5; 10; 15; 20; 25; 30; 35] and ρk y = 0.1 (simulations performed with ρ = 0.5 produced results that were very similar to those presented). We fixed √ ηω * /ρ 2 = 0.12, but we checked that our results were not changing for smaller values of this parameter (we calculated F for √ ηω * /ρ 2 as small as 0.04 without finding differences). crosses a finite critical value. In addition, surprisingly, we find windows of L n /L s in which some of the even modes can become more unstable (i.e. F > 1). The amount of this shear induced destabilization can be significant and depends on the mode number, on the size of the system (i.e. L x /ρ) and weakly on ρk y . It is also interesting to note that small n even and odd modes can invert their stabilization trend and form a second unstable branch with relatively small growth rate (especially for the even modes) but less effected by the magnetic shear (see the odd n = 1 mode in Fig.4). In general, as L n /L s is increased and the even modes approach their marginally stable state, their rotation frequency matches the diamagnetic frequency, as can be seen in Fig.5. On the other hand, the second unstable branch described above is associated to quickly rotating modes with ℑ(γ) > ω * , and has therefore a different character. The situation is different for the odd modes, which do not stabilize at a specific rotation frequency. On the other hand, also the odd modes enter the second unstable branch when their frequency is higher than |ω * |.   perturbations (i.e. small and large k x ) are more effectively damped than intermediate mode numbers. In addition, the mode numbers of the most resilient modes shift to larger values as L x /ρ increases. Note also that the windows of shear destabilization discussed above occur only for the even modes and that are more effective for L x /ρ ≈ 10 and less significant for smaller and larger values of L x /ρ. The new branch of the odd modes is clearly displayed in the low n region of Fig.7, where it becomes more prominent as L x /ρ is increased.
It is interesting to determine the behaviour of the most unstable mode as a proxy of the overall stability of the equilibrium. In Figs.8 and 9 we plot ℜ(γ) max /ω * as a function of L n /L s and L x /ρ for the even and odd modes. In the plane L n /L s = 0 we have added a dashed line representing the theoretical predictions given at the end of Secs.IV E and IV E 1, finding in both cases a good agreement. For the even modes, the maximum growth rate shows an approximately linear decay as a function of L n /L s , while for the odd modes the damping is well represented by a cubic function. To facilitate the comparison between these curves, we duplicated them on the L x /ρ = 5 plane.
The critical L n /L s for the complete stabilization of the modes as a function of the system size is well represented with a power law decay. In particular, we find that our numerical results are well matched by:  for the even modes and: for the odd modes. Note the relatively weak exponent in both expressions, suggesting that the critical value remains significant even in large systems. These two expressions are probably the most important result of the present work as they suggest that the unconditional drift wave stability obtained in [7,11] is exclusive to infinite systems in which all the energy is dispersed through the boundary conditions. Equations 54 and 55 and the numerical data they represent are plotted in the ℜ(γ) max = 0 plane of Figs.8 and 9.

VI. DISCUSSION
In all the limits discussed, a combination of density gradients, finite resistivity and finite system size can destabilize a spectrum of unstable modes. The effect of the magnetic shear is generally stabilizing, although relatively small values of it can make some of the modes more unstable (this effect never occurs for the fastest growing mode).
The destabilization is due to the fact that the reflecting boundary conditions that we are imposing trap the energy of the waves. The same effect would therefore occur in periodic configurations, which are often used in numerical simulations (see the Appendix). It is not uncommon to find in literature works that employ numerical domains of a few tens of ρ and relatively small L n /L s , compatible with a significant spectrum of unstable modes. One example is [20] in which turbulence is driven by the nonlinear version of the modes presents here. Our work suggests that reflecting or periodic boundary conditions can destabilize or enhance the growth rate of the drift waves and therefore affect the anomalous transport estimated in numerical simulations.
A real plasma, however, is an open system in which the energy is, in principle, not bounded to a specific region. Nevertheless, the destabilizing mechanism that we described could play a role in a number of realistic situations. For example, plasma inhomogeneities can reflect part of the energy back to its source (i.e. the resonant surface), similarly to what would happen to a pulse travelling in a rope consisting of two sections of different thickness.
As the reflection would not be perfect, we expect that the calculations we presented would describe the worst case scenario (i.e. the highest limit for the growth rate). In addition, in toroidal systems, different poloidal modes are connected with each other through curvature coupling. This would allow modes that resonate at different positions to exchange energy so that an outgoing wave emitted on one resonant surface might look like an ingoing wave for a neighbouring surface, as originally proposed by Taylor [10]. Using this picture, our calculations in the simpler and more intuitive slab geometry helps to clarify the mechanism which destabilizes drift waves in toroidal configurations [11]. In this case, L x would represent the distance between two resonant surfaces.
It is important to notice that, for the instability to occur, the energy reflection must take place before the wave reaches the region where the ion Landau damping is strong.
This region is characterized by the requirement that k v th,i /ω ∼ 1, where k = k y x/L s and v th,i = (T i /m i ) 1/2 is the ion thermal velocity (see, e.g. [4]). The nature of the problem is therefore determined by the length scale L Ld ≡ L s ω/(k y v th,i ) and the condition L x < L Ld is necessary for instability. Note that assuming ω ∼ ω * , we have that L Ld = ρ s (L s /L n ) T e /T i .
In the particular case of cold ions treated in this paper L Ld → ∞, thus assuring that our calculation is consistent.
In our model, the role played by the resistivity is non trivial. Indeed, it is essential in order to drive the unstable modes and appears in the equations in the form of a singular perturbation. In this respect, the unstable drift waves we discussed are similar to the resistive tearing modes [21]. On the other hand, the resistivity can also provide a dissipation mechanism and, if too large, it can even damp the perturbations. As noticed at the end of Section III A, in the β = ρ = 0 regime the growth rate vanishes for both η = 0 and η → ∞. The damping effect becomes also visible in the case of infinite systems, such as those examined in [7,11], where the resistivity is stabilizing the modes.
Our work is limited by the fact that we used a simple fluid model. Finite Larmor radius effects are restricted by the conditions kρ ≪ 1 and for ρk x ≈ (ρ/L x )πn ≪ 1 so that the modes analysed are correctly described if their mode number is relatively small. Electron inertia is neglected, together with electron wave particle interactions. This prevents us from properly treating collisionless regimes which, within our model are stable but might become unstable upon reintroduction of these effects. The calculation is electrostatic an approximation that is not justified in the edge region of the fusion devices [22]. Interesting physics might occur when tearing modes are coupled to the modes described in this paper. In particular, we expect density driven (not current driven) magnetic islands, resulting from the electromagnetic version of the odd parity modes (which would have even electromagnetic flux at the resonant surface). This ∆ ′ independent modes might be related to the microtearing modes. This problem will be addressed and discussed in a subsequent article.

VII. SUMMARY AND CONCLUSIONS
We investigated the linear stability of finite size resistive inhomogeneous plasmas in a fluid approximation. The eigenmodes of the unstable drift waves, or universal instabilities, that we have analysed are standing waves generated by the reflective boundary conditions that are applied at a finite distance from the resonant surface. The study of the effect of these boundary conditions on the stability of the equilibrium was the main scope of the work presented. We found that the wave reflection provides a robust destabilization mechanism, which can persist also in relatively large systems.
The general dispersion relation characterising the modes was given in Eq.38. This expression, however, is difficult to interpret and requires simplifications in order to make its physics more transparent. Four dimensionless parameter, √ ηω * /ρ 2 , L n /L s , kρ and L x /ρ, govern the problem and determine the complex frequency of the modes.
Several limits of Eq.38 were investigated by exploiting the smallness of some of the above mentioned parameters. In particular, we were able to find exact analytic solutions in the regime L x /ρ ≪ min(L s /L n , γ 1/2 √ ηω * /ρ 2 ), which were reported in Eqs.20 and 21 and in the regime L x /ρ ≫ γ 1/2 √ ηω * /ρ 2 , described in Eqs. 22. These simple studies allowed us to determine that the system can indeed be unstable when is bounded at a finite distance L x and that its drive mechanism is the density gradient, combined by a finite resistivity. It is important to notice that, the in the limit L x → ∞, Eq.38 correctly becomes the dispersion relation derived in [7,11] for infinite systems with outgoing wave boundary conditions, which predicts stability.
For small but finite values of √ ηω * /ρ 2 , Eq.38 becomes Eq.41 and Eq.42 for even and odd instabilities respectively. These expressions are still implicit in the complex frequency and the stabilization and destabilization mechanisms are difficult to identify. Finally, by taking their limit for small L n /L s we arrived to Eq.45 and Eq.50 which describe a spectrum of unstable odd and even modes. The important effect of the magnetic shear could not be cast in a simple analytic form (but it is rigorously contained in Eqs.41 and 40) and we therefore resorted to a numerical characterisation (Sec.V B 1). In particular, we found that the growth rate of the most unstable mode is reduced by a finite magnetic shear, the critical value of which is finite and depends on the position of the wave reflection (see Eqs.54 and 55). This is probably our most important result as it shows that unconditional stability is a peculiar feature of the infinite systems and is not reproduced in more general configurations.
Under the condition that the at least partial wave reflection occurs before the mode can be damped by wave-particle interactions, we therefore expect to find unstable modes driven by density gradients in the presence of finite (albeit relatively small) magnetic shear. In other words, our conclusion is that the universal instability can exist in realistic plasmas and definitely in numerical simulations. The extension of the calculation presented in the text to a case with periodic boundary conditions is relatively straightforward. In particular, the Dirichlet boundary conditions φ(±L x ) = 0 applied to Eq.5 should be replaced with φ(L x ) = φ(−L x ). The nature of the problem would not change significantly as the energy trapping, responsible for the unstable modes, would persist with the new boundary conditions. Although we did not solve the analytical problem and obtained a dispersion relation for the periodic case, we performed several numerical simulations confirming that the unstable mechanism is present. As an example, we show in Fig.10 the unstable part of the spectrum for the same problem with Dirichelt and periodic boundary conditions (with √ ηω * /ρ 2 = 0.12, L n /L s = 0.1, kρ = 0.1 and L x /ρ = 15). L n /L s