# Solving a complicated, non-polynomial dispersion relation

I am quite new to Mathematica, so I apologise if the following question is very basic; I can’t seem to find an appropriate answer anywhere!

I am trying to solve a quite complicated dispersion relation (see below) for some frequency $$\omega$$ which, in general, may be complex, and is a function of the variables $$k_x$$, $$k_y$$, $$k_z$$, as well as the numerical constants $$T$$, $$B$$ and $$\tau$$. I would like to be able to plot the real and imaginary parts of the solution $$\omega$$ for a large range of values of $$k_x$$, $$k_y$$ and $$k_z$$ given the values of $$T$$, $$B$$, $$\tau$$, e.g., create a contour plot of $$\omega$$ in the $$(k_x,k_y)$$ plane.

When I have been tackling simpler problems of this type with Mathematica, I have been using Solve[(* dispersion relation *) ==0, [\omega]] to obtain [\omega] = [\omega][kx,ky,kz,T,B,\[tau]], which then will allow me to produce the desired output by simply plotting this function with respect to the given variables for values of $$T$$, $$B$$, etc. However, I immediately run into a problem in the more complicated case that I am considering, as Mathematica is unable to solve the dispersion relation using Solve or Reduce given its complexity, and non-polynomial nature.

What is the best way to approach this problem while keeping the solution method sufficiently general that I can still plot $$\omega$$ as a function of $$k_x$$, $$k_y$$ and/or $$k_z$$ for various values of $$B$$, $$T$$ and $$\tau$$? From what I have gathered by searching online, an appropriate way would be to use NSolve to create a table for $$\omega$$ for various values of $$k_x$$, $$k_y$$, and $$k_z$$. However, I would like to also retain the dependence numerical constants $$T$$, $$B$$ and $$\tau$$, but I can’t see a good way of doing this without simply trying to use Solve, as above.

Any help would be greatly appreciated!

For completeness, the dispersion relation is as follows: \begin{align} 0=&\zeta^2\left\{k_x^2 + k_y^2 -2\left( \zeta^2 – \frac{1}{2} \right)\left[1+\zeta \mathcal{Z} – \zeta \zeta_* – \zeta_* \left( \zeta^2 – \frac{1}{2} \right) \mathcal{Z}\right]\right\} \ &- \left[ \frac{1}{2}(k_x^2 + k_y^2) – 2 \zeta_d \zeta_* \right]\left\{ 1- 2\tau\left( \zeta^2 – \frac{1}{2} \right)\left[1+\zeta \mathcal{Z} – \zeta \zeta_* – \zeta_* \left( \zeta^2 – \frac{1}{2} \right) \mathcal{Z}\right]\right\} \ &+\zeta \zeta_d \left[(k_x^2 + k_y^2)\tau -1 \right]\left\{1 – \zeta_* \mathcal{Z} +2\left( \zeta^2 – \frac{1}{2} \right)\left[1+\zeta \mathcal{Z} – \zeta \zeta_* – \zeta_* \left( \zeta^2 – \frac{1}{2} \right) \mathcal{Z}\right] \right\} \end{align} where $$\zeta = \frac{\omega}{2k_z}, \quad \zeta_* = \frac{k_y T}{2 k_z}, \quad \zeta_d =\frac{k_y B}{2 k_z},$$ with $$T$$, $$B$$ and $$\tau$$ simply numerical constants, while $$\mathcal{Z}$$ is the plasma dispersion function $$\mathcal{Z}(\zeta) = \frac{1}{\sqrt{\pi}} \int \text{d}x \: \frac{e^{-x^2}}{x-\zeta},$$ which, in Mathematica, can be expressed, for x$$= \zeta$$, as I*Sqrt[Pi]*Exp[-x^2]*(1+I*Erfi[x]). The presence of this function is the reason that the dispersion relation is non-polynomial, and is also badly-behaved for certain values of $$\zeta$$