Invariant submanifolds Let $M$ be a smooth manifold, and let $N$ be a submanifold.  Let $V$ be a smooth vector field on $M$ which generates a flow $\Phi_t$ on $M$.  My intuition tells me (perhaps modulo some technical assumptions) that the following is true:


  If $V(p)$ is tangent to $N$ for all $p\in N$, then $N$ is an invariant submanifold of $\Phi_t$.


Is this true? What sorts of technical assumptions would I need to worry about to make the statement rigorous?  I imagine, for example, that there could be global topological issues so that perhaps the statement only holds locally.

Is there a good (basic) reference on invariant submanifolds?
I am a physicist, not a mathematician, so my answer may lack some of the rigor you expect.  Still, since no one else has taken a stab at this problem in the last year, I will shed what little light I can.  

The reference you are looking for may be Sophus Lie's 1884 Differential Invariant Paper.  I was introduced to it in Oak Ridge in 1992 by Dr. Lawrence Dresner, an applied mathematician in Magnetics Division of the Y-12 lab.  Lie's work was being used to solve nonlinear PDE pertaining to superconductor stability problems in the Tokamak fusion reactor.

Sophus Lie was an 18th century Norwegian mathematician, but his style was rooted in the 17th century, so his work can be somewhat dense.  I recommend the translation by M Ackerman with the comments and additional material by Robert Hermann.  It was published in English in 1976 by MATH SCI PRESS, ISBN 0-915692-13-9.  The original paper is "On Differential Invariants", S. Lie, Math. Annallen, Vol. 24 (1884), 537-578.  I am just a student, and not really qualified to speak with authority on this issue, but as I understand it Lie's basic premise is precisely your statement.  

Sophus Lie was trying to do for differential equations what Evariste Galois did for polynomials.  A Lie Group is a group that preserves the structure of the smooth manifold.  Lie Group stabilizers, or "differential invariants" as he called them, form an embedded, invariant submanifold.  Because of this, DEQ's may be rewritten in terms of group stabilizers.  Because a function of stabilizers is itself a stabilizer, this puts the DEQ into the kernel of the map with the Lie algebra of the group which is used to find a solution. As Hermann says in the preface, "The key idea is that one should study the structure of the orbit space of a symmetry group on the space of solutions."

Again, pardon the lack of rigor, but perhaps you might benefit from a practical example. If a projectile is fired into a fluid the force of friction is proportional to the square of the velocity, so deriving an equation for how far a projectile penetrates a fluid as a function of time can be a challenge.
$$ F=-\alpha v^2 \rightarrow m\ddot{y}=-\alpha \dot{y}^2$$Here y is the penetration distance, m is mass and $\alpha$ is the drag coefficient, a unitless constant dependent on projectile shape and the density and viscosity of the fluid. ($\dot{y}=\frac{dy}{dt}$ and $\ddot{y}=\frac{d^2y}{dt^2}$) Also, y(0)=0 and $\dot{y}(0)=v_o$, the initial velocity of the projectile.  

