How are composite hadron fields related to elementary quark fields? (This question is related to: A pedagogical exposition of the hadron physics?)
I'm a mathematician who has been trying to learn quantum field theory for a while. I've gone through large parts of quite a few books, and there's a conceptual issue that's been bothering me for a while.
It's common to see statements like "the proton is composed of two up quarks and a down quark," or even equations like $\pi^0=\frac{1}{\sqrt{2}}(u\bar u-d\bar d)$. The reading I've done suggests that one builds an effective field theory containing fields corresponding to all the hadrons, but that the exact relationship between the hadron fields and the quark fields is just poorly understood, so, sad as it might be, there's no way to say what the hadron fields have to do with the quark fields.
But all the talk about hadrons being "made of" particular combinations of quarks and antiquarks makes it seem like this can't be the whole story. In particular, one finds these composite particles by looking for copies of the trivial representation in a tensor product of $SU(3)$ representations, so it seems like a "proton field" should somehow be related to some corresponding product of the up and down quark fields. How does this work? When someone says that a proton is made of two up quarks and a down quark, what are they saying about the relationship between the proton field and the up and down quark fields?
 A: Here is my basic understanding of field theory, applied to quantum mechanical systems, from nuclei to elementary particles.
For particle physics usage, one starts with the particles and the quantum mechanical equation which these particles fit as free waves : Klein Gordon for bosons,Dirac for fermions and quantized Maxwell for photons . One then has the basis to build a field theory: each particle has  a  representative field in all space  time, on which quantum mechanical creation and annihilation operators act, creating a particle at (x,y,z,t) or annihilating a particle  and thus generating the propagation of particles in space time. The use is that they allow many particle interactions and the Feynman diagrams predict measurabe quantities for these many body interactions.
Thus one always has to start with the ground state of "particle" in particle physics.
In your example, you are building a field theory where the particle is the pi0, a boson, and has its free particle wavefunction representing the pi0 field. It has no meaning to ask :

When someone says that a proton is made of two up quarks and a down quark, what are they saying about the relationship between the proton field and the up and down quark fields?

because if you are building a field theory with protons, the proton field is just the plane wave wave function solution of the dirac equation with the mass of the proton. Any connection with its compositenes is not reflected in this plane wave.
Maybe as a mathematician, you think that all these fields are really there in space time? They are just as a complicated coordinate system, they are not sitting there in space so as to have to have definable a relationship.
If you want to describe a pi proton scattering using a field theory that only has hadrons  as a basis, you can compare it with the data , and will find  that the comparison fails at high energies ( quark and gluon jets) and a different field theoretical model is needed, based on the composite state of protons and pions, the fields have to be the plane wave solutions of quarks and antiquarks etc. The new model will describe the data where also the hadron model did, so it dominates, but due to the different wavefunction basis for the two field theories, to demonstrate how the simpler hadron  one emerges is not analytically possible.  
A: 
When someone says that a proton is made of two up quarks and a down quark, what are they saying about the relationship between the proton field and the up and down quark fields?

This won't be anything close to a complete answer, but it highlights one context in which this question has been addressed relatively directly: numerical calculations using lattice QCD. 
One of the major goals of these numerical calculations is to predict the spectrum of hadrons. One way to do this is to consider a vacuum expectation value like $G(x-y)=\langle 0|A(x)A(y)|0\rangle$ where $A$ is, say, the effective field operator for some hadron. After evaluating $G(x-y)$ numerically as a function of $|x-y|$, we could infer the hadron's mass $m$ by fitting the result to $G(x-y)\sim \exp(-m|x-y|)$. I'm glossing over lots of details, of course; but the main point here is that this only works if we have some idea of which operator we should use for $A$. Since the model is formulated in terms of quark and gluon fields, this means that we need to have some idea of how to express a hadron field operator $A$ in terms of the quark and gluon fields. 
We don't know exactly how to do that, but for some applications we know enough. Things like the hadron's spin, charge, and other conserved quantities (like isospin in pure QCD) provide some constraints, and as usual, in the absence of additional information, the simplest choice is a good place to start. This is addressed in section 5.2 in [1], which says this on page 261:

