I think I may have solved this myself, although my method is rather long-winded. Please correct me if I have made a mistake anywhere.
I believe the answer is $\frac{5}{4} \rho \sigma e^{\frac{1}{4} \sigma^2}$.
Here is the proof.
By definition,
\begin{equation*}
E[\vartheta U] = \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} \vartheta \chi e^{\frac{1}{2} \vartheta} \, \text{pdf}_{\chi, \vartheta}(\chi, \vartheta) \, d\chi \, d\vartheta,
\end{equation*}
where $\text{pdf}_{\chi, \vartheta}$ is the joint pdf of $\chi$ and $\vartheta$, which is well known for the bivariate normal distribution. By Fubini's theorem, we can rewrite this as
\begin{equation*}
E[\vartheta U] = \int_{-\infty}^{+\infty} \vartheta e^{\frac{1}{2} \vartheta} I_1 \, d\vartheta,
\end{equation*}
where
\begin{align*}
I_1 &:= \int_{-\infty}^{+\infty} \chi \, \text{pdf}_{\chi, \vartheta}(\chi, \vartheta) \, d\chi \\
&= \int_{\chi=-\infty}^{\chi=+\infty} \frac{1}{2\pi \sigma \sqrt{1 - \rho^2}} \chi \exp\left( -\frac{z}{2 (1 - \rho^2)} \right) \, d\chi
\end{align*}
with $z := \chi^2 - \frac{2\rho}{\sigma} \chi \vartheta + \frac{\vartheta^2}{\sigma^2}$, whence $\chi = \frac{1}{2} \frac{dz}{d\chi} + \frac{\rho}{\sigma} \vartheta$. We substitute this for $\chi$ in the integral above to obtain
\begin{align*}
I_1 &= \int_{\chi=-\infty}^{\chi=+\infty} \frac{1}{2\pi \sigma \sqrt{1 - \rho^2}} \left( \frac{1}{2} \frac{dz}{d\chi} + \frac{\rho}{\sigma} \vartheta \right) \exp\left( -\frac{z}{2 (1 - \rho^2)} \right) \, d\chi \\
&= \int_{\chi=-\infty}^{\chi=+\infty} \frac{1}{2\pi \sigma \sqrt{1 - \rho^2}} \frac{1}{2} \frac{dz}{d\chi} \exp\left( -\frac{z}{2 (1 - \rho^2)} \right) \, d\chi + \int_{\chi=-\infty}^{\chi=+\infty} \frac{1}{2\pi \sigma \sqrt{1 - \rho^2}} \frac{\rho}{\sigma} \vartheta \exp\left( -\frac{z}{2 (1 - \rho^2)} \right) \, d\chi \\
&= \frac{1}{4\pi \sigma \sqrt{1 - \rho^2}} I_2 + \frac{\rho}{\sigma} \vartheta \int_{\chi=-\infty}^{\chi=+\infty} \text{pdf}_{\chi, \vartheta}(\chi, \vartheta) \, d\chi
\end{align*}
where
\begin{equation*}
I_2 := \int_{\chi=-\infty}^{\chi=+\infty} \exp\left( -\frac{z}{2 (1 - \rho^2)} \right) \, dz.
\end{equation*}
Since $z(\chi)$ is a parabola with a minimum attained at $(\chi, z) = \left( \frac{\rho}{\sigma} \vartheta, \frac{1 - \rho^2}{\sigma^2} \vartheta^2 \right)$,
\begin{align*}
I_2 &= \int_{\chi=-\infty}^{\chi=\frac{\rho}{\sigma} \vartheta} \exp\left( -\frac{z}{2 (1 - \rho^2)} \right) \, dz + \int_{\chi=\frac{\rho}{\sigma} \vartheta}^{\chi=+\infty} \exp\left( -\frac{z}{2 (1 - \rho^2)} \right) \, dz \\
&= -\int_{z=\frac{1 - \rho^2}{\sigma^2} \vartheta^2}^{z=\frac{\rho}{\sigma} \vartheta} \exp\left( -\frac{z}{2 (1 - \rho^2)} \right) \, dz + \int_{z=\frac{1 - \rho^2}{\sigma^2} \vartheta^2}^{\chi=+\infty} \exp\left( -\frac{z}{2 (1 - \rho^2)} \right) \, dz \\
&= 0.
\end{align*}
Thus
\begin{equation*}
I_1 = \frac{\rho}{\sigma} \vartheta \int_{\chi=-\infty}^{\chi=+\infty} \text{pdf}_{\chi, \vartheta}(\chi, \vartheta) \, d\chi = \frac{\rho}{\sigma} \vartheta \, \text{pdf}_{\vartheta}(\vartheta)
\end{equation*}
and
\begin{align*}
E[\vartheta U] &= \frac{\rho}{\sigma} \int_{-\infty}^{+\infty} \vartheta^2 e^{\frac{1}{2} \vartheta} \, \text{pdf}_{\vartheta}(\vartheta) \, d\vartheta \\
&= \frac{\rho}{\sigma} \int_{-\infty}^{+\infty} \vartheta^2 e^{\frac{1}{2} \vartheta} \, \frac{1}{\sqrt{2\pi \sigma^2}} e^{\frac{-\vartheta^2}{2\sigma^2}}\, d\vartheta \\
&= \frac{\rho}{\sigma} \int_{-\infty}^{+\infty} \vartheta^2 \, \frac{1}{\sqrt{2\pi \sigma^2}} e^{\frac{-\vartheta^2 + \vartheta \sigma}{2\sigma^2}}\, d\vartheta \\
&= \frac{\rho}{\sigma} \int_{-\infty}^{+\infty} \vartheta^2 \, \frac{1}{\sqrt{2\pi \sigma^2}} e^{\frac{-(\vartheta - \sigma/2)^2 + \sigma^2/4}{2\sigma^2}}\, d\vartheta \\
&= \frac{\rho}{\sigma} e^{\frac{1}{4}\sigma^2} \int_{-\infty}^{+\infty} \vartheta^2 \, \frac{1}{\sqrt{2\pi \sigma^2}} e^{\frac{-(\vartheta - \sigma/2)^2}{2\sigma^2}}\, d\vartheta,
\end{align*}
where we recognise the integral as the second moment of a $N(\sigma/2, \sigma^2)$ random variable, which is $\frac{5}{4} \sigma^2$.