A vector field chooses one tangent vector at every point. In the same way, we can choose a plane through every point: at each point p of \mathbb{R}^3, a two-dimensional subspace D_p of the tangent space T_p\mathbb{R}^3. Such a choice is called a plane field, or a distribution of rank 2 [1].
We call D smooth when, near every point, it is spanned by two smooth vector fields V_1 and V_2: at each point q their values are linearly independent and span D_q. The same definition gives distributions of any rank k on manifolds of any dimension, but planes in \mathbb{R}^3 already show everything that matters here.
In the figure, D is spanned by V_1 = \partial_x + \frac{z\varphi_x}{\varphi}\partial_z and V_2 = \partial_y + \frac{z\varphi_y}{\varphi}\partial_z, where \varphi(x, y) = 1 + 0.4\sin(1.5x)\cos(1.5y). Each small square lies in the plane at its center, and as the point p moves, its square tilts smoothly.
An integral surface of D is a surface N whose tangent plane at every point q of N is D_q. We say that D is integrable when every point lies on an integral surface.
For the field of the previous figure, the surfaces z = c\varphi(x, y) are integral surfaces. The tangent plane of z = c\varphi is spanned by \partial_x + c\varphi_x\partial_z and \partial_y + c\varphi_y\partial_z, and on this surface c = z/\varphi, so these two vectors are V_1 and V_2. The point (x, y, z) lies on the surface with c = z/\varphi(x, y), so these surfaces fill all of space.
In the figure the point p moves up and down, and the surface through p changes with it. For c = 0 it is the plane z = 0, and for c < 0 it is the mirror image of the surface for -c. At each point of the surface, the small square of D lies in the surface.
If we choose a line instead of a plane at every point, the question has a quick answer. Near any point, a line field is spanned by a vector field X that does not vanish, and the integral curves of X are tangent to the lines. These curves always exist, by the existence theorem for ordinary differential equations [1].
The vector field is not unique. For any positive function f, the field fX spans the same lines, and its integral curves follow the same paths at different speeds. In the figure two points move along one curve, and each velocity arrow lies on the line chosen at its point. The line field determines the curves, but not the speed along them.
For planes the same argument does not work. Following one vector field in D gives only a curve, and nothing guarantees that the curves of two fields fit together into a surface. The next figure shows a plane field with no integral surface at all.
Let D be spanned by X = \partial_x - \frac{y}{2}\partial_z and Y = \partial_y + \frac{x}{2}\partial_z. On the z-axis the plane is horizontal. At a point at distance r from the z-axis, the plane contains the horizontal direction pointing away from the axis, since xX + yY = x\partial_x + y\partial_y. It also contains -yX + xY = -y\partial_x + x\partial_y + \frac{r^2}{2}\partial_z, so it is tilted around that direction by the angle whose tangent is \frac{r}{2}.
Every horizontal half-line from the origin is tangent to D. On the half-line in the direction (\cos\alpha, \sin\alpha, 0), the vector field (\cos\alpha)X + (\sin\alpha)Y has no \partial_z part, so the half-line is one of its integral curves. An integral curve of a vector field in D that starts on an integral surface stays on it, so an integral surface through the origin would contain all these half-lines near the origin, that is, a disk in the plane z = 0.
But on that disk the plane of D is horizontal only at the center. Away from the center the squares tilt more and more, so the disk is not tangent to D, and no integral surface passes through the origin. The Lie bracket shows that the same happens at every point.
Suppose N is an integral surface and X and Y are vector fields with values in D. On N they are tangent to N, and the Lie bracket of two vector fields tangent to a surface is again tangent to it. So [X, Y] lies in D at every point of N.
For the twisting field, [X, Y] = (X(\frac{x}{2}) - Y(-\frac{y}{2}))\partial_z = \partial_z. The z-component of Y grows at the rate \frac{1}{2} in the x-direction, the z-component of X falls at the rate \frac{1}{2} in the y-direction, and the bracket adds the two rates.
But \partial_z is never in D. A vector aX + bY in D has horizontal part a\partial_x + b\partial_y, which is 0 only when a = b = 0. So no point lies on an integral surface of the twisting field.
A curve (x(t), y(t), z(t)) is tangent to D exactly when its velocity is \dot x X + \dot y Y, that is, when \dot z = \frac{1}{2}(x\dot y - y\dot x). The right side is the areal velocity of the point (x(t), y(t)) around the origin, the rate at which the segment from the origin to the point sweeps area in the xy-plane.
So a curve tangent to D can follow any path in the xy-plane, and its height changes by the area swept. Going once around a loop counterclockwise, it returns above its starting point, higher by the area the loop encloses. In the figure the loop is a circle of radius 1 through the origin, and z grows to \pi.
This gives another proof that there is no integral surface. Near a point, an integral surface would be the graph of a function over the xy-plane, since D never contains \partial_z. A curve tangent to D that starts on the graph stays on it, because its height is determined by its path in the xy-plane, so after going around a small loop it would return to its starting point instead of rising.
We call D involutive when the bracket of any two vector fields with values in D again has values in D. If V_1 and V_2 span D, it is enough to check that [V_1, V_2] lies in D, because for functions f and g we have [fV_1, gV_2] = fg[V_1, V_2] + f(V_1 g)V_2 - g(V_2 f)V_1, and every bracket of fields in D expands into terms of this kind. The previous scenes show that an integrable D is involutive.
The Frobenius theorem says that the converse is true: every involutive plane field is integrable [1]. More precisely, near every point there are coordinates (u, v, w) in which D is spanned by \partial_u and \partial_v, and each surface w = c is an integral surface. These surfaces do not meet each other, and together they fill a neighborhood of the point.
For example, away from the z-axis, let D be spanned by the horizontal direction pointing away from the axis and the vertical direction. The fields \frac{1}{r}(x\partial_x + y\partial_y) and \partial_z that span it have bracket 0, so D is involutive. Its integral surfaces are the half-planes bounded by the z-axis, and in the cylindrical coordinates (r, \theta, z) they become the planes \theta = c, as in the figure.
The theorem also decides when the equations u_x = F(x, y, u) and u_y = G(x, y, u) have a solution for every initial value. The graphs of solutions are the integral surfaces of the plane field spanned by \partial_x + F\partial_u and \partial_y + G\partial_u, and involutivity becomes the condition F_y + GF_u = G_x + FG_u.
Now we prove the theorem for an involutive plane field D in \mathbb{R}^3. Near a point p, after a linear change of coordinates, we may assume that D_p does not contain \partial_z. Then near p the projection (x, y, z) \mapsto (x, y) maps each plane D_q isomorphically onto the xy-plane.
So there is exactly one vector field V_1 with values in D whose projection is \partial_x, and exactly one V_2 whose projection is \partial_y. They have the form V_1 = \partial_x + a\partial_z and V_2 = \partial_y + b\partial_z, and together they span D. In the figure the planes above one point of the xy-plane tilt differently, but all of them project onto the same square spanned by \partial_x and \partial_y.
Their bracket is [V_1, V_2] = (V_1 b - V_2 a)\partial_z, so its projection is 0. Since D is involutive, [V_1, V_2] lies in D, and the only vector of D_q whose projection is 0 is the zero vector. Hence [V_1, V_2] = 0.
Move p to the origin, and let \Phi_s and \Psi_t be the flows of V_1 and V_2. Since [V_1, V_2] = 0, the flows commute: \Phi_s \circ \Psi_t = \Psi_t \circ \Phi_s [1].
Define F(s, t, c) = \Phi_s(\Psi_t(0, 0, c)). Then \partial F/\partial s = V_1, and because the flows commute, also \partial F/\partial t = V_2. At the origin \partial F/\partial c = \partial_z, and the vectors V_1, V_2 and \partial_z are linearly independent, so (s, t, c) are coordinates near the origin. In these coordinates V_1 = \partial_s and V_2 = \partial_t.
So D is spanned by \partial_s and \partial_t, and each surface on which c is constant is an integral surface. The figure constructs one of them: from the point q = (0, 0, c) it follows V_2, and then V_1 from each point reached. Following them in the other order (dashed) passes through the same points. With k commuting vector fields, the same argument proves the theorem for distributions of any rank on any manifold.
The opposite case is a field like the twisting one, in which X, Y and [X, Y] together span every direction. By the Chow-Rashevskii theorem, any two points of a connected manifold can then be joined by a curve tangent to D [2].
For the twisting field this is what the previous scenes suggested: a curve tangent to D can move in any horizontal direction, and going around loops changes its height. The figure shows the points reached from the origin by curves tangent to D of length at most L, measuring length with X and Y as orthonormal vectors. The region is a solid around the origin. Its width grows like L and its height like L^2, so for small L it is flat, and it thickens as L grows.
For an involutive field, a curve tangent to D never leaves the integral surface on which it starts, so the points it reaches form a piece of one surface. Which of the two cases occurs depends on the bracket.