EDIT: I corrected some typos in my solution and wrote some details more precise (I was a bit sloppy first as I originally did not plan to discuss the full solution, but rather some ideas whether a solution can be constructed somewhat more naturally. E.g. using the Lie algebra and Lie group operations for $\mathbb{S}^3$.)

I am trying to solve the equation $$ \mathrm{d}G(x,y) = \mathrm{Vol}(x)+(-1)^{n}\mathrm{Vol}(y):= H $$ for $G\in \Omega^{n-1}\bigl((\mathbb{S}^n\times \mathbb{S}^n)\backslash \Delta=:M\bigr)$. Here $\mathrm{Vol}$ is the standard volume form on $\mathbb{S}^n$ and $\Delta$ the diagonal. The solution for $n=1$ is $$ G(x,y)= - \alpha(x,y)\quad \text{for all }x,y\in \mathbb{S}^1 $$ where $\alpha(x,y)$ is the counterclockwise angle from $x$ to $y$. I can find a solution for general $n\in \mathbb{N}$ as $$ G = \int_{[0,1]} \psi^*H $$ for a contraction $\psi: [0,1]\times N \rightarrow N$ of an open thickening $N$ of $M$ onto the antidiagonal $\overline{\Delta}:=\{(x,-x)\}$ (on $\overline{\Delta}$ holds namely $H=0$). Here I use a closed extension of $\mathrm{Vol}$ to $\mathbb{R}^{n+1}_{\neq 0}$. However, the solutions I obtain (in coordinates $x^i$ on $\mathbb{R}^{n+1}$) look too complicated and it is computationally hard to work with it further.

Does anybody have any idea how to solve this equation "nicely" e.g. by some natural construction of $G(x,y)$ similar to the case $n=1$?

I am mainly interested in the case $n=3$. Can the Lie group structure be used to construct a solution?

Thanks for any ideas!

Big picture: $H$ is a smooth integral kernel of the orthogonal projection to harmonic forms and $G$ is a (singular) integral kernel of the Green operator (:=the inverse of $\mathrm{d}$)

UPDATE1: As Robert Bryant suggests, the sign is indeed $(-1)^n$.

UPDATE2: Reading the nice answer of Robert Bryant, I just typed my complicated solution for a comparison. Perhaps it is the same except Robert using a clever notation: For every $n$ a solution is

$$ G(x,y)=\sum_{k=0}^{n-1} g_k(x\cdot y) \omega_k(x,y) $$ where $$ g_k(u) = \int_{0}^1 \frac{t^k(t-1)^{n-1-k}}{(2t(t-1)(1+u)+1)^{\frac{n+1}{2}}}\mathrm{dt} $$ can be solved inductively and $$ \omega_k(x,y) = \frac{1}{k!(n-1-k)!} \sum_{\sigma\in \mathbb{S}_{n+1}} (-1)^\sigma x_{\sigma_1}y_{\sigma_2}\mathrm{d}x_{\sigma_3}\ldots\mathrm{d}x_{\sigma_{2+k}}\mathrm{d}y_{\sigma_{3+k}}\ldots\mathrm{d}y_{\sigma_{n+1}}$$

P.S. I am still tempted to say something like $G(x,y)((v_1,w_1),\ldots,(v_{n-1},w_{n-1}))$ is the angle between the hyperplane $\langle w_1,\ldots,w_{n-1}\rangle$ at $y$ and the parallel transport of the hyperplane $\langle v_1,\ldots,v_{n-1}\rangle$ from $x$ to $y$ along the shortest arc on the unique great circle through $x$ and $y$. This might be a complete nonsense as well... If my computations are correct then for $n=2$ holds $$ G(e_1,\cos(\alpha)e_1 + \sin(\alpha)e_2) = \frac{\sin(\alpha)}{\cos(\alpha)-1}(\mathrm{d}y_3 - \mathrm{d} x_3 ). $$

linearmap $$G_{(u,v)}:T_u\mathbb{S}^2\oplus T_v\mathbb{S}^2\to\mathbb{R},$$ which your proposed formula is not. [This comment was made before your most recent edit, where you removed your originally proposed formula for $G$ when $n=2$.] $\endgroup$ – Robert Bryant Jan 27 '18 at 15:33