Solution (source code)

= Solution

Let $K_{ij}=\partial_j u_i$, $C_{ij}=\langle R_iR_j\rangle$, $\alpha=2H/\zeta$ and $D_R=2k_BT/\zeta$. The configuration equation of the <Hookean dumbbell model> becomes
$$
\partial_t\psi=-\nabla_{\mathbf R}\cdot[(K\mathbf R-\alpha\mathbf R)\psi]+D_R\nabla_{\mathbf R}^2\psi.
$$
Here the paper's $k$ denotes the <Boltzmann constant> $k_B$. Integrating the <second moment> equation against $R_iR_j$ and discarding boundary terms at infinite separation gives
$$
\dot C_{ij}=K_{i\ell}C_{\ell j}+K_{j\ell}C_{i\ell}-2\alpha C_{ij}+2D_R\delta_{ij}.
$$
For example the spring drift differentiates each factor $R_i,R_j$, while the <Laplacian> gives $\nabla_{\mathbf R}^2(R_iR_j)=2\delta_{ij}$. The <Hookean dumbbell conformation equation> is
$$
\boxed{\dot C-KC-CK^T=\frac{4k_BT}{\zeta}I-\frac{4H}{\zeta}C.}
$$
In the row-vector convention of $\mathbf R\cdot\nabla\mathbf u$, the displayed gradient tensor is $K^T$, giving exactly the transposes in the requested equation. The dimensionless <conformation tensor> $A=HC/(k_BT)$ consequently has <upper-convected derivative> $-(A-I)/\lambda$, with <polymer relaxation time> $\lambda=\zeta/(4H)$. This normalization matters: the unscaled $C$ has equilibrium value $(k_BT/H)I$, rather than $I$.