Let $D$ be a $2$ dimensional distribution of $\mathbb{R}^3$.
Is there a $1$ dimensional foliation of $\mathbb{R}^3$ with Frenet curves such that for every leaf $\gamma$ of the foliation we have $\mathrm{span}(\gamma'(t), \gamma ''(t))=D(\gamma(t))$ where $\gamma(t)$ is the unit speed parametrization of the leaf $\gamma$?


One can generalize the question when $\mathbb{R}^3$ is equiped with an arbitrary Riemannian metric and we require a foliation with $\mathrm{span}(\gamma'(t), \nabla\_{\gamma'(t)} \gamma '(t))=D(\gamma(t))$.



Can we foliate the space $\mathbb{R}^3$ with Frenet curves whose tangent and normal vectores span a given $2$ dimensional distribution?
The answer is 'no' in general for an arbitrary Riemannian metric $g$ on a $3$-manifold $M$ and $2$-plane field $D\subset TM$.


I'll give the argument for the flat metric on $\mathbb{R}^3$ and leave the (easy) generalization to arbitrary metrics for the interested.


Let $M$ be $\mathbb{R}^3$ endowed with the flat metric, also fix an orientation on $M$. Every $2$-plane field on $M$ is orientable, so when one fixes a $2$-plane field $D$ on $M$, one might as well fix an orientation of the $2$-plane field, which means that $D$ (and its orientation) can be specified by specifying its oriented unit normal vector field, say $\mathbf{b}:M\to S^2$. What the OP is asking is whether one can find an orthogonal unit vector field $\mathbf{t}:M\to S^2$ such that the integral curves of $\mathbf{t}$ are curves that have $\mathbf{b}$ as their Frenet bi-normal, i.e., that the principal normals of these integral curves are multiples of $\mathbf{n} = \mathbf{b}\times\mathbf{t}$.


Suppose that $\mathbf{b}:M\to S^2$ is given and let $B\subset M\times S^2\times S^2\times S^2$ consist of the oriented orthonormal frames $(x,e\_1,e\_2,e\_3)$ that satisfy $e\_3 = \mathbf{b}(x)$. Then the projection $x:B\to M = \mathbb{R}^3$ makes $B$ into a circle bundle over $M$, and we have the usual *first structure equations*
$$
\mathrm{d}x = e\_i\,\omega\_i = e\_1\,\omega\_1 + e\_2\,\omega\_2 + e\_3\,\omega\_3
\qquad\text{and}\qquad 
\mathrm{d}e\_i = e\_j\,\omega\_{ji}\,,
$$
where $\omega\_{ij}=-\omega\_{ji}$ and the $\omega\_i$ satisfy the *second structure equations*
$$
\mathrm{d}\omega\_i = -\omega\_{ij}\wedge\omega\_j
\qquad\text{and}\qquad
\mathrm{d}\omega\_{ij} = - \omega\_{ik}\wedge\omega\_{kj}\,.
$$
By construction, the $1$-forms $\omega\_1,\omega\_2,\omega\_3,\omega\_{12}$ are a basis for the $1$-forms on $B$, while we must have functions $a\_{ir}$ (where $i=1,2$ and $r=1,2,3$) such that 
$$
\omega\_{i3} = a\_{ir}\,\omega\_r\,
$$
Now, the choice of a frame field $(\mathbf{t},\mathbf{n},\mathbf{b}):M\to S^2\times S^2\times S^2$ with the above properties, i.e., 
$$
\mathrm{d}\mathbf{t}\equiv 0
\quad \text{modulo}\quad \mathbf{n},\omega\_2,\omega\_3\,,
$$  

is really a choice of a section of $B$ that falls into the locus where 
$$
\omega\_{31} = -\omega\_{13} = -a\_{1r}\,\omega\_r \equiv -a\_{11}\,\omega\_1 \equiv 0
\quad\text{modulo}\quad \omega\_2,\omega\_3\,,
$$
i.e, in the locus where $a\_{11}=0$. Thus, we need to know how $a\_{11}$ varies on the fibers of $B\to M$ in order to understand its zero locus.


Now, a straightforward calculation using the structure equations shows that
$$
\mathrm{d}(a\_{11}+a\_{22}) \equiv \mathrm{d}(a\_{12}-a\_{21}) \equiv 0
\quad\text{modulo}\quad \omega\_1,\omega\_2,\omega\_3\,,
$$
while
$$\left.
\begin{aligned}
\mathrm{d}(a\_{11}{-}a\_{22}) \equiv -2(a\_{12}{+}a\_{21})\,\omega\_{12}\\
\mathrm{d}(a\_{12}{+}a\_{21}) \equiv \phantom{-}2(a\_{11}{-}a\_{22})\,\omega\_{12}
\end{aligned}\right\} \quad\text{modulo}\quad \omega\_1,\omega\_2,\omega\_3\,.
$$


It follows that there are unique functions $r$ and $s$ on $B$ that are constant on the fibers of $B\to M$ (and hence are well-defined on $M$) that satisfy $r = \tfrac12(a\_{11}{+}a\_{22})$ and $s\ge0$
$$
4s^2 = (a\_{11}{-}a\_{22})^2 + (a\_{12}{+}a\_{21})^2.
$$


It follows, in particular, that, over the open set $U$ in $M$ where $|r|>s$, the function 
$$
a\_{11} = r + \tfrac12(a\_{11}{-}a\_{22})
$$
never vanishes. Hence, there cannot be a section of $B\to M$ over $U$ with the desired properties. It is not difficult to write down examples of $2$-plane fields on $M$ where this inequality holds, and hence such $2$-plane fields will not have the desired properties.


Meanwhile, over the open set $V$ in $M$ on which $|r| < s$, one finds that the equation $a\_{11}=0$ defines four distinct sections of $B\to M$ over $V$ that have the desired property.


The boundary cases, where $|r|=s$ require further analysis.


Thus, whether there exists a foliation of the desired kind depends at least on a certain inequality in the first derivatives of the $2$-plane field being satisfied.


**Added Remark:** A 'geometric' interpretation of this obstruction can be made as follows: On $B$ (which is determined by the unit vector field $\mathbf{b}:M\to S^2$) quadratic form
$$
Q = a\_{11}\,{\omega\_1}^2 + (a\_{12}{+}a\_{21})\,\omega\_1{\circ}\omega\_2 + a\_{22}\,{\omega\_2}^2
$$ 
turns out to be the pullback under the submersion $x:B\to M$ of a quadratic form $Q\_{\mathbf{b}}$ on the $2$-plane field $D = \mathbf{b}^\perp\subset TM$, and what one finds is that, if $\mathbf{t}:M\to S^2$ is actually a section of $D$ such that its integral curves have their principal normals also lying in $D$, then $\mathbf{t}$ must be a null vector for the quadratic form $Q\_{\mathbf{b}}$. Thus, if $Q\_{\mathbf{b}}$ is definite (positive or negative), then there is no solution, while if $Q\_{\mathbf{b}}$ is of split type, there are two distinct null directions at each point, and these define the two possible line fields with the desired property of their principal normals.


The interesting case is what happens when $Q\_{\mathbf{b}}$ vanishes identically, which, of course, can happen, but is a very special case (in particular, the plane field $D$ has to be Frobenius, but even this is not sufficient). Then there is much more that can happen with the foliations.


