I’m curious if there are recommended settings for the degree and order of spherical harmonics? When I did some testing with degree and order set to 30 we had residuals of roughly 6 meters, and when degree and order were set to 100 our residuals were over 100 meters. I thought this was strange since it seems like the larger the degree and order, the more accurate the fit should be. So I’m curious if there is a recommended max value for the degree and order.

The maximum value depends on both the accuracy you need and the altitude of the orbit.
For LEO satellites (or satellites with highly eccentric orbits and low perigee), I have seen people use degree/order 60 for classical missions, and 120 for high accuracy missions.
For MEO satellites (like GNSS), order 12 is often sufficient, but you may want to go slightly beyond that.
For GEO satellites, order 8 is also often sufficient.

Accuracy should improve if you use higher degrees. If yo do not see that, then there is a problem in the modelling. A classical error is to use a too high tolerance for the integrator. We made this mistake in many of our tutorials (we should correct them) by using a tolerance of 10m, but now we recommend a tolerance of 0.001m.

Do your residuals come from real measurements or simulated measurements? In the later case, you should check the settings of the generator too.

My integrator is set here to be:
integratorBuilder = DormandPrince853IntegratorBuilder(prop_min_step, prop_max_step, prop_position_error)

I checked and our prop_position_error was set to 10 meters. We changed it to 0.001 meters and this seemed to decrease our residuals. Is this the value you mentioned for the tolerance for the integrator? And yes, this comes from real measurements.

Would you recommend always setting the tolerance (which I assume is what I called prop_position_error) to be 0.001 meters, or is this something that we should change based on our measurements?

Yes, this is the value I mentioned for the tolerance.

You should keep 0.001 meters even if your measurements are less accurate than that. This setting is an internal one used for propagation. It is only an estimated error and it is only for one step, the global error will be far greater that that because errors have a tendency to build up during propagation. This is the reason why we now recommend a very low tolerance so the propagation remains realistic.