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 ). $$



A primitive to $\mathrm{Vol}(x)\pm \mathrm{Vol}(y)$ on $(\mathbb{S}^n \times \mathbb{S}^n)\backslash \Delta$
**NB:** *I have edited my answer to remove the comment at the beginning about correcting the original sign of the OP in the formula for $H$, since the OP has now corrected the erroneous sign in the question. I have also added one remark at the end about the case when $n$ is even.*


Now, one has
$$
H = x^\*\Omega + (-1)^n\,y^\*\Omega
\tag1
$$
where $x:\mathbb{S}^n\times \mathbb{S}^n \to \mathbb{S}^n\subset \mathbb{R}^{n+1}$ and $y:\mathbb{S}^n\times \mathbb{S}^n \to \mathbb{S}^n \subset \mathbb{R}^{n+1}$ are the projections on the first and second factors and $\Omega$ is the $\mathrm{SO}(n{+}1)$-invariant volume form on $\mathbb{S}^n$ corresponding to a fixed orientation and inner product on $\mathbb{R}^{n+1}$.


Thus, to remove any sign ambiguities, fix an orientation and inner product on $\mathbb{R}^{n+1}$, and let $(u\_0,\ldots,u\_n)$ be an oriented orthonormal basis of $\mathbb{R}^{n+1}$, so that $U = u\_0\wedge\cdots\wedge u\_n$ is a basis of $\Lambda^{n+1}(\mathbb{R}^{n+1})$. Thus, for example, one has 
$$
\tfrac1{n!}\,x\wedge (\mathrm{d}x)^n = x^\*\Omega\,U
\quad\text{and}\quad \tfrac1{n!}\,y\wedge (\mathrm{d}y)^n = y^\*\Omega\,U
\tag2
$$
To save writing, for a form $\phi$ with values in $\Lambda^{n+1}(\mathbb{R}^{n+1})$, I will henceforth write $[\phi]$ for the scalar-valued differential form such that $\phi = [\phi]\, U$. Thus, $x^\*\Omega = [\tfrac1{n!}\,x\wedge (\mathrm{d}x)^n]$, etc. 


Using this notation, for example, one has, when $n=1$, that $\mathrm{d}G = H$
where 
$$
G = [\,f(x{\cdot}y)\,x\wedge y\,],\tag3
$$
where $f:[-1,1)\to\mathbb{R}$ is the (unique) smooth function that satisfies
$$
f(\cos t) = \frac{(\pi-t)}{\sin t}\quad\text{for}\quad 0<t<\pi,
$$
as one easily verifies by parametrizing as $x =(\cos\theta,\sin\theta)$ and $y=(\cos\phi,\sin\phi)$. Note that, since $x{\cdot}y<1$ on $M$, $G$ is a smooth function on $M$ and that $G$ vanishes when pulled back to the anti-diagonal, where $x+y=0$.


When $n=2$, calculation yields that 
$$
G = \frac{[x\wedge y\wedge \mathrm{d}(x{-}y)]}{1-x{\cdot}y}\tag4
$$
satisfies $\mathrm{d}G=H$, and, when $n=3$ (the case that the OP is particularly interested in), one finds that
$$
G = \left[\,x\wedge y \wedge \bigl(\,f(x{\cdot}y)\,(\mathrm{d}x)^2
-\left(\tfrac12 + (x{\cdot}y)\,f(x{\cdot}y)\right)\mathrm{d}x\wedge\mathrm{d}y
+f(x{\cdot}y)\,(\mathrm{d}y)^2\bigr)\,\right]\tag5
$$
satisfies $\mathrm{d}G = H$, where $f:[-1,1)\to\mathbb{R}$ is the (unique) smooth function that satisfies
$$
f(\cos t) = \frac{\pi-t+\sin t\cos t}{4\sin^3t} \quad\text{for}\quad 0<t<\pi.\tag6
$$
(It may not be apparent that $f$ is smooth at $-1$, but it is, once one sets $f(-1) = \tfrac16$.)


The general pattern is clear, i.e., that, for general $n$, there will be an expression for $G$ of the form
$$
G = \left[x\wedge y\wedge \left(\,f\_0(x{\cdot}y)\,(\mathrm{d}x)^{n-1} 
 + f\_1(x{\cdot}y)\,(\mathrm{d}x)^{n-2}{\wedge}\mathrm{d}y + \cdots 
 + f\_{n-1}(x{\cdot}y)\,(\mathrm{d}y)^{n-1}\,\right)\right]\tag7
$$
where the $f\_i$ are smooth functions on $[-1,1)$ that are determined by the condition that $\mathrm{d}G = H$. These coefficients can be determined recursively, but the general answer does not seem to be particularly illuminating, so I won't include it here.


However, I will remark that when $n$ is even, the functions $f\_i$ turn out to be rational functions of $x{\cdot}y$. For example, when $n=4$, one finds
$$
G = \frac{\left[\,x\wedge y\wedge\bigl(\,(\mathrm{d}x-\mathrm{d}y)^3
+(1-x{\cdot}y)\,((\mathrm{d}x)^3-(\mathrm{d}y)^3)\bigr)\,\right]}
{18\,(1-x{\cdot}y)^2}
$$


