VRchat now has “IK 2.0”, which is (mostly) able to avoid self-collisions with the avatar. It's not a full collision detection system, just enough to put elbows and knees in reasonable places when you have torso, wrist and foot tracking only. Without this, you get arms going through the torso and such.

Any good theory papers on this? Open source code? Thanks

Maybe what they do is just supporting joint limits? Doing this prevents arms from penetrating torso, but one knee could still penetrate the other knee.
Capsule based collision detection wouldn't be any harder, though.

Personally i support joint limits. For natural limits, there exist many conventions.
A cone limit with nonuniform scale, so we get an elliptical cross section to restrict swinging motion. (Jolt physics engine supports this, as i've recently seen.)
Or a boolean combination of multiple cones. I remember Havok had this, although i found that somewhat strange.
Personally a also using a combination of cones. Here a picture of the shoulder region:

Advantage of cones is that we can calculate a closed form analytical solution. And i'm pretty happy with the behavior. Still need to tweak the limits to be more natural.

This models the swinging limits. Beside that we need to limit twisting motion with two angles. (Visualized by yellow arcs.)
To calculate the current twist angle, i use such code:

inline float QuatRotAngleX (const sQuat &q0, const sQuat &q1)
	{
	// factor rotation about x axis: sQuat qt = q0.Inversed() * q1; float halfTwistAngle = atan (qt[0] / qt[3]);
		return 2.0f * atan (
			( ( ( (q0[3] * q1[0]) + (-q0[0] * q1[3]) ) + (-q0[1] * q1[2]) ) - (-q0[2] * q1[1]) ) /
			( ( ( (q0[3] * q1[3]) - (-q0[0] * q1[0]) ) - (-q0[1] * q1[1]) ) - (-q0[2] * q1[2]) ) );
		//sQuat qt = q0.Inversed() * q1; 
		//return 2.0f * atan (qt[0] / qt[3]);
	}
	inline float QuatRotAngleY (const sQuat &q0, const sQuat &q1)
	{
		sQuat qt = q0.Inversed() * q1; 
		return 2.0f * atan (qt[1] / qt[3]);
	}
	inline float QuatRotAngleZ (const sQuat &q0, const sQuat &q1)
	{
		sQuat qt = q0.Inversed() * q1; 
		return 2.0f * atan (qt[2] / qt[3]);
	}

I guess that's pretty standard (optimized the X one a bit).
This can still cause surprising results if the swing goes > 90 degrees (which never happens with natural limits), but better than any other method i have tried.

I have included those limits recently into my IK solver in a way so the solver finds a solution not violating limits if it exists.
It was possible without a need for iterations.
Before that, i did only ‘slide’ along the limits. Meaning it does not violate limits, but does not find the solution in one step if they are violated.
That's much simpler, not hard at all, and good enough for very most applications i would say.

Joint limits are not enough. Elbow position is an underdetermined problem if you only have torso and wrist position. Sit back and grab a table. Now move your elbow around.

Phia, of the Virtual Reality Show, tries out VRChat's IK 2.0.

Notice her putting her arms across her body and one hand behind her back. The IK gets that right and keeps the arm outside the torso, although the tracker is not providing elbow info.

Nagle said:
Joint limits are not enough. Elbow position is an underdetermined problem if you only have torso and wrist position.

I insist, because the problems you mention are the exact same which motivated me to improve the IK solver to respect joint limits.

Notice, preventing collisions may sound important because it avoids failure cases, but it does not help you to solve the actual problem, which is to get natural poses.
And also notice: The human body structure has evolved to avoid self collisions or other such problems. If we achieve natural poses, we won't run into them at all. At least not so often it remains a big problem.

That said, the remaining question is: How do joint limits help us to get natural poses?

The primary goal of IK is to drive an an effector to some target. If we implement this, we observe we deal with an underdetermined system. Usually many solutions are possible (or none at all).
And we realize the open variable is just one: The swivel angle.
If you grab a door handle, you can rotate your elbow on a circle around the line from shoulder to wrist, while still holding the door knob at desired position and rotation. That's the swivel angle.
Thus we can rephrase our problem of getting a natural pose to the question: What's a good choice for the swivel angle?

To answer this, there are multiple ideas, which i call secondary objectives for the IK solver:
Minimize energy to reach the pose (easy: just minimize the change form the initial pose)
Minimize energy to hold the pose (easy: just make the elbow hang down with gravity)
Maximize agility (hard: calculate the pose which tires to keep most distant from all joint limits)

I have never seen a IK solver in animation programs which can do the agility objective. Implementing it was not easy, but i am very happy with the results.
It finally allows me to automate the swivel angle. Before that, i always had to control it manually, at least in cases. No matter if i worked on animations or on robotic ragdolls.
Now i can set weights for those 3 objectives, and giving agility the highest weight always gives me good results in any situation.

So that's probably something nobody else could tell you. It may not be worth the effort for your application, but you can consider it.

Here's how it works:

You can see the circle for the elbow hinge joint.
This circle needs to be clipped by the limits of both the spherical shoulder and wrist joints.
For the swing limits that's intersection of cones with a circle. But the circle always has the same plane than the axis of the cones, so analytical solution is possible (but still pretty difficult).
For the twist limits i also found an analytical solution from intuition and trial and error (i do not really understand how it works, but it does).
Various visualizations of colored arc segments show those limits here. After that i do a boolean union of those arc segments to finally get the potential solution space.
The maximized agility then simply is the center of said space, which mostly is just a single arc.
But it can also be two arcs, which is a problem, as we need to choose. And our choice my jump between the two options as the pose dynamically changes.
This means that my advanced solver actually generates more flipping cases than a naive solver which only slides on the limits.
However, knowing the full solution space is the only way to calculate the agility objective respecting limits.
Besides, the advanced solver as is much harder to make smooth and practical in cases where no solution exists. It can find the closest point to a potential solution, but again - this point is more likely to jump around than with a naive, sliding solver.

As i work on robotic ragdolls, it's worth for me to deal with such complexity.
For something like VR Chat, i would say it's overkill.
But really depends on your needs, so you should know.
ML surely is an alternative and eventually easier way to solve this.