In order to project out the channels with different quantum numbers the appropriate hadronic operators have to be constructed from the quark and gluon fields. The choice of the composite operators is to a large extent arbitrary. In fact, ...one has to find the optimal operator, which has a strong enough copuling to the hadron in question and, at the same time, can be evaluated without too great difficulties. 

The word interpolating operator is sometimes used to refer to a choice of the operator $A$. The book goes on to show examples like $A\sim \overline{d}\gamma_5 u$ for the charged pion and $A\sim$[a particular combination of $uud$] for a proton. The idea is that when one of these operators $A$ is applied to the vacuum state $|0\rangle$, the resulting state-vector will be a superposition in which at least have one term corresponds to a single-particle state of the hadron of interest. Choosing $A$ thoughtfully can help minimize the contributions of other terms, or at least help isolate the hadron with the smallest mass, which is the one that dominates $G(x-y)$ at large $|x-y|$. This is acknowledged on page 266 in [1], which says:

The quark-gluon composite operators for the calculation of hadron masses in lattice QCD simulations have to be chosen carefully, in order to minimize the errors of the results. Besides the quantum number structure discussed in the previous subsection, the other ingredient is the coordinate dependence of the trial wave functions for mesons and baryons. For a strong overlap resulting in a high signal to noise ratio, the quark-gluon distributions in space have to bear some qualitative resemblance to the true wave funtions.

The paper [2] includes a relatively concise overview of these ideas. As noted in the OP and illustrated by this answer, the precise relationship between hadrons and the quark/gluon fields is still poorly understood. No short list of references can fairly represent all of the thought that has gone into this, but a few other examples include [3], [4], and [5].
One of the most interesting insights comes from the large-$N_c$ limit, where $N_c$ is the number of colors. This limit might seem far removed from reality (where $N_c=3$), but numerical calculations suggest that in some ways it is a surprisingly good approximation. Assuming that QCD is still confining for large $N_c$, the analysis published in [6] concludes that mesons are pure quark-antiquark in the large-$N_c$ limit (that is, the valence-quark approximation becomes exact), along with a number of other interesting conclusions. Regarding baryons in the large-$N_c$ limit, section 38.7 in [7] says this:

Large-N QCD can be treated as a weakly coupled field theory of mesons. It is a theory of effective local meson fields, with effective local interactions, in which the three-meson coupling scales as $1/\sqrt{N}$, the four-meson as $1/N$, and so on. At large $N$ all coupling constants are weak. As we know already, many weakly coupled field theories possess, in addition to elementary excitations, heavy solitonic states whose masses diverge at weak coupling as the inverse of the coupling. Are there such states in QCD and its effective mesonic counterpart? The answer is positive. In QCD we have $N$-quark states — baryons — whose mass is proportional to $N$. As a reflection of this fact, the low-energy mesonic theory must have solitons with nonvanishing baryon numbers and masses scaling as $N$. These are the Skyrmions ... some implications of the Skyrmion model are model-independent; they follow from QCD in the ’t Hooft limit...

The point of this excerpt is that whereas QCD describes baryons in terms of quark and gluon fields (in principle), the low-energy effective theory describes baryons in a completely different way (namely as solitons). This illustrates just how non-trivial the relationships between the assumed fields and the predicted particles can be in QFT.

References:
[1] Montvay and Münster (1994), Quantum Fields on a Lattice (Cambridge Univeristy Press)
[2] "Group-theoretical construction of extended baryon operators in lattice QCD," https://arxiv.org/abs/hep-lat/0506029
[3] Selem and Wilczek (2006), "Hadron Systematics and Emergent Diquarks," https://arxiv.org/abs/hep-ph/0602128
[4] Lorcé and Liu (2016), "Quark and gluon orbital angular momentum: Where are we?," https://arxiv.org/abs/1601.05282
[5] Greensite (2011), An Introduction to the Confinement Problem (Springer)
[6] Witten (1979), "Baryons in the $1/N$ expansion," Nuclear Physics B 160: 57-115
[7] Shifman (2012), Advanced Topics in Quantum Field Theory: A Lecture Course (Cambridge University Press)
