On-shell symmetry from a path integral point of view Normally supersymmetric quantum field theories have Lagrangians which are supersymmetric only on-shell, i.e. with the field equations imposed. In many cases this can be solved by introducing auxilary fields (field which don't carry dynamical degrees of freedom, i.e. which on-shell become a function of the other fields). However, there are cases where no such formulation is known, e.g. N=4 super-Yang-Mills in 4D.
Since the path integral is an integral over all field configurations, most of them off-shell, naively there is no reason for it to preserve the on-shell symmetry. Nevertheless the symmetry is preserved in the quantum theory.
Of course it is possible to avoid the problem by resorting to a "Hamiltonian" approach. That is, the space of on-shell field configurations is the phase space of the theory and it is (at least formally) possible to quantize it. However, one would like to have an understanding of the symmetries survival in a path integral approach. So:

How can we understand the presence of on-shell symmetry after quantization from a path integral point of view?

 A: 
How can we understand the presence of on-shell symmetry after quantization from a path integral point of view?

One can derive a Schwinger-Dyson equation associated with the current conservation, also known as a Ward identity; see e.g. Peskin and Schroeder, An Introduction to Quantum Field Theory, Section 9.6; or Srednicki, Quantum Field Theory, Chapter 22.
A: Quantum field theoretic analog of the classical equations of motion (EOM):
$$ -i\, \delta S (x)\, \phi_{1}(x_1) ... \phi_{n}(x_n) = \sum _{j=1}^n  \phi_{1}(x_1)... \delta \phi_j (x)\delta(x-x_j)... \phi_{n}(x_n)=  \delta\phi_{1}(x_1)\delta(x-x_1)\phi_2(x_2)... \phi_{n}(x_n)+\phi_1(x_1)\delta\phi_2(x_2)\delta(x-x_2)\phi_3(x_3)...\phi_n(x_n)+...+\phi_1(x_1)...\phi_{n-1}(x_{n-1})\delta\phi_n(x_n)\delta(x-x_n)$$
where
$$ \delta S\equiv \sum_{i=1}^n{\delta S\over \delta \phi_{i} (x)}\delta \phi_i(x) $$
that is, $${\delta S\over \delta \phi_{i} (x)}$$ is the Euler-Lagrange EOM for the field $\phi_{i} (x)$.
Answer to comment: In case you have local operators $O_i (x_i)$, this is the equation:
$$ -i\, {\delta S (x)\over \delta\phi (x)} O_{1}(x_1) ... O_{n}(x_n) =  {\delta O_{1}(x_1)\over \delta \phi(x)}O_2(x_2)... O_{n}(x_n)+...+O_1(x_1)...O_{n-1}(x_{n-1}){\delta O_n(x_n) \over \delta \phi (x)}$$
A: The point of view you give of quantizing phase space is not so good. If you have a classical phase space, there will be more than one quantum analog, differing by commutator terms that vanish classically. I never found the idea that you start with a classical system particularly lucid--- the quantum system is defined by the path integral, and the answer must always be found in the path integral itself.
In 4d, as you said in the comments, the problem is sidestepped by using N=1 superspace plus the appropriate R symmetry, which rotate the different SUSY's into each other, ensuring that the full SUSY is there. In order for the theory to be SUSY, you just have to check that the theory is N=1 and that the R-symmetry is good. This is usually manifest, since the R symmmetry is simple global rotation of fields into each other. But this doesn't work in 10 or 11 dimensions.
A part of the answer suggested by Qmechanic, the Schwinger Dyson equation (Heisenberg equation of motion) is obeyed by the field operators, so the variation really does vanish as an operator, in a certain sense described precisely below. But this is not a complete answer because it isn't true that everything that vanishes classically when the equations of motion are satisfied also vanishes quantum mechanically, in the path integral. The reason is that the quantities in the path integral can be multiplied at the same point in space and time, and there are singularities of the coinciding operator products which make the naive equations of motion fail in exactly the type of expressions that occur in the action.
For a trivial example, consider the product of two X operators
$$\langle x(t) x(t')\rangle = 0 $$
In the trivial Euclidean free-particle (Brownian motion) 1-d (quantum mechanics, not field theory) path integral $S= \int \dot{x}^2$. The operator X obeys the equation of motion $\ddot{X}=0$, so the second derivative is zero, but of course it isn't at the coinciding point (or else propagator would vanish at all times with the reasonable boundary conditions). The correct second derivative of the correlation function is
$$ \langle \ddot{x}(t) x(t')\rangle = \delta(t-t')$$
The equation of motion holds except at coincident points. The same is true in field theories for the same reason--- the propagator is sourced at the point where the two operators coincide by a delta-function with unit amplitude (if the fields are canonically normalized).
Further, since $\ddot{X}=0$, you might think that the following transformation is a symmetry
$$ \dot x + \epsilon f(t)\ddot{X} $$
Since the variation of S to first order in $\epsilon$ is $\int \dot{x}\ddot{x}f(t)$. So naively, using the equation of motion $\ddot{x}=0$, you would naively find that this is a symmetry. But this is obvious nonsense--- if f is nonconstant (if f is constant, the transformation is a time translation) this is obviously not a symmetry of the action. since it is a time translation by a different amount at different times, and there is a unit of time defined by the diffusion.
The reason is that the $\ddot{x}$ is multiplied by a coinciding $\dot{x}$, and when there are coinciding quantities, the equations of motion are not necessarily satisfied. The actual variation in S is
$$ \int f(t) {d\over dt} {\dot{x}^2\over 2}$$
and you see that it is a perfect derivative if f is constant, assuming Stratonovich convention for time derivatives--- centered differences--- and this is the symmetry of time translation in this case), but it is not a symmetry if f is not constant.
When are equations of motion satisfied in a path integral?
When you have a path integral:
$$ \int e^{iS} D\phi $$
The integral is invariant under a shift of integration variables $\phi(x)+\delta\phi(x)$, where $\delta\phi$ is an arbitrary function of x, since each integral is separately translation invariant. There is no determinant for this transformation--- it's just a shift in the integration variable by a constant at each time. 
The change in the integrand when you do the shift is, given by expanding the thing to first order.
$$ \int e^{iS(\phi+\delta\phi)} D\phi = \int e^{iS} D\phi + i\int (\int {\delta S\over \delta \phi(x)}\delta\phi(x) d^dx) e^{iS} D\phi $$
and from this, you learn that the equations of motion are satisfied. This is the correct demonstration of the equations of motion.
$$ {\delta S\over \delta\phi} = 0 $$
This is the equation of motion for the $\phi$ field only, only one field. Now suppose you have some other insertions (composite operators $O_k$):
$$ \int O_1(x_1) O_2(x_2) ... O_n(x_n) e^{iS} D\phi $$
Now doing the same shift, to first order, you find
$$ \int (O_1 +\delta O_1(x_1)) ... (O_n+\delta O_n(x_n)) e^{iS} (1+i\delta S) D\phi = \int e^{iS} D\phi $$
So you find that the equation of motion for the one field $\phi$ is satisfied even with composite operator insertions except at coinciding points:
$$ -i\, {\delta S\over \delta \phi(x)} \, O_{1}(x_1) ... O_{n}(x_n) = {\delta O_1 (x_1)\over \delta \phi(x)} O_2(x_2) ... O_n(x_n) + O_1(x_1){\delta O_2(x_2)\over \delta\phi(x)}...O_n(x_n) + ... + O_1(x_1)O_2(x_2) ... O_{n-1}(x_{n-1}){\delta O_n(x_n)\over \delta \phi(x)}.$$
This is the operator equation for the equation of motion multiplying arbitrary expressions of other fields. The expression is known as the Schwinger-Dyson equation.
It is not just true as a vacuum expectation value, because it is satisfied independent of the path integral boundary conditions, which don't have to be vacuum, and you don't need to do a path integral from time minus infinity to time infinity, the identity holds in an arbitrary state.
This Schwinger-Dyson argument for insertions of any type with lots of fields running around implies the following rules for using the equation of motion:


*

*The equations of motion derived from varying $\phi$ are satisfied as operator equations away from insertions. 

*the equation of motion derived from varying $\phi$ are satisfied even at cross-insertions of other fundamental fields (independent path-integrations) other than $\phi$.

*the equations of motion fail only at insertions of composite or elementary operators which vary when you shift the field whose equation of motion you are talking about, and the failure is of the ward-identity type. 


So it isn't true that anything that vanishes classically has no effect. But when you have a product of the equation of motion multiplying functions of fields which are different from the one you vary to get the equation of motion, these are still identically zero.
In the SUSY case, your variation involves products of fields with their SUSY partners, which are independent integration variables. The equations of motion are satisfied in the products that you reduce, so using the equations of motion in the SUSY closure is justified, but you check it theory by theory by hand, by looking to see which equation of motion you use, and what it is multiplying.
The interpretation is that when you do the SUSY transformation, you also have to slide the fields around a little bit in a measure preserving way to finish the transformation.
A: First of all there is an error in the Schwinger-Dyson equation that was derived in the previous post.  There will be n-1 fields on the right-hand side.  In any event, the SD equations don't save you because the supercurrent has fields sitting at the same point, so these contact terms come in and spoil the equations of motion.  I.e., you get the equation of motion when you differentiate the supercurrent, but it is multiplied by fields at the same points, so all these delta functions come in.  Try it at home, you'll see what I mean.
