I. The eﬀect of symmetric and spatially varying equilibria and ﬂow on MHD wave modes: Slab geometry

Realistic theoretical models of magnetohydrodynamic wave propagation in the diﬀerent solar magnetic conﬁgurations are required to explain observational results, allowing magnetoseismology to be conducted and provide more accurate information about local plasma properties. The numerical approach described in this paper allows a dispersion diagram to be obtained for any arbitrary symmetric magnetic slab model of solar atmospheric features. This proposed technique implements the shooting method to match necessary boundary conditions on continuity of displacement and total pressure of the waveguide. The algorithm also implements fundamental physical knowledge of the sausage and kink modes such that both can be investigated. The dispersion diagrams for a number of analytic cases which model magnetohydrodynamic waves in a magnetic slab were successfully reproduced. This work is then extended by considering density structuring modelled as a series of Gaussian proﬁles and a sinc ( x ) function. A further case study investigates properties of MHD wave modes in a coronal slab with a non-uniform background plasma ﬂow, for which the governing equations are derived. It is found that the dispersive properties of slow body modes are more greatly altered than those of fast modes when any equilibrium inhomogeneity is increased, including background ﬂow. The spatial structure of the eigenfunctions is also modiﬁed, introducing extra nodes and points of inﬂexion which may be of interest to observers.