This DEQ is invariant to the Lie group G(t,y)=$(\lambda t, \lambda^\beta y)\lambda_o=1$  Note that if $t'=\lambda t$ and $y'=\lambda^\beta y$, $\dot{y}'=\frac{dy'}{dt'}=\frac{\lambda^\beta dy}{\lambda dt}=\lambda^{\beta -1}\dot{y}$ and $\ddot{y}'=\lambda^{\beta -2}\ddot{y}$.  Applying these to the DEQ, 
$$m\lambda^{\beta -2}\ddot{y}=-\alpha \lambda^{2\beta -2}\dot{y}^2$$
For invariance, $\beta =0$.  Now find the infinitestimal transformations of the primed variables and use the method of characteristics to find the stabilizers (differential invariants).

$$\bigg(\frac{dt'}{d\lambda}\bigg)_{\lambda _o=1}=t$$
$$\bigg(\frac{dy'}{d\lambda}\bigg)_{\lambda _o=1}=\beta y$$
$$\bigg(\frac{d\dot{y}'}{d\lambda}\bigg)_{\lambda _o=1}=(\beta -1)\dot{y}$$
$$\bigg(\frac{d\ddot{y}'}{d\lambda}\bigg)_{\lambda _o=1}=(\beta -2)\ddot{y}$$
$$d\lambda=\frac{dt}{t}=\frac{dy}{\beta y}=\frac{d\dot{y}}{(\beta -1)\dot{y}}=\frac{d\ddot{y}}{(\beta -2)\ddot{y}}$$
$$\frac{dt}{t}=\frac{dy}{\beta y}\rightarrow \beta ln t=ln y + \mu \rightarrow \mu=\frac{y}{t^\beta}\bigg|_{\beta =0}=y$$
$$\frac{dt}{t}=\frac{d\dot{y}}{(\beta -1)\dot{y}} \rightarrow \nu=\frac{\dot{y}}{t^{\beta -1}}\bigg|_{\beta =0}=t\dot{y}$$
$$\frac{dt}{t}=\frac{d\ddot{y}}{(\beta -2)\ddot{y}} \rightarrow \eta=\frac{\ddot{y}}{t^{\beta -2}}\bigg|_{\beta =0}=t^2\ddot{y}$$
These constants of integration, $\mu$, $\nu$, and $\eta$, are group stabilizers for the group of differential equations, a subset of the group of polynomials.  There are an infinite number of them but we only need the first three for our DEQ.  These differential invariants form an embedded submanifold in the solution space, which means the DEQ may be rewritten in terms of the invariants.  Multiplying the DEQ by $t^2$, we have
$$mt^2\ddot{y}=-\alpha(t\dot{y})^2 \rightarrow m\eta=-\alpha \nu^2$$
Since $$t\frac{d\nu}{dt}=\eta-(\beta-1)\nu$$ we can use this expression of the Lie algebra to solve our DEQ.  $$t\frac{d\nu}{dt}=\nu -\frac{\alpha}{m}\nu^2 \rightarrow \frac{d\nu}{\nu-\frac{\alpha}{m}\nu^2}=\frac{dt}{t}$$
This separation of variables is a direct result of the group invariance.  With a little fractional decomposition, some integration, substitution of $\nu=t\dot{y}$, more integration and application of initial conditions, you get
$$y=\frac{m}{\alpha}ln\bigg(1+\frac{\alpha}{m}v_o t\bigg)
$$This solution is easily verified and provides accurate results.

Again, the applied math may not be exactly what you are looking for, but hopefully it will spark some thoughts.  On page 102 of the afore-mentioned book is Lie's Theorem 4.4.1.  "Every infinite continuous group determines an infinite sequence of differential invariants, which can be defined as the solutions of complete systems."  In the following comments by Hermann, he states "As far as I can tell, Lie (in the spirit of the 18th century) assumes that everything is suitably 'general', and that these facts are self evident.  As far as I know, they are not proved to this day!"  What was true for Hermann in 1976 may still hold.  The method obviously works so the proof should exist, but since nobody has answered your post in a year it may be possible that such a valuable proof is still out there waiting to be discovered.
I am a physicist, not a mathematician, so my answer may lack some of the rigor you expect.  Still, since no one else has taken a stab at this problem in the last year, I will shed what little light I can.  

The reference you are looking for may be Sophus Lie's 1884 Differential Invariant Paper.  I was introduced to it in Oak Ridge in 1992 by Dr. Lawrence Dresner, an applied mathematician in Magnetics Division of the Y-12 lab.  Lie's work was being used to solve nonlinear PDE pertaining to superconductor stability problems in the Tokamak fusion reactor.

Sophus Lie was an 18th century Norwegian mathematician, but his style was rooted in the 17th century, so his work can be somewhat dense.  I recommend the translation by M Ackerman with the comments and additional material by Robert Hermann.  It was published in English in 1976 by MATH SCI PRESS, ISBN 0-915692-13-9.  The original paper is "On Differential Invariants", S. Lie, Math. Annallen, Vol. 24 (1884), 537-578.  I am just a student, and not really qualified to speak with authority on this issue, but as I understand it Lie's basic premise is precisely your statement.  

Sophus Lie was trying to do for differential equations what Evariste Galois did for polynomials.  A Lie Group is a group that preserves the structure of the smooth manifold.  Lie Group stabilizers, or "differential invariants" as he called them, form an embedded, invariant submanifold.  Because of this, DEQ's may be rewritten in terms of group stabilizers.  Because a function of stabilizers is itself a stabilizer, this puts the DEQ into the kernel of the map with the Lie algebra of the group which is used to find a solution. As Hermann says in the preface, "The key idea is that one should study the structure of the orbit space of a symmetry group on the space of solutions."

Again, pardon the lack of rigor, but perhaps you might benefit from a practical example. If a projectile is fired into a fluid the force of friction is proportional to the square of the velocity, so deriving an equation for how far a projectile penetrates a fluid as a function of time can be a challenge.
$$ F=-\alpha v^2 \rightarrow m\ddot{y}=-\alpha \dot{y}^2$$Here y is the penetration distance, m is mass and $\alpha$ is the drag coefficient, a unitless constant dependent on projectile shape and the density and viscosity of the fluid. ($\dot{y}=\frac{dy}{dt}$ and $\ddot{y}=\frac{d^2y}{dt^2}$) Also, y(0)=0 and $\dot{y}(0)=v_o$, the initial velocity of the projectile.  

This DEQ is invariant to the Lie group G(t,y)=$(\lambda t, \lambda^\beta y)\lambda_o=1$  Note that if $t'=\lambda t$ and $y'=\lambda^\beta y$, $\dot{y}'=\frac{dy'}{dt'}=\frac{\lambda^\beta dy}{\lambda dt}=\lambda^{\beta -1}\dot{y}$ and $\ddot{y}'=\lambda^{\beta -2}\ddot{y}$.  Applying these to the DEQ, 
$$m\lambda^{\beta -2}\ddot{y}=-\alpha \lambda^{2\beta -2}\dot{y}^2$$
For invariance, $\beta =0$.  Now find the infinitestimal transformations of the primed variables and use the method of characteristics to find the stabilizers (differential invariants).

$$\bigg(\frac{dt'}{d\lambda}\bigg)_{\lambda _o=1}=t$$
$$\bigg(\frac{dy'}{d\lambda}\bigg)_{\lambda _o=1}=\beta y$$
$$\bigg(\frac{d\dot{y}'}{d\lambda}\bigg)_{\lambda _o=1}=(\beta -1)\dot{y}$$
$$\bigg(\frac{d\ddot{y}'}{d\lambda}\bigg)_{\lambda _o=1}=(\beta -2)\ddot{y}$$
$$d\lambda=\frac{dt}{t}=\frac{dy}{\beta y}=\frac{d\dot{y}}{(\beta -1)\dot{y}}=\frac{d\ddot{y}}{(\beta -2)\ddot{y}}$$
$$\frac{dt}{t}=\frac{dy}{\beta y}\rightarrow \beta ln t=ln y + \mu \rightarrow \mu=\frac{y}{t^\beta}\bigg|_{\beta =0}=y$$
$$\frac{dt}{t}=\frac{d\dot{y}}{(\beta -1)\dot{y}} \rightarrow \nu=\frac{\dot{y}}{t^{\beta -1}}\bigg|_{\beta =0}=t\dot{y}$$
$$\frac{dt}{t}=\frac{d\ddot{y}}{(\beta -2)\ddot{y}} \rightarrow \eta=\frac{\ddot{y}}{t^{\beta -2}}\bigg|_{\beta =0}=t^2\ddot{y}$$
These constants of integration, $\mu$, $\nu$, and $\eta$, are group stabilizers for the group of differential equations, a subset of the group of polynomials.  There are an infinite number of them but we only need the first three for our DEQ.  These differential invariants form an embedded submanifold in the solution space, which means the DEQ may be rewritten in terms of the invariants.  Multiplying the DEQ by $t^2$, we have
$$mt^2\ddot{y}=-\alpha(t\dot{y})^2 \rightarrow m\eta=-\alpha \nu^2$$
Since $$t\frac{d\nu}{dt}=\eta-(\beta-1)\nu$$ we can use this expression of the Lie algebra to solve our DEQ.  $$t\frac{d\nu}{dt}=\nu -\frac{\alpha}{m}\nu^2 \rightarrow \frac{d\nu}{\nu-\frac{\alpha}{m}\nu^2}=\frac{dt}{t}$$
This separation of variables is a direct result of the group invariance.  With a little fractional decomposition, some integration, substitution of $\nu=t\dot{y}$, more integration and application of initial conditions, you get
$$y=\frac{m}{\alpha}ln\bigg(1+\frac{\alpha}{m}v_o t\bigg)
$$This solution is easily verified and provides accurate results.

Again, the applied math may not be exactly what you are looking for, but hopefully it will spark some thoughts.  On page 102 of the afore-mentioned book is Lie's Theorem 4.4.1.  "Every infinite continuous group determines an infinite sequence of differential invariants, which can be defined as the solutions of complete systems."  In the following comments by Hermann, he states "As far as I can tell, Lie (in the spirit of the 18th century) assumes that everything is suitably 'general', and that these facts are self evident.  As far as I know, they are not proved to this day!"  What was true for Hermann in 1976 may still hold.  The method obviously works so the proof should exist, but since nobody has answered your post in a year it may be possible that such a valuable proof is still out there waiting to be discovered.