Upon reading this answer to this question, the last paragraph mentions the following. "Requiring the [almost complex] structure to be integrable corresponds to a certain PDE for this map." The map in question, I presume, can be taken to be some $f: M \to Gr(3, \mathbb{C}^\infty)$ for the almost complex $6$-dimensional manifold $M$, which corresponds to $TM$ by pullback of the tautological bundle.
I'm no expert in these topics so I am struggling to unpack this claim. What is this PDE and how do we obtain it? Furthermore, is there a reference to Yau's mentioned program or a survey of attempts/results in that direction?
EDIT
I've emailed Professor Yau and to my delight he responded. Here is what I've learned. The map $f$ above determines an almost complex structure $J$ on $TM$. Then, as Moishe mentioned in the comments, we have that the Newlander-Nirenberg theorem gives the integrable condition through the vanishing of the Nijenhuis tensor $N_J$, which takes as inputs two vector fields $X,Y$. In particular, we may write $$ N_J : \Gamma(TM) \times \Gamma(TM) \to \Gamma(TM) $$ which is the partial differential operator in question. Yau mentions that $N_J$ can be written in terms of the differential of $f$. This part eludes me.
So my updated question is, how can $N_J$ be written in terms of $df$?
EDIT 2
Perhaps the PDE in question relating to the Nijenhuis tensor is different than I explained, and it's "simply" a question of working locally. I'll try to sketch out the idea here. First, we have
$$ N_J (X,Y) = [X,Y] + J([JX,Y]+[X,JY]) - [JX,JY], $$
and writing it in local coordinates, this gives
$$ X = \sum_{i=1}^{n} X^i \partial_i, ~~~~ Y = \sum_{i=1}^{n} Y^i \partial_i ,$$
$$ N_J (X,Y) = \sum_{i=1}^{n} \sum_{j=1}^{n} \Big[ (X^j \partial_j Y^i - Y^j \partial_j X^i)\partial_i $$ $$ +J \Big( ((JX)^j \partial_j Y^i - Y^j \partial_j (JX)^i) \partial_i \Big) $$ $$ +J \Big( (X^j \partial_j (JY)^i - (JY)^j \partial_j X^i) \partial_i \Big) $$ $$ -\Big( (JX)^j \partial_j (JY)^i - (JY)^j \partial_j (JX)^i \Big)\partial_i \Big]. $$
The first and last term of the big equation are in terms of $\partial_i$ while the two middle terms are in terms of $J(\partial_i)$. Assuming that $\partial_i$ and $J(\partial_i)$ are linearly independent for all $i$, setting $N_J(X,Y) = 0$ gives us a system of two first order PDEs with $2n$ variables $X^i, Y^i,$ namely
$$ \sum_j (X^j \partial_j Y^i - Y^j \partial_j X^i) - ( (JX)^j \partial_j (JY)^i - (JY)^j \partial_j (JX)^i ) = 0$$ and $$ \sum_j ((JX)^j \partial_j Y^i - Y^j \partial_j (JX)^i) + (X^j \partial_j (JY)^i - (JY)^j \partial_j X^i) = 0. $$
We can (hopefully) solve this locally and then use a partition of unity argument to look for a global set of solutions.
The question still stands for how $N_J$ (or in particular the linear $J$) can be written in terms of $df$, but does the above make sense?
