Can we do path integrals in gauge theories without fixing a gauge? I am aware that when quantizing gauge theories with a path integral, one needs to add a gauge fixing term to avoid over-counting gauge related field configurations. From an aesthetic perspective, I find this procedure distasteful. I would like to know if there is any proposal to circumvent adding this term in the Lagrangian, and to be able to do the path integral without fixing a gauge.
 A: As of today, nobody knows how to canonically quantise a classical theory with gauge symmetries. The standard approach (Dirac's algorithm) where one replaces the canonical brackets by (anti)commutators is meaningless if the symplectic form is degenerate. See Quantization of Gauge Systems, by Marc Henneaux & Claudio Teitelboim for a full discussion of this. In practice, in order to formulate a consistent theory in the canonical formalism one must first eliminate the gauge symmetries, either by turning them into (second class) constraints or by more elaborate methods.
A second, more direct approach is to follow Feynman's quantisation, where we postulate that the matrix elements can be calculated from a functional integral,
$$
A\sim\int a(\varphi)\ \mathrm e^{iS[\varphi]}\ \mathrm d\varphi
$$
Attempts for formalise the integral above in as much generality as needed have failed. A possible approach, to discretise the space of field configurations, has two possible outcomes: the lattice formulation either breaks gauge invariance (in which case we have essentially fixed the gauge by means of the regularisation), or it doesn't (in which case the integral diverges, inasmuch as we are integrating over $\mathbb R^n$ a function  that does not decay in some directions). In either case, we see that a naïve implementation of Feynman's approach cannot work either.
Even in the most pragmatical sense, the quantum theory is ill-defined in the presence of gauge symmetries: if we convene to sidestep all the formal manipulations and define the theory  through its Feynman rules (formally speaking, through Hori's formula),
$$
Z[J]\sim \mathrm e^{iS_\mathrm{int}[\delta]}\mathrm e^{-\frac i2 J\cdot \Delta\cdot J}
$$
where $\Delta$ is the inverse of the quadratic part of the Lagrangian, the programme fails, because
$$
\mathcal L_0\equiv\frac 14 F^2
$$
is not invertible.
None of these approaches seems to work. The problem can be traced back to the representations of the Poincaré Group. One may show, using the properties of the Poincaré group but nothing about Lagrangians or path integrals, that the propagator of an arbitrary vector field is
$$
\Delta(p)=\frac{-1+pp^t/m^2}{p^2-m^2}-\frac{pp^t/m^2}{p^2-\xi m^2}
$$
where $m$ is the mass of the spin $j=1$ particles created by the vector field, and $\xi\equiv m^2/m_L^2$, where $m_L$ is the mass of the spin $j=0$ particles created by the vector field.
It's easy to check that the limits $\xi\to\infty$ and $m\to 0$ are both separately well-defined, but you cannot take both limits at the same time. This means that you cannot have, at the same time, a vector field that creates massless spin $j=1$ particles and no longitudinal states. So you must either


*

*use massive particles, as in the Proca Lagrangian,

*accept that there can be negative norm states, as in $R_\xi$ QED,

*or that the field that creates particles is not a vector, as in QED in the Coulomb gauge.


In the first case the term $\frac 12 m^2 A^2$, and in the second case the term $\frac 12\xi^{-1}(\partial\cdot A)^2$, breaks the gauge invariance of the Lagrangian. In the third case the gauge is fixed by a constraint. In neither of these cases is the Lagrangian gauge invariant.
A: In lattice gauge theory, on a finite lattice, the volume $vol(\mathcal{G})$ of the group of group transformations is finite, since $\mathcal{G}$ is a finite product of copies of the gauge group $G$.   The integral 
$\int_{\mathcal{F}} \mathcal{O}(\phi) e^{-S(\phi)} d\phi$ over the space of lattice connections is also finite.   Consequently, one can compute expectation values without doing any gauge-fixing, just by computing
$$
\frac{1}{vol(\mathcal{G})} \int_{\mathcal{F}} \mathcal{O}(\phi) e^{-S(\phi)} d\phi
$$
which is equal to $\langle \mathcal{O} \rangle = \int_{\mathcal{F}/\mathcal{G}} \mathcal{O}(\phi) e^{-S(\phi)} d\phi$, as long as the observable $\mathcal{O}$ is gauge-invariant.  
Gauge-fixing is computationally convenient, especially for matching up with short distance perturbation theory, but not really necessary.  
A: You are misunderstanding what a gauge theory is if you think we shouldn't get rid of the gauge symmetry at some point. A gauge symmetry is not like other symmetries, it does not relate configurations of the dynamical variables that are physically distinct - instead, it relates configuration of the dynamical variables which are physically indistinguishable. There is no detectable difference between any configuration and its gauge-transformed version at all. Unlike, say, a rotational symmetry where a vector pointing in one direction is distinct from its rotated version, in this case, there is really no physically meaningful distinction between configurations related by gauge symmetries. See also, for instance, this question, this question, this question and more. 
Gauge symmetries reflect redundancy in the variables we have chosen to describe the system, they are entirely features of a particular theoretical choice and not inherent properties of the physical system under consideration, like e.g. rotational symmetry. There is therefore no need to try to preserve this symmetry - if it gets lost in an equivalent but more convenient description of the system, we shouldn't hesitate. It is a curious fact that rather often the gauge theoretical description turns out to be the most convenient.
Except, of course, when we want to do things like the path integral. To take the naive path integral over an action with gauge symmetry that has not been fixed is manifestly absurd physically: You are integrating over a space of dynamical variables, where each configuration of them has infinitely many different configurations that describe the exact same state of the exact same physical system, and you're integrating over all of them. What is this supposed to be? It's certainly not the integral over all possible physical paths, it's massively overcounting them and you have no way to control the manner in which it does that. 
The natural physical path integral is one integrating over each physically distinct configuration once. When we completely fix a gauge, this is exactly what the gauge fixing does: From all the possible equivalent configurations, the gauge condition picks one and only one representant, and we then wish to integrate over this space of representants, as it is the space of physically distinct configurations. Unfortunately, Gribov ambiguities mean that we can usually not do that throughout all of field configurations space and may be stuck defining the path integral only over a subset of physical configurations, a so-called Gribov region.
Therefore, it is unreasonable to expect there to be a path integral without fixing a gauge. The path integral, by its very purpose, must integrate over the space of all physically distinct configurations, and the way to achieve that in a gauge theory is some manner of gauge fixing, there is no way to evade this fact.
