I am busy working on a 3D $n$-body galaxy simulator and I am having some difficulty wrapping my head around determining the initial velocity components of each particle to ensure a circular orbit around a central mass.
I know that the magnitude of the velocity to ensure a circular orbit is calculated from: $$ v = \sqrt{\frac{GM}{r}} $$
How then can I determine what the individual velocity components in 3D will be for the particle?
If the bodies are noninteracting (and only the central mass is providing the interaction needed for the bodies to rotate) you can simply do the following.
For each body, choose a distance from the mass $M$, call it $r$: from this you can calculate the absolute velocity, as you did. Place the object at any point that is distant $r$ from the origin.
Now, the motion will happen in the plane that contains the initial position of the object, its initial velocity and the center of rotation. This is essentially a 2D problem. In this case, have a look at the expression for velocity and acceleration in polar coordinates: $$ v=\dot r\hat e_r+r\dot \theta\hat e_\theta,\\ a=(\ddot r-r\dot \theta^2)\hat e_r+(r\ddot\theta+2\dot r\dot\theta)\hat e_\theta, $$ where $\hat e_r$ and $\hat e_\theta$ are polar directions in the plane of motion. The acceleration is determined by gravitational force, that is purely radial. If you try a test motion where $r(t)=r$ (the one you choosed) and $\theta(t)=\frac{v}{r}t$ (the $v$ you computed from $r$), you see that this motion has the exact acceleration to solve Newton's problem and is compatible with boundary conditions on the position. You just have to set your initial conditions on the velocity, by saying that at time $0$ there is no radial velocity, and the velocity is purely radial: motion equations will ensure that you get a purely circular motion, with $\dot r=0$.
Returning to 3D space, it is straightforward to generalize our analysis: placing a point at position $r$, its initial velocity in spherical coordinates must be of the magnitude that you calculated, and must have no radial component. This will ensure a circular motion in the initial velocity/center of attraction plane. You can return to Cartesian coordinates by applying elementary transformation laws.
The velocity in $(x,y)$ plane at the circular orbit can be viewed best in polar coordinates (spherical for $\theta \to \pi/2$ if you have 3D model). Then the initial velocity must have $v_r = 0$ and $\vec v = v_\phi \vec e_\phi$. You can calculate the $v_\phi$ through \begin{equation} v = \sqrt{\vec v \cdot \vec v} = \sqrt{v_r^2 + r^2 v_\phi^2} = r v_\phi. \end{equation}
If you want to work in Cartesian coordinates (which I assume are better for you), you have \begin{equation} \renewcommand{\vec}[1]{\boldsymbol #1} \vec v = \frac{\textrm d \vec r}{\textrm d t} = (-r\sin(\phi)\frac{\textrm d \phi}{\textrm d t},r \cos(\phi)\frac{\textrm d \phi}{\textrm d t}) = \Omega r (-\sin(\phi),\cos(\phi)), \end{equation} where $\vec s = (-\sin \phi, \cos \phi)$ is the direction you are looking for. $v = \Omega r$ is then general case of what you have written for gravitation field, so you can use your formula for $v$ and use this vector $\vec s$ for your direction. If you have net in Cartesian coordinates, you have to realize that \begin{equation} \tan\left(\phi\right) = \frac y x, \end{equation} so if you start in place $(x,y)$, your direction is \begin{equation} \vec s = \left( -\sin\left( \arctan\left(\frac yx \right) \right),\cos\left( \arctan\left(\frac yx \right) \right) \right). \end{equation}
Edit: In this derivation I assumed radial velocity is equal to zero, that's why the vector is so simple. And this assumption is why it works for your circular orbits.
There is possibility to use spherical coordinates instead and get velocity vector \begin{align} \vec v = r\left( \begin{array}{c} v_\theta \cos (\theta) \cos (\phi)- v_\phi \sin (\theta) \sin (\phi) \\ v_\theta \cos (\theta) \sin (\phi)+ v_\phi \sin (\theta) \cos (\phi) \\ - v_\theta \sin (\theta) \\ \end{array} \right) \end{align} and remember that \begin{align} \phi = \arctan\left(\frac yx\right) \qquad \text{and} \qquad \theta = \arccos\left( \frac{z}{x^2 + y^2 + z^2} \right). \end{align}
The last step is to calculate the velocities $v_\theta, v_\phi$. These are equal to angular velocities $\Omega_\theta, \Omega_\phi$, so you has to project your velocity into these two components. I think the best here is to realize, that \begin{align} v = \sqrt{\vec v \cdot \vec v} = r^2 \Omega_\theta^2 + r^2 \sin^2(\theta) \Omega_\phi^2 = \sqrt{\frac{GM}{r}}, \end{align} then you can define $G,M$ and $x,y,z$. By postulating these constants you have well defined $r,\theta$, so the previous equation has only two unknowns $\Omega_\theta, \Omega_\phi$. Now you have one degree of freedom, you can chose one of these numbers and the second one must satisfy the equation, for exapmle you can set $v_\phi = \Omega_\phi$ and get \begin{align} v_\theta = \Omega_\theta = \sqrt{\frac{GM}{r^5}} - r^2\sin^2(\theta)\Omega_\phi = \sqrt{\frac{GM}{r^5}} - r^2\sin^2(\theta)v_\phi. \end{align}
Summary: If you are in 3D, you can chose the rotating direction by choice of $v_\theta, v_\phi$ (that's up to you) and than you have the initial velocity given as \begin{align} \vec v = r\left( \begin{array}{c} v_\theta \cos (\theta) \cos (\phi)- v_\phi \sin (\theta) \sin (\phi) \\ v_\theta \cos (\theta) \sin (\phi)+ v_\phi \sin (\theta) \cos (\phi) \\ - v_\theta \sin (\theta) \\ \end{array} \right) \end{align}
The underlying physics says, that you have rotation in 3D around a fixed point on the circle. That you have to choice "how much it rotates in $\phi$-direction = $(x,y)$ plane, and the program has all information now".
