Real and Imaginary Frequency Responses of a single complex pole
If I understood you correctly you want a formula for the magnitude of the frequency response of a filter with a single complex pole. As you've already noticed, the system function is

$$H(z)=\frac{1}{1-az^{-1}}\tag{1}$$

which is valid regardless of whether $a$ is real-valued or complex-valued. Let's assume that $a=a_r+ja_i$ is a complex-valued number. If we also assume that $|a|&lt;1$ (i.e. the system is causal and stable), then we get the systems frequency response by replacing $z$ by $e^{j\omega}$ in (1):

$$H(e^{j\omega})=\frac{1}{1-ae^{-j\omega}}\tag{2}$$

Now you just have to compute the magnitude of (2). Using

$$|1-u|^2=1-2\text{Re}\{u\}+|u|^2$$

for any complex $u$, we get with $u=ae^{-j\omega}$ from (2)

$$|H(e^{j\omega})|=\frac{1}{\sqrt{1-2\text{Re}\{ae^{-j\omega}\}+|a|^2}}=
\frac{1}{\sqrt{1-2(a_r\cos\omega+a_i\sin\omega)+|a|^2}}\tag{3}$$

EDIT:

This is how you can derive the transfer functions relating the real-valued input to the real part and the imaginary part of the output, respectively. Let $y(n)=y_r(n)+jy_i(n)$. Then you have the following two difference equations for the two real-valued output signals:

$$y_r(n)=x(n)+a_ry_r(n-1)-a_iy_i(n-1)\\
y_i(n)=a_ry_i(n-1)+a_iy_r(n-1)$$

Transforming these equations to the $\mathcal{Z}$-domain gives

$$Y_r(z)=X(z)+a_rz^{-1}Y_r(z)-a_iz^{-1}Y_i(z)\\
Y_i(z)=a_rz^{-1}Y_i(z)+a_iz^{-1}Y_r(z)$$

Now you have two equations for the two unknowns $Y_r(z)$ and $Y_i(z)$. After some algebra you arrive at

$$Y_r(z)=H_r(z)X(z)=\frac{1-a_rz^{-1}}{1-2a_rz^{-1}+(a_r^2+a_i^2)z^{-2}}X(z)\\
Y_i(z)=H_i(z)X(z)=\frac{a_iz^{-1}}{1-2a_rz^{-1}+(a_r^2+a_i^2)z^{-2}}X(z)\tag{4}$$

Note that even though

$$H(e^{j\omega})=H_r(e^{j\omega})+jH_i(e^{j\omega})$$

$H_r(e^{j\omega})$ and $H_i(e^{j\omega})$ are not the real and imaginary parts of the total frequency response $H(e^{j\omega})$ of the complex-valued system because they are both themselves complex-valued. You can obtain the magnitudes of the frequency responses by computing the magnitudes of $H_r(e^{j\omega})$ and $H_i(e^{j\omega})$ as given in (4), i.e. by replacing $z$ by $e^{j\omega}$ and computing the magnitudes of the two complex functions.

Another option is to express the frequency responses $H_r(e^{j\omega})$ and $H_i(e^{j\omega})$ using the total frequency response $H(e^{j\omega})$ (2):

$$H_r(e^{j\omega})=\frac12[H(e^{j\omega})+H^*(e^{-j\omega})]\\
H_i(e^{j\omega})=\frac{1}{2j}[H(e^{j\omega})-H^*(e^{-j\omega})]\tag{5}$$

From (5), the magnitudes can be written as

$$|H_r(e^{j\omega})|=\frac{1}{2}\sqrt{|H(e^{j\omega})|^2+2\text{Re}\{H(e^{j\omega})H(e^{-j\omega})\}+|H(e^{-j\omega})|^2\}}\\
|H_i(e^{j\omega})|=\frac{1}{2}\sqrt{|H(e^{j\omega})|^2-2\text{Re}\{H(e^{j\omega})H(e^{-j\omega})\}+|H(e^{-j\omega})|^2\}}$$

where $|H(e^{j\omega})|$ is given by (3).

EDIT 2:

Here's some Matlab/Octave code using equation (4) to plot the different transfer functions. $|H_r(e^{j\omega})|$ and $|H_i(e^{j\omega})|$ are of course symmetrical because they correspond to real-valued systems, whereas $|H(e^{j\omega})|$ is not due to the complex pole.

w = linspace(0,2*pi,300);          % frequency vector
ar = .4; ai = -.3; a = ar+j*ai;    % complex pole
u = exp(-j*w);                     % u = z^(-1)
Hr = (1-ar*u)./(1-2*ar*u+(ar^2+ai^2)*u.^2);
Hi = (ai*u)./(1-2*ar*u+(ar^2+ai^2)*u.^2);
H = 1./(1-a*u);
plot(w/pi,abs(Hr),w/pi,abs(Hi),w/pi,abs(H))
legend('Hr', 'Hi', 'H');




EDIT 3:

If you need to split the numerators and denominators of (4) into real and imaginary parts, here's what you get:

$$H_r(e^{j\omega})=\frac{1-a_r\cos\omega+ja_r\sin\omega}{\left[1-2a_r\cos\omega+(a_r^2+a_i^2)\cos(2\omega)\right]+j\left[2a_r\sin\omega-(a_r^2+a_i^2)\sin(2\omega)\right]}\\
H_i(e^{j\omega})=\frac{a_i\cos\omega-ja_i\sin\omega}{\left[1-2a_r\cos\omega+(a_r^2+a_i^2)\cos(2\omega)\right]+j\left[2a_r\sin\omega-(a_r^2+a_i^2)\sin(2\omega)\right]}$$