Do gravitational waves radiate equally in all directions?
If so, is this an inherent property of all gravitational waves or is it just due to how they are normally produced? In other words, in a theoretical scenario, could we have gravitational waves that don't radiate isotropically?
If they do not, in a theoretical scenario, could gravitational waves that radiate isotropically exist?
No, gravitational waves are not emitted isotropically. In the weak-field limit (i.e. far from the sources), the radiation emitted by a gravitational system is determined by the third time derivative of its quadrupole moment, which, being a tensor, needs to be projected along the line of sight to yield a (scalar) energy flux. This projection is what gives the angular dependence to the energy flux.
That you need to consider the quadrupole moment is easy to understand: the dipole moment of a mass distribution can always be set to exactly $0$ if the origin of coordinates is chosen in the system center of mass. The monopole moment is the total system mass, which is conserved (in the weak-field limit), hence the monopole field at large distances is always $-GM/r$, a time constant, hence no radiation. The first physically relevant term is thus the quadrupole moment, $D_{ab}$.
We can use the exact expression to make the above more precise. Calling $X_{ab}$ the third time derivative of $D_{ab}$, the energy flux averaged over one wave period in an element of solid angle $d\!\Omega$ in the direction $\vec n$, where $\vec n$ is a unit vector, is (Landau & Lifshitz, Field Theory, vol. 2, Ch.13) is:
$$d\!I = \frac{G}{144\pi c^5}\left((X_{ab}n_an_b)^2 + 2X_{ab}X_{ab} -4 X_{ab}X_{ac}n_b n_c \right) d\!\Omega $$
Obviously, $n_a$ are the components of $\vec n$. The tensor $D$ is defined as
$$D_{ab} \equiv \int (3x_a x_b -r^2 \delta_{ab}) \rho d\!V$$
which shows $D$ to be symmetric and traceless, $D_{aa} = 0$. Both properties carry over to $X \equiv d^3 D/dt^3$, so that we can always choose a system of axes such that $X$ is (at least instantaneously) diagonal, and traceless:
$$ X =\left(\begin{matrix}X_1 & 0 & 0 \\ 0 & X_2 & 0 \\ 0 & 0 & -(X_1+X_2) \end{matrix}\right) $$
The most symmetric that $X$ can be is by having $X_1 =X_2$, but there is no way that it can have $X_3=X_1=X_2$, because of its being traceless. If we now choose $\vec n$ is the direction identified by $X_1$, we have
$$d\!I = \frac{G}{144\pi c^5}\left(4X_2^2+4X_1X_2 \right) d\!\Omega $$
while if we take $\vec n$ in the direction of $X_3$ we find:
$$d\!I = \frac{G}{144\pi c^5}\left( 2X_1^2 + 2 X_2^2 \right) d\!\Omega $$
which differs from the previous formula whether $X_1 = X_2$ or not.
