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$