How do I quantize a classical field theory I have not been able to find any information about this on the Internet. I am a middle-schooler, 14, who self-studies physics, and I know up to and including ODEs, and some of the calculus of variations, enough so that I can take a variation of the action. However, I have absolutely zero budget, and no ability to get textbooks. So, I was wondering how a physicist would quantize a classical field theory, such as Gauss's law for gravity, which I'm currently trying to create a quantum field theory for, as part of my plan to create a possible theory of quantum gravity, by adding in the relativistic corrections later.
 A: First off, props for even knowing the proper use of the terminology. I dare say most middle-schoolers probably think classical field theory is about listening to Mozart in a meadow.
Now, if I may start at the beginning...
You seem to know this already, but I'll restate it for pedagogical reasons. When someone says they have a classical field theory they want to quantize, that means they usually have a Lagrangian density in terms of a scalar, gauge, vector, spinor, or other type of field.  Additionally, the Lagrangian usually obeys a known (that is known to them) type of symmetry (eg rotational invariance) and they usually have the equations of motion for it as well. Specific to your example, Gauss' Law for gravity is the equation of motion for the Newtonian gravity Lagrangian density (Wikipedia can provide this equation). In this Lagrangian, the scalar field $\phi$ is the important component, not a gravity field.
Now, assuming a physicist has a Classical field theory, one of the standard methods of quantizing it (promoting it to a Quantum Field Theory) is this:
1. Solve for a general case solution to $\phi$ if possible, this greatly eases the process.
2. Write the operator $\phi$ in terms of annihilation and creation operators (let's call them $a_k$ and $a_k^\dagger$ respectively).
As an example, if I managed to find a general solution for my $\phi$ was:
$$\phi(x)=e^{ikx}+e^{-ikx}$$
Then I could rewrite it as:
$$\phi(x)=A(a_ke^{-ikx}+a_k^\dagger e^{ikx})$$
3. Impose the cannonical (anti-)commutation relations on the theory. If the theory is written using bosons (and I believe yours is), then use commutators. If it uses fermions, then use anti-commutators.  For your case, the cannonical commutation relations should be:
$$[\phi_a(\mathbf x),\Pi_b^0(\mathbf y)]~=~i\delta_{ab}\delta^3(\vec x-\vec y)$$
$$[a_{k},a_{k'}^\dagger]~=~\delta^3(\vec k-\vec k')$$
When I say "impose these relations" what I mean is that in step 2, you need to put normalization coefficients (coefficients like "A" simply used to get the right outcome in step 3) in front of the creation and annihilation operator terms. Then, when imposing the relations, you need to set LHS=RHS and solve for these coefficients such that the relations hold. It is important to note that these coefficients need not be constant, just that they be c-numbers (not functions of x nor operators).  In fact, one of the most basic coefficients ends up with a $\omega_k$ in it, which is a function of the momentum, k.
Then you're done! If you can do steps 2 and 3, then your resulting theory is a quantum field theory.  I'm sure others will let me know I've missed a method or two. If you think of something I've forgotten, let me know and I'll add it in.
Edit
It has been pointed out that I (somehow) forgot to mention the path integral formulation method of quantizing a theory. This is a conceptually harder method that requires a more advanced understanding (i.e. formal training) of theoretical physics. I am not including a full description of it here because I don't want to increase the level of complexity of this answer too much. In general, path integral is found to work better than the one I listed, but it is certainly not as straight forward (debatable).  It also is inherently relativistic, which is an advantage over the more basic method (unless you intend to add in relativity later).
A: "I am a middle-schooler, 14, who self-studies physics"
Yeah, I was there, doing that kind of stuff, too at around that age. You're not 14 anymore, but 22 or 23 at the time of writing. So, now we have better perspective and we can address this matter more cogently.
"I'm currently trying to create a quantum field theory for, as part of my plan to create a possible theory of quantum gravity, by adding in the relativistic corrections later"
You need to worry more about how you actually do that second step, there, not the first one! There is no known non-relativistic version of General Relativity cast in a form of a Lagrangian field theory whose Lagrangian is the non-relativistic limit of the Einstein-Hilbert or Einstein-Cartan Lagrangians or any of the others used for General Relativity or its variants. And it may very well be impossible for there to be any.
Moreover, quantizing the non-relativistic theory of gravity as a hybrid field (for gravity) and mechanics (for particles) is going to crash, too. You still have the divergence problem - every particle experiences its own gravity potential as a self-force. And you still have the same issues with "infinity subtraction" to try and devise some consistent non-ad-hoc way to account for it. That's the very problem you have for any quantum field theory too, never mind a quantum theory of gravity. So, you're not getting away from anything. It'll hound you and come after you, even all the way back to Newtonian gravity.
You could cheat and treat the combined hybrid system as just pure mechanics, with the "potential" reduced to particle coordinates, but then you'd be missing the whole point of the exercise ... unless you somehow find a way to replace field theories by particle-based mechanics theories in Relativity (and make that the thing you 'relativistic-correct' to); but pure particle dynamics in Relativity are also hounded by various road-blocks, like the Leutwyler Theorem.
As for field theories, there is Newton-Cartan gravity, but that is formulated in terms of field equations - one that's very shoddily cobbled together - rather than in terms of any action principle. And it's incomplete: there is no field law or dynamics for the mass density, itself.
A field theory is quantized, only when it is first formulated in terms of an action principle. Attempts at coming up with a Lagrangian field theory of Newtonian gravity, as both a theory for gravity and fluid dynamics, have never been done. Remember: the Einstein equations for General Relativity include both the law of gravity and fluid dynamics in them - effectively being a formulation of fluid dynamics, itself, in the presence of gravity, with the gravity being represented chrono-geometrically as a warping of the space-time continuum (with most of the effect we feel as gravity being a warping of the time dimension).
When you lay out Newton-Cartan and General Relativity side-by-side and try to match the parts, they don't match up at all; and there's no fluid dynamics in Newton-Cartan at all.
The main reason for the discrepancy is that there is an even deeper discrepancy between relativistic and non-relativistic physics (in the form they are currently formulated as) which seems to block any prospect of the "start with non-relativistic theory, add in relativistic corrections later". The precise term for that last step is called "deformation", and the venue for addressing this is what's known as deformation theory.
It is quite likely that one or both theories have not been formulated quite right, as they miss each other in the forwards direction as a "deformation" or backwards direction as a "non-relativistic limit". Something more needs to be added to each; which, in the case of non-Relativistic physics, means "retro-added"; so that they line up right.
Non-relativistic physics is purportedly governed by the symmetry group known as the Galilean Group, while relativistic physics is governed by the Poincaré group and its homogeneous version, the Lorentz group. General Relativity actually incorporates neither in its formulation, per se; e.g. some solutions to the Einstein Field equations may permit the space-time signature to change, so that a locally Minkowski space-time in one region is connected across a singular surface to a locally Euclidean 4-dimensional timeless space in another - the "signature changing" solutions.
Historically, General Relativity was formulated on the Einstein-Hilbert action, which lives on a pseudo-Riemannian geometry, but is actually agnostic on what the signature of the underlying geometry is.
It can be formulated on a Riemann-Cartan geometry, which directly embodies a localized version of the Lorentz group into it. That leads to either Einstein-Cartan theory, if using the Einstein-Cartan action, or General Relativity itself reformulated in Riemann-Cartan geometries by using the Palatini action in place of the Einstein-Hilbert action. So, either of these, or something similar or upwardly-compatible are what you should be aiming for by "relativistic correction to (...)".
The main discrepancy is that for non-relativistic theory, the inhomogeneous Galilei group is not the non-relativistic limit of the Poincaré group, at all; though the homogeneous Galilei group is the limit of the Lorentz group, and the Lorentz group can be presented as a deformation of the homogeneous Galilei group.
Part of this is because the correct symmetry group for non-relativistic theory is neither the homogeneous nor inhomogeneous Galilei group, but the Bargmann group - which is the "central extension" of the inhomogeneous Galilei group. It does not have the Poincaré group as a deformation, but a trivial central extension of Poincaré, as one.
The main punchline is that while the geometry naturally associated with the Lorentz and Poincaré groups is the 3+1 dimensional Minkowski geometry, the geometry associated with Bargmann is 4+1 dimensional (the Bargmann geometry). You have a total mis-match here.
The most direct way to see this is that while the 3-components of momentum and of the total energy - or 4 components in all - transform together in Relativity as a 4-vector under the Lorentz or Poincaré, in non-relativistic theory, you don't have any such thing as "total energy", but instead you only have kinetic energy and mass. That's 5 components, in all, not 4. They transform together, under Bargmann as a 5-vector.
You can't deform this to any "relativistic" version, except by treating kinetic energy and mass, once again, as independent components, even though they continue to have non-trivial transformations into each other under Lorentz or Poincaré ... but now as components of a 5-vector.
All of this is central to the whole enterprise of defining a "stress tensor". If you line up the components in the stress tensor in a matrix, the rows give you the components for the continuum-version of each of the above items. So, in Relativity, the top row would be the 4-components for a "Energy current", and each of the next three rows would be the 4-components for each of the components of the "momentum current". That's a 4 x 4 matrix in all.
For non-relativistic theory, you'd have 5 rows. In order to make the matrix symmetric - as the stress tensor should be - it has to also have 5 columns, which means 5 coordinates. That's a 5 x 5 stress tensor matrix. The top row is for mass density, and gives you the fluid dynamics balance equation for it. The next three rows are for each of the three components of momentum, giving you the three balance equations for momentum density. Finally, there is the bottom row that gives you the density for kinetic energy.
It is quite possible to take the 4+1 dimensional version of the Einstein-Hilbert action in a locally Bargmann geometry, or its Riemann-Cartan geometric versions (a 4+1 dimensional Palatini or Einstein-Cartan action). Additional constraints then need to be imposed to narrow the results down to Newtonian gravity, while retaining some ability to deform that into a Relativistic version for General Relativity; but now with its space-time geometry embedded in 4+1 dimensions. Nonetheless, even with all that, there still continues to be a mismatch.
The closest I've seen to it being accomplished is "Bargmann Structures and Newton-Cartan theory", by Duval, Burdet, Künzle and Perrin in Physical Review D 31(8), 1985 April 14. It's not open access, you have to go to the library for it or log onto a machine on a campus that carries the journal.
To quote them, at the head of their section "VI. Newton's Field Equations": "Strangely enough, Newton's field equations [...] cannot be easily derived from a specific space-time Lagrangian density. It just seems that there might exist some purely geometric obstruction to the existence of a well-defined variational problem in the four-dimensional picture."
The reformulation on a locally Bargmann geometry to follow helps, or as they say "We will show here that the introduction of Bargmann structures improves the situation", but it does not entirely resolve the problem. There's still something more that's missing; and I think that extra missing item is on the Relativity side of the mis-match, not the Newtonian side. That is: we might not actually have the full version of what the correct relativistic theory should be, but only a reduced and constrained version of what it ought to be.