The understanding of excitation and propagation of magnetohydrodynamic (MHD) waves is an important field of research in solar and plasma physics. It is thought that these waves may contribute to coronal heating due to their observed ubiquity throughout the solar atmosphere. These waves may be able to transfer sufficient amounts of energy to the upper layers of the solar atmosphere and undergo dissipation, consequently heating up localised plasma. Dissipation processes may include, but not be limited to, resonant absorption (Hollweg & Yang 1988;Ruderman & Roberts 2002;Howson et al. 2019), Alfvén wave phase mixing (Heyvaerts & Priest 1983;Ofman & Aschwanden 2002) and transformation of waves via parametric processes (Voitenko & Goossens 2002;Fedun et al. 2004;Voitenko & Goossens 2004;Vásconez et al. 2015). Complete understanding of the properties of these waves including their dispersive nature is essential in determining their contribution to the energy budget of the solar atmosphere.
The dispersion diagram provides a useful tool in understanding the properties of various MHD waves. It indicates the permittable phase speeds at which trapped waves propagate and provides information about their dispersive characteristics. Previously, there ⋆ E-mail: sjskirvin1@sheffield.ac.uk has been a number of studies devoted to the analysis of MHD wave properties in different magnetic configurations. Roberts (1981) and Edwin & Roberts (1982) investigated MHD waves in Cartesian geometry modelled as a magnetic interface and a magnetic slab respectively. The magnetic slab has been investigated in both a magnetic and a field free environment. A visual representation of the magnetic slab is shown in Figure (1). These theoretical studies have provided the fundamental models which describe wave propagation for magnetic waveguides found in the solar atmosphere. Further extensions of these models were provided by introducing steady background plasma flows (Nakariakov & Roberts 1995a), linear background flows (Zaqarashvili 2011) and also moving to cylindrical geometry (Edwin & Roberts 1983), with magnetic twist (Erdélyi & Fedun 2007, 2010 and curvature of the waveguide (Verwichte et al. 2006a,b), to name but a few. Generalising the traditional slab model of a solar waveguide to a more realistic case has been an extensive area of research. For example, coronal loops have been modelled in planar geometry with smooth density profiles (Lopin & Nagorny 2015a) although these studies consider a continuous background plasma profile and do not match any boundary conditions, simply because there are no distinct boundaries in the model. The period ratio of the fundamental mode to twice the first overtone P 1 /2P 2 was investigated by Macnamara & Roberts (2011)  found that there were striking similarities to that of the step function considered traditionally, suggesting that the step function profile for density may actually be a useful model for wave investigations.  derived a generalised dispersion relation for fast waves in a cold coronal slab with a finite plasma-β (ratio of plasma pressure to magnetic pressure). However, in this investigation continuous transverse profiles are considered split into three regions where the middle regime is allowed to be arbitrary. Previous works have developed a number of techniques to numerically solve the differential equation with carefully chosen profiles for plasma density or Alfvén speeds (Oliver et al. 1993;Verwichte et al. 2006a;Soler et al. 2017;Thackray & Jain 2017). These studies also focus on coronal structures, namely coronal loops, where plasma pressure is assumed negligible and the plasmaβ is zero. The cold plasma-β limit restricts these investigations to fast magnetoacoustic waves only as it is assumed that the Lorentz force dominates. The behaviour of MHD waves propagating in coronal loops with specific inhomogeneous density profiles has also been studied before in both planar and cylindrical geometries, see e.g. Edwin & Roberts (1988); Nakariakov & Roberts (1995b); Lopin & Nagorny (2015a,b); Li et al. (2018). The specific choice of the density profile could result in an analytical derivation of the dispersion relation. However, in the majority of these cases a continuous density profile is considered, such that there are no boundaries of the waveguide and as a result, no boundary conditions to be matched. Studies of wave behaviour in the context of different realistic background plasma profiles are important as they provide information on group velocity as a function of density structure (Edwin & Roberts 1988), and, therefore, on how the waves may transfer their energy. Whilst this research has greatly advanced understanding of wave phenomena in coronal structures, the dispersion relation is derived in all cases through carefully chosen plasma profiles allowing an analytical solution to be obtained. To model wave propagation in a realistic magnetic configuration with a spatially varying plasma equilibrium in the presence of background flows, a numerical approach has to be used. Very recently Claes et al. (2020) have developed a numerical code to solve the full MHD spectrum for any given 1D equilibrium. In the approach presented by the authors a finite element method (FEM) was implemented to identify the permittable eigenvalues.
In this paper, the numerical methodology for obtaining solutions located on the dispersion diagram based on the shooting method approach will be introduced in Section 2 to a case of a magnetic slab in Cartesian (planar) geometry. In Section 3 the results of the proposed method are directly compared to previous studies of wave propagation and analytically derived dispersion diagrams. These comparisons are done under the context of both static and steady background plasma flows. In Section 4 the developed method will be applied to inhomogeneous density profiles which take the form of a Gaussian distribution and a sinc(x) function, better known as sin(x)/x, which have not been previously investigated and can not be solved analytically. Finally in Section 5 the properties of MHD waves in a coronal slab in the presence of a Gaussian non-uniform background flow is investigated.

METHOD
The set of ideal, compressible MHD equations used in the investigation are: Figure 1. The magnetic slab model traditionally used to model features observed in the solar atmosphere. Red arrows represent magnetic field lines and intend to show the different magnetic configurations in each regime internal and external to the waveguide. Blue arrows represent the vertical equilibrium velocity field which may be non-uniform in spatial coordinate x.
The slab is assumed to be homogeneous and unbounded in the y-dimension. Parameters ρ, P, U 0 and B denote plasma density, pressure, background flow and magnetic field respectively. Subscripts i and e relate to internal and external properties respectively. The special scenarios in this work investigate a spatially varying density profile ρ 0i (x) and vertical background flow U 0i (x). In reality any/all equilibrium background plasma variables are allowed to be inhomogeneous and the numerical approach presented is capable of handling such a scenario. The specific spatial profiles considered in this work are shown in Figures (6), (12) and (16) for inhomogeneous equilibrium density and also Figure (21) for an inhomogeneous background flow.
where ρ is the plasma density, v is the velocity field, p is the plasma pressure, B is the magnetic field, µ 0 and γ are the magnetic permeability of free space and the ratio of specific heats respectively.
∂ ∂z is the total (material) derivative. By linearising the system of Equations (1)-(5) and seeking wave solutions of the form e i(kz−ωt) where k is the wavenumber and ω its frequency, a governing differential equation can be derived for either the transverse velocity perturbationv x or displacement component,ξ x and also an expression for the total pressure perturbationP T . Total pressure is defined as plasma plus magnetic pressure given by p+ B 2 /2µ 0 . After matching boundary conditions, namely the continuity ofv x andP T , a transcendental dispersion equation for ω and k is retrieved. The solutions of such equation will give the permittable 'trapped' physical wave modes which can propagate through the waveguide. This approach previously works due to the analytical solution of the governing differential equation being known, however, breaks down when the format of the Magnetoacoustic waves in inhomogeneous waveguides 3 differential equation becomes more complicated and no known analytical solution exists. Such reasons that the differential equation inv x may become more complicated could be from adding extra inhomogeneities into the original equilibrium, such as not assuming the plasma properties are uniform in space and unchanging in time.
The physical nature, such as size and shape, of the inhomogeneities are important as they influence propagation of MHD waves along waveguides and may significantly effect their dispersive nature. The plasma properties inside and outside the slab are assumed to be known, dependant upon whether photospheric or coronal conditions are being investigated. The general differential equation takes the following form for a magnetic slab: where F and G are arbitrary functions, c 0 , v A and c T are the sound, Alfvén and tube speeds respectively. A prime denotes a derivation with respect to the spatial coordinate x. Note that equation (6) may be with respect toP T also. The functions F and G are usually constant or zero in previous studies (see e.g. Roberts 1981;Edwin & Roberts 1982;Nakariakov & Roberts 1995a) due to the uniform plasma being investigated. However, these functions become more complicated and non constant if, for example, an inhomogeneous plasma is considered. When the differential equation is with respect to horizontal velocity perturbation then the corresponding total pressure perturbation expression for a magnetic slab takes the form: with A and B again arbitrary general functions. In this paper, a numerical method to obtain the dispersion diagram is presented without requiring the corresponding dispersion relation, usually obtained analytically. This method takes advantage of the shooting method, a numerical analysis technique designed to solve a boundary value problem (BVP) by reducing it to an initial value problem (IVP). It also requires that all quantity profiles are symmetric about the centre of the waveguide.
Let's consider a magnetic slab in two-dimensional planar geometry such that the internal plasma parameters are different to the external plasma parameters but both regimes are uniform. There are certain conditions which need to be met to identify the permittable wave modes in this model. Firstly, the waves are required to be laterally evanescent outside the waveguide, such that the local wave amplitude tends to zero infinitely far from the waveguide, confining waves to exist only within the slab. Secondly,v x andP T must be continuous at the boundaries of the waveguide. To satisfy the first condition, it can be assumed that the wave amplitude,v x is extremely small far away from the slab hencev ′ x must too be very small.
Numerically, it is not required to consider the point at infinity so long as the domain is restricted to be considerably larger than the width of the waveguide. The location of infinity must be chosen such that it supports a sufficient number of wavelengths for the wave amplitude to decay. The boundaries of the slab provide a BVP, however, here a BVP solver cannot be used because there is not enough information aboutv x orP T at the boundary. In the shooting method, the value of either function at an initial point is known, with the aim to effectively guess the format of the function derivatives until a solution is found that goes through the second boundary value. Fortunately, under the assumption that the perturbations decay far away from the waveguide, it is known that at a location covered by enough wavelengths bothv x andv ′ x are very small, almost zero. Therefore using these initial conditions, Equation (6) can be solved numerically exterior to the slab up to the boundary. By assuming the solutions are evanescent outside the waveguide, in this paper we seek trapped modes and ignore leaky waves which may undergo damping due to wave leakage away from the waveguide (Stenuit et al. 1998). Focusing solely onv x for the moment, the condition restrictingP T will be implemented later.
The solution forv x is found numerically exterior to the slab up to and including the boundary using a numerical solver for ordinary differential equations (ODEs). Therefore, the second order differential equation inv x is reduced to first order by introducing a new variable and the initial values of bothv x andv ′ x are assumed very small (≪ 1). Solving this set of equations provides values at each x location for bothv x andv ′ x . At the boundary, a value ofv x has now been obtained, however, the value ofv ′ x inside the slab is not known. Therefore, the shooting method is implemented to 'shoot' to the opposite boundary depending on what wave mode is being investigated, for example sausage or kink.
Magnetoacoustic wave modes are split into sausage and kink depending if the boundaries of the waveguide oscillate out of or in phase with each other, respectively. The sausage mode is the antisymmetric solution about x = 0 and the kink mode is the symmetric solution forv x . These wave modes can be further separated into surface and body respectively, if the wave amplitude exists only at the boundary or oscillates throughout the waveguide (Edwin & Roberts 1982;Priest 2014). Therefore, enough information is known at specific locations in the domain to find solutions for either the kink or sausage modes by applying the aforementioned conditions. If the sausage mode is being investigated, the algorithm is design to 'shoot' from the one boundary to the other requiring the values to be anti-symmetric. When investigating the kink mode, the algorithm finds a solution on the opposite boundary with a value ofv x equal to the known boundary value. Assuming symmetry of the model in this way means that the solution does not need to be found in the other external region. The method of shooting for the internal region, where in this paper the plasma density and background flow is allowed to be non-uniform, reduces Equation (6) to first order and obtains values forv x andv ′ x at every spatial location. It should be noted that in a non-uniform plasma Equation (6) can possess singularities which correspond to solutions lying inside the slow/Alfvén continuum.
The total pressureP T , Equation (7), must also be continuous across the slab boundaries and an ω, k pair will only exist if both conditions are met. The equation describingP T can be derived from the momentum Equation (2) and contains a term multiplied byv ′ x . When reducing the differential Equation (6) to first order, the solution forv ′ x is obtained trivially at all locations in the domain. Therefore, this value can be used to calculate the value forP T at the boundary for both the internal and external plasma parameters of the slab. A match forP T is obtained if the internal and external values are within a given small accuracy parameter ǫ of each other. If this is the case then the corresponding phase speed, from the ω, k solution, is plotted on the dispersion diagram. It is highly unlikely that an ω, k pair from the initial sampling will be an exact solution, therefore a numerical bisection method is adopted alongside the shooting method. The wavenumber is fixed and the algorithm iterates through different wave frequencies. The value of the internal eigenfunction solution at the boundary and the external eigenfunction solution at the boundary for each frequency is stored. The value of eigenfunction difference at the boundary changes sign when an exact solution is passed (i.e. goes from positive to negative). When this condition is true, the algorithm goes back to the previous frequency sampled and the step size is reduced in order to locate the exact eigenvalue. The solutions are split into their separate branches on the dispersion diagram and fitted with a high order polynomial curve.

Magnetic slab case
In this section the results of the method will be compared to results previously obtained by Edwin & Roberts (1982) in which a magnetic slab embedded in a stationary magnetic environment is investigated. The graphical representation of this problem is shown in Figure 1. The analytic dispersion relation is derived with a closed form solution by assuming that all plasma quantities including sound and Alfvén speed are uniform in all regions. The known analytic dispersion relation is given by Equation (11) in Edwin & Roberts (1982) and solutions of this transcendental equation give the resulting dispersion diagram. Here, the functions F and G are constant and therefore Equation (6) can be written as: where The corresponding total pressure perturbation given by Equation (7) can be written as: which is proportional to the derivative in the velocity perturbation.
Equations (8) and (9) providing the two boundary conditions which require to be continuous. All variables outside of the slab will have indexes e but take the same form as Equations (8) and (9). Following the procedure outlined in Section 2 the shooting method is applied to solve equation (8). A solution will only be obtained if a value of frequency and wavenumber simultaneously satisfies Equations (8) and (9) forv x andP T . It can be seen in Figure 2a that the obtained solutions for surface waves and body waves fits well those results found previously by Edwin & Roberts (1982, see Figure 3) under photospheric conditions (v Ae < c i < c e < v Ai ). Typically, obtaining solutions for body modes, located between c T i < v ph < c i (v ph = ω/k) for the photospheric case, tends to be more difficult due to the reduced step size required when numerically solving the dispersion relation by the bisection method, however this method finds the exact solutions with no alterations in the algorithm. More MHD wave modes could be retrieved by increasing the number of samples in the domain, however this comes at a cost of increased numerical intensity, so a balance must be found. This also includes higher harmonics of body modes which will be retrieved also with increasing resolution.
The results of the calculation for the scenario of a magnetic slab under coronal conditions (c e < c i < v Ai < v Ae ) are shown in Figure 2b and can be compared to Figure 4 in Edwin & Roberts (1982).
Similarly to the photospheric case the phase speed solutions are recovered well for both the sausage and kink modes. Under coronal conditions as stated in Edwin & Roberts (1982) only body modes exist, indicating here the power of this method to recover body mode solutions with no additional steps required in the algorithm. It can be seen only specific branches of the fast body waves are recovered, presumably these are the first harmonics and higher harmonics will be identified with greater resolution.

Magnetic slab with steady flow
The magnetic slab model can be further extended to investigate observed features by the inclusion of a background steady flow. Previous studies have included the addition of a steady background plasma flow, see e.g. Nakariakov & Roberts (1995a); Zaqarashvili (2011); Ebadi et al. (2011). It is well known that the introduction of a background flow into the model increases the amplitude of the wave perturbations. The flow also gives the waves an observed phase shift when compared to the no flow magnetic slab model, proportional to the speed of the flow. Furthermore, introducing a background flow supports development of the Kelvin-Helmholtz instability at the boundary of the flux tube.
In the presence of a steady background plasma flow, the structure of Equation (6) remains the same, but the wave frequency is now phase shifted by a magnitude proportional to the plasma flow speed U 0 , that is ω becomes ω − kU 0 . As a result the coefficient in Equation (6) becomes: and where Ω is the phase shifted frequency. Now that a flow has been introduced, two conditions at the boundaries of the slab still required to be satisfied. Whereas before, the continuity of total pressure and horizontal velocity perturbation were the conditions, the latter is replaced by the continuity of the horizontal displacement perturbation. This is due to the additional background flow U 0 , which can locally amplify the wave displacement at the boundary and needs to be accounted for. Therefore the new boundary conditions which need to be satisfied are given by: as explained in Nakariakov & Roberts (1995a) where ±x 0 are the locations of the slab boundaries (see Figure 1). Equations (6) and (7) along with condition (10) are used in order to attempt to retrieve ω and k pairs on the dispersion diagram. Nakariakov & Roberts (1995a) have found that the properties of magnetoacoustic waves in a magnetic slab with a background flow were not drastically different to those of the static model, although for specific values of external flow, the dispersion diagram was slightly different as certain wave modes disappeared from the diagram. Figure 4a shows the results of the methodology under photospheric conditions with a downward steady flow external to the slab. The phase velocity axis is normalised relative to the internal Alfvén speed so that direct comparison can be made with the results of Figure 2b in Nakariakov & Roberts (1995a). The obtained results show a good agreement of those first retrieved by Nakariakov & Roberts Magnetoacoustic waves in inhomogeneous waveguides 5  (1995a) and it is possible to sample regions between specific characteristic speeds such that higher resolution can be obtained in the regions cut off by the presence of the flow. The corresponding coronal solutions are shown in Figure 4b with a steady internal flow of v 0 = 0.35v Ai and no external flow. Again, the same solutions are obtained as those shown in Figure 1c in Nakariakov & Roberts (1995a) including the backward propagating body modes at small v ph .
The next section extends these uniform cases and adds inhomogeneity into the model by changing the form of Equations (6) and (7).

NON-UNIFORM DENSITY PROFILES
The properties of MHD waves in a non-uniform plasma have been investigated before and the general formalism for any 1D inhomogeneity, including stratification along the direction of non-uniformity, has previously been presented in (see e.g. Goedbloed et al. 2010;Goedbloed et al. 2019;Roberts 2019).
In this section, a spatially dependant internal plasma density is introduced into the model. This significantly changes the expression given by Equation (6)   v where: The equivalent expression forP T is given by equation (9) where all quantities are now a function of x. Consider a magnetic slab embedded in a stationary, uniform and magnetised environment under coronal conditions (c e < c i < v Ai < v Ae ). Background plasma flow is ignored and the characteristic speeds are chosen to match those given in the coronal case in section 3.1. Inside the magnetic slab, a Gaussian profile in plasma density is introduced, see Figure 6, given by the expression: where x c is the centre of the Gaussian located at x = 0, W is the standard deviation (i.e. the width) of the density distribution and ρ 0i is the maximum internal density given by the value in section 3.1. The other equilibrium density structuring investigated in this work considers a sinc(x) function. The motivation behind modelling a sinc(x) profile comes from observations of magnetic bright points (MBP's) which have been observed to have spatial intensity distributions similar to such a profile (see e.g. Jess et al. 2010). To fit the numerical domain, the normalised sinc(x) function is modelled in this work using: for a coronal slab and is normalised such that the maximum is comparable to the uniform slab case. Pressure balance is accounted for by a change in temperature inside the slab. The width of the Gaussian profile determines the gradient of the inhomogeneity. The characteristic speeds are therefore also spatially dependant. This specific case corresponds to a cool magnetic flux tube.

Inhomogeneous magnetic slab under coronal conditions
The solar atmosphere is highly inhomogeneous and types of inhomogeneity could arise from non-uniform density and magnetic field structuring, or unsteady flows. Investigating the trapped wave modes of a solar waveguide within a non-uniform background plasma is relevant to study from a theoretical point of view such that models can be created which more accurately represent those seen in observations. Typically this investigation is done analytically, however, when plasma variables are non-uniform in space, the governing MHD equations become more complicated to solve. Some specific plasma profiles have been previously extensively investigated, chosen such to allow the derivation of an analytic dispersion relation and as such retrieve solutions on the dispersion diagram. A review of the density profiles which have been studied before are given in Table 4 of Li et al. (2018). Here, the proposed numerical approach will be applied to investigate density profiles which can not be analysed analytically. For a large Gaussian width, the inhomogeneity inside the slab is weak and the plasma is similar to the uniform case described in Section 3.1, therefore corresponding to the same results as the uniform case. This case is shown in Figure 7a, where the characteristic speeds at the boundary (subscript 'B') correspond to the uniform speeds in Figure 2b. The width of the profile here is chosen to be W ≫ 2x 0 , i.e. many orders of magnitude larger than the width of the waveguide. By changing the width of the inhomogeneity (i.e. the standard deviation W), the value of ρ 0i decreases at the boundary compared to the uniform case and alters the trapped modes of the system. Seen in Figures 7b, 7c, 7d are the resulting dispersion diagrams for the density profiles in Figure 6. The fast body solutions are still bounded between v Ae and v Ai however are cut off by the internal Alfvén speed at the boundary for a large Gaussian inhomogeneity (W < 1) as it is not possible for trapped global eigenfrequencies to enter this continua. Slow body waves in an inhomogeneous coronal slab are bounded between the tube speed at the boundary of the waveguide (c T B ) and the maximum internal sound speed due to the density structuring for small inhomogeneity. This can be seen in Figures 7b and 7c, unlike the uniform case, where the slow body waves are trapped between c i and c T i . Furthermore, the band located between c T B and maximum c T i is a resonant region known as the cusp (slow) continuum. This band is due to the singularity in Equation (11), when v ph = c T i which provides great interest as resonant absorption can occur here and has been subject to previous analytical investigation. Firstly by Keppens (1996) who studied this effect in a cylinder with an unmagnetised surrounding and later by Yu et al. (2017a,b) in a photospheric slab with a weakly magnetised surrounding. It is worth further noting that the work by Keppens (1996) also investigated the leaky modes, which can radiate energy away from the waveguide. This study however is beyond the scope of the current paper and will be investigated in future work. Turning focus now to the spatial profile modelled as a sinc(x) function. This profile is shown by the blue line in Figure 4. The value of ρ 0i at the boundary for the sinc(x) profile is similar to that of a Gaussian profile with W = 1.5, however the structuring inside the slab is much different. The dispersion diagram shown in Figure  9 has similar characteristics to those shown in Figure (7c) with slight change in the positioning of slow body modes such that the majority of these branches now lie inside the green shaded region which represents the spatial sound speed band.
Comparisons of the eigenfunctionsP T andv x for all possible modes are shown in Figure 11. For the fast modes, both sausage and kink, equilibrium inhomogeneity has a minor effect on the physical properties of the wave mode. It can be seen in Figure 10a that as the Gaussian inhomogeneity is increased, the anti-nodes of the fast sausage mode shift towards the center of the waveguide, an effect which has been shown in coronal loop analysis by Verth et al. (2007). The amplitude of the total pressure perturbation is also locally increased at the centre of the waveguide as the Gaussian inhomogeneity is increased. Figure 10b indicates the nodes of the total pressure perturbation become more pronounced as the inhomogeneity increases, a similar albeit more minor effect can be seen in the spatial structure of thev x eigenfunction. It is worth mentioning here that the profile proportional to sinc(x) does not appear to affect the physical spatial distribution of the eigenfunctions for fast modes in a coronal plasma, suggesting that these modes may not be a suitable choice to use for spatial coronal-seismology.
However for the slow modes, the inhomogeneity has a much greater effect. Figure 10c shows the perturbed eigenfunctions for the slow body sausage mode. It is obvious that increasing the inhomogeneity away from a uniform plasma has a clear effect on the physical properties of this mode. Decreasing the Gaussian width creates extra nodes and anti-nodes in the resultingv x perturbation, these extra anti-nodes may be misinterpreted in observational data as an entirely different mode. Increased inhomogeneity also has an effect onP T , changing the center of the waveguide from having a maximum to having a minimum at this location for the slow sausage mode. The slow body kink mode is also greatly affected by small inhomogeneity compared to the uniform slab. As Gaussian inhomogeneity is increased, the maximumv x perturbation is achieved closer to the boundaries of the waveguide, rather than obtaining a single maximum at the centre in the uniform scenario. The total pressure perturbation shown in Figure 10d is still zero at the center of the waveguide as expected for the kink mode however displays an anti-symmetric behaviour to the uniform scenario as the Gaussian profile in density becomes more pronounced. The slow body modes in a coronal plasma may be a good indicator into the underlying plasma density structure. Thev x distributions in Figure 10c appear to be proportional to the derivative of the corresponding equilibrium density profiles shown in Figure 4.

Inhomogeneous magnetic slab under photospheric conditions
For photospheric conditions (i.e. v Ae < c i < c e < v Ai ), the density profiles investigated are shown in Figure 12 to model an evacuated magnetic slab. Under these conditions, the Gaussian density maintains the same expression as that for the coronal case, however the sinc(x) function must be normalised to the numerical domain and now takes the form: This model may better represent those conditions found in sunspot umbrae and penumbrae, with a continuous internal density profile. The case when the width of the profile is large compared to the width of the waveguide is shown in Figure 13a and is comparable to the uniform case shown in Figure 2a as expected. Adding extra density inhomogeneity into the internal region alters the plasma properties at the boundary compared to the centre. Figures 13b, 13c and 13d show this and the shaded regions denote the area covered by the inhomogeneity for each characteristic speed. Slow surface waves are trapped at speeds below c T B and above v Ae . This result is expected from theory presented by Edwin & Roberts (1982). Shown in Figure 13c the value of c T B surpasses v Ae and as a result slow surface waves cease to exist. An interesting region to note is the area contained within c B < v ph < c i . As extra inhomogeneity is added into the equilibrium, this region becomes larger. Due to the presence of the boundary layer in the model provided by the width of inhomogeneity, this region varies continuously between the boundary c B and minimum value of internal sound speed c i . Likewise compared to the coronal case this introduces Alfvén and cusp resonant singularities into the differential equation and provides the possibility for dissipation processes to occur. Figure 15 is the resulting dispersion diagram for a photospheric equilibrium with a density structure modelled as a sinc(x) function. The algorithm finds eigenvalues very similar to the Gaussian scenario with a width equal to 0.9 as shown in Figure 13d. Interestingly, the shape of the equilibrium inhomogeneity, does not appear to have a great affect on the eigenvalues of the equilibrium system. Similar to the coronal case this is true for fast modes, which may not feel the inhomogeneity as much as slow modes as they are able to travel across magnetic field lines and as such travel more freely across any inhomogeneity. Figure 18 shows the perturbed eigenfunctions for the fast surface sausage and fast surface kink mode under photospheric conditions. It is clear that inhomogeneity in the form of a Gaussian or a sinc(x) function has very little effect on the physical behaviour of the total pressure and velocity perturbation.
Due to the decreasing density ratio at the boundary for Gaussian widths with larger inhomogeneity, the cusp continuum band becomes larger, therefore cutting off the slow surface and body modes. To further investigate these further, we consider initial density profiles with a similar Gaussian width but smaller inhomogeneity, see Figure 16.
Displayed in Figure 20 are the perturbed eigenfunctions for body sausage and body kink modes under photospheric conditions and equilibrium density structure shown in Figure 16. Again, it is clear that any equilibrium inhomogeneity, even minor changes in the form of a Gaussian profile, has a significant effect on the properties of slow modes. Thev x perturbation for the sausage mode has antinodes which shift again towards the centre of the waveguide, an affect mirrored by theP T perturbation of the kink mode, due to the asymmetry of the kink and sausage modes.

CORONAL SLAB IN PRESENCE OF INHOMOGENEOUS BACKGROUND FLOW
In this section, an analysis of magnetoacoustic wave properties in the case of a magnetic slab of uniform plasma in the presence of a non-uniform background flow will be conducted. The governing equations are derived under the context that the background plasma flow is symmetrically-arbitrary and spatially varying. This scenario is very applicable to the solar atmosphere. Jet-like features and waveguides with a background flow are routinely observed in the solar atmosphere, including spicules, fibrils and prominences to name a few (de Pontieu et al. 2007;Berger et al. 2010;Pereira et al. 2011). Due to current spatial resolution limits of ground and space based telescopes, the spatial profiles of these flows are still unknown however it is common in fluid dynamics that flows are not steady (e.g. Orszag & Kells 1980  is shown in Figure (1). In our model, we assume that the magnetic field aligned plasma flow is present only inside the magnetic slab e.g. v 0i = (0, 0, U 0i (x)). The solution for the case exterior to the slab is the same as that shown in Section 3.1. The general formulation for a non-uniform equilibrium, including background plasma flow, has been previously derived in Frieman & Rotenberg (1960); Goedbloed et al. (2019). Here we present a specific case of Equation (13.11) from Goedbloed et al. (2019), in which the plasma properties are uniform expect for a symmetric inhomogeneous background plasma flow, which unlike previous studies is allowed to be discontinuous across the slab boundary. This chosen configuration of an inhomogeneous vertical background flow aligned with the magnetic field eliminates any magnetic shear effects, including those associated with the background flow.
The set of linearised, Fourier-decomposed MHD Equations (1)-(5) for each perturbed quantity in the presence of an inhomogeneous background plasma flow described above are: −iΩ(x)P 1 + c 2 Density profile modelled as a Gaussian with a varying inhomogeneity in a photospheric slab. Width of inhomogeneity given by W = 10 5 (black), W = 3 (yellow), W = 1.5 (green), W = 0.9 (red). A spatial profile proportional to a sinc(x) function (blue) is also modelled. In all cases the density is discontinuous at the waveguide boundary and tends towards ρ e (shown) at the boundary. The boundaries of the slab are indicated by the red dashed lines.
where Ω(x) = ω − kU 0i (x) is the spatially Doppler shifted frequency and a prime denotes a differentiation with respect to the spatial coordinate x.
Equations (12)-(17) can be combined to eliminate all perturbed quantities other thanv x and as such derive a governing equation for velocity amplitude: where, where, with Ω 2 A (x) = k 2 v 2 Ai − Ω 2 (x). Equation (18) has no known closed form analytical solution, due to it's complicated nature caused by the spatially varying coefficients. It can be seen that if the spatial dependence on flow is removed (i.e. an initial constant flow) that Ω(x) → ω − kU 0i = Ω, as such m 2 i (x) is also no longer a function dependant on space. Furthermore, Ω ′′ (x) = Ω ′ (x) = 0 such that Equation (19) becomes equal to zero. Equation (18) now becomes: which is the same result as shown in Equation (8) but with a Doppler shift with respect to the flow velocity for the steady flow scenario. Likewise, if magnetic field is neglected (B = 0) such that v Ai , c T i = 0 and assume an incompressible plasma (c i → ∞), where in this limit the sound speed becomes unbounded, the governing Equation as given in Timofeev (2000). Equation (21) is a form of Rayleigh's equation (see e.g Chandrasekhar 1961), a well known previously obtained expression in hydrodynamics concerning inviscid shear flows (Rayleigh 1879; Hirota et al. 2014). The corresponding expression for the total pressure perturbationP T inside a magnetic slab with a spatially varying background plasma flow is given by: The set of Equations (18) and (22) provides the required expressions for the numerical analysis. Similar to the analysis conducted in Section 4, an internal Gaussian flow profile is considered with the spatial profiles analysed shown in Figure 21 given by the expression: where A is the velocity amplitude, x c = 0 is the centre of the slab, W is the width of the Gaussian and assuming that there is no plasma flow outside the slab. The resulting dispersion diagrams for a magnetic slab under coronal conditions with an inhomogeneous Gaussian internal background plasma flow for certain cases of Figure 21 are shown in Figure 23. As the maximum of the Gaussian flow located at x c = 0 remains constant for all cases, the maximum Doppler shift on the waves in all cases also remains constant. There is not much difference in the dispersive properties between the fast body modes -both forward and backward propagating, as the inhomogeneity is increased. The dispersive properties of the slow body modes however, are much more affected by the nonuniformity of the background plasma flow. The region bounded between −c T i + U B < v ph < −c T i + U 0i is a band in which flow related resonances are present. This region is not shown in Figure  23 and although this phenomena is not investigated further in this paper, it has been studied before analytically (Taroyan & Erdélyi 2002).
Displayed in Figure 25 are enlarged regions of the dispersion diagrams in Figure 23 to display the region of forward and backward propagating slow body modes in detail. Seen in Figure 24a  The eigenfunctionsP T andv x in the presence of an inhomogeneous background flow are shown in Figure 27. Similar to the study of inhomogeneous density profiles, the equilibrium inhomogeneity appears to have little effect on the spatial properties of fast propagating wave modes. Figure 29 are the eigenfunctions for the slow forward and backward propagating sausage and kink modes for the model shown in Figure 21. Figures 28a and 28b show the spatial behaviour of the eigenfunctions for a forward propagating slow body mode.

Shown in
The perturbation ofv x corresponding to the sausage mode obtains an extra point of inflexion caused by the inhomogeneity of the background plasma flow. For both thev x perturbation which corresponds to the sausage mode and theP T perturbation corresponding to the kink mode, the maximum value of the eigenfunction antinode shifts towards the centre of the waveguide, where the flow speed is a maximum and it's gradient is zero. If it were possible to model a Gaussian flow stretching to infinity, the corresponding spatial eigenfunctions would become discontinuous at the centre of the waveguide.

CONCLUSIONS
In this paper we have developed a generalised mathematical model that represents some possible equilibria cases in the solar atmosphere. A new numerical approach, based on the shooting method, has been developed for obtaining the eigenvalues and eigenfunctions  for magnetoacoustic wave modes present in different equilibria and allows for an arbitrarily symmetric spatial plasma structuring inside the waveguide. The numerical method is utilised along with the physical properties of MHD sausage and kink wave modes. The technique has been tested by comparison with known results in e.g. a magnetic slab under different magnetic regimes. Furthermore, a constant steady field aligned plasma flow is introduced for which the known solutions are obtained under both photospheric and coronal conditions. The method is then applied to a waveguide in which the internal plasma structuring is modelled as a Gaussian profile. The analytical function used to describe this density distribution cannot derive a dispersion relation analytically. It is found that the cut off values for slow body modes are dependant upon the values of the sound speed at the boundary c B and cusp speeds c T B and the size of internal inhomogeneity in both coronal and photospheric conditions. Also, this technique can be useful to identify the bands in which resonance can occur and potentially lead to dissipation processes such as resonant absorption and phase mixing. The new properties of MHD wave modes in a non-uniform magnetic slab have been studied. Whilst fast surface and body modes are not greatly affected by the equilibrium inhomogeneity, the physical eigenfunctions of slow body modes are significantly altered. We have found that even for an equilibrium Gaussian density distribution of sufficient width, an extra node and point of inflexion appears in the spatial structure of the eigenfunctions. This result is important for the interpretation of observational signatures of MHD wave modes as it shows non-uniform equilibria may lead to the misunderstanding of their spectral patterns.
To further show the strength of this technique, an example of a coronal slab with a non-uniform background flow was investigated. The general governing equations for the perturbations of total pressure and velocity were derived, which reduce to the known expressions when inhomogeneity is ignored. These equations were then solved numerically to obtain eigenvalues and plot the dispersion diagrams for a number of different non-uniform background flow profiles. The non-uniform flow created an asymmetry between the bands in which forward and backward propagating slow body waves are trapped within the waveguide. A similar behaviour to a non-uniform plasma density is observed with the resulting eigenfunctions also. As the inhomogeneity of the background flow is increased, the spatial structure of slow body waves -both forward and backward propagating -become distorted, while the fast modes remain unaltered.
The major benefit of the presented methodology is that more complicated plasma equilibrium can be introduced into the original slab model. Non-uniform plasma flow profiles which more accurately reflect those observed in e.g. sunspots can be included and the resulting wave analysis can be conducted. It is well known that the resonant absorption and phase mixing mechanisms rely upon the presence of a non uniform plasma or inhomogeneous boundary layer -for which this method could provide a better understanding of in terms of wave properties. This would previously not have been achievable due to the complicated mathematics involved and the numerous simplifications and assumptions needed to be able to obtain an analytical solution to the MHD equations.
Immediate future steps involve extending the current work to a cylindrical model, which perhaps models solar features such as sunspots, jets and coronal loops more appropriately. This will be possible using the introduced methodology as the physics remains the same however the differential operators take a different form in cylindrical geometry. Further extensions could include modelling asymmetric profiles, for example non-steady asymmetric background plasma flows, in which the physics for MHD wave modes would be modified slightly. Another very important aspect which can be studied is the behaviour of complex wave frequencies.
Complex frequencies would provide more information on the physical behaviour of wave damping, leaky modes and any instabilities that are present in a static or steady equilibrium. Implementation of asymmetry, extra geometries and complex frequencies will help move a step closer to developing a technique which can be used in wave analysis for a general arbitrary model. The power of this being that it can be used alongside observational data to conduct magnetoseismology with a realistic model to best fit observational results and provide a greater insight into the physical properties of