Generalized Wick-Isserlis' theorem Let $n\ge 1$ be an integer and let $\vec{x} := \left( x_j \right)_{j=1}^n$ be normal variables with zero mean and with a correlation matrix ${\bf C}$. The question is to compute the following expectation value:
\begin{equation}
\mu_T(n):=E\left[ x_1 \cdot \prod\limits_{\xi=2}^n \theta_{T}(x_\xi) \right] = ?
\end{equation}
where 
\begin{equation}
\theta_T(x) := \tanh(\frac{x}{T})
\end{equation} 
is a "soft-sign" function  which switches between a Heaviside function and an idenity function when $T$ goes to zero and to infinity respectively. We compute the result for $n=2$. We have:
\begin{eqnarray}
\mu_T(2) &=& \int\limits_{{\mathbb R}^2} x_1 \theta(x_2) 
\frac{\exp(-\frac{1}{2} \vec{x}^T \cdot {\bf C}^{-1} \vec{x} )}{ \sqrt{ (2\pi)^2 \det({\bf C})}} \prod\limits_{j=1}^2 dx_j \\
&=& \int\limits_{{\mathbb R}^2} [\frac{\partial}{-\imath k_1} \delta(k_1)] \cdot [(-\frac{\imath T}{2}) \cdot \frac{1}{\sinh(k_2 \frac{T \pi}{2})}] \exp(-\frac{1}{2} \vec{k}^T \cdot {\bf C} \cdot \vec{k}) \prod\limits_{j=1}^2 dk_j \\
&=& {\bf C}_{1,2}\int\limits_{{\mathbb R}^1} [(\frac{T}{2})\frac{  k_2}{\sinh(k_2 \frac{T \pi}{2})}] \exp(-\frac{1}{2} k_2^2  )  dk_2 \\
&=& {\bf C}_{1,2} \cdot \frac{2}{T \pi^2} \int\limits_{{\mathbb R}} \frac{k_2}{\sinh(k_2)} \cdot \exp(-\frac{1}{2} \frac{ k_2^2}{\frac{\pi^2 T^2}{2^2}}) dk_2 \\
&=& {\bf C}_{1,2} \sqrt{\frac{2}{\pi}} \cdot  \int\limits_{{\mathbb R}} \frac{k_2}{\sinh(k_2)}\cdot \frac{\exp(-\frac{1}{2} \frac{k_2^2}{\sigma_T^2})}{\sqrt{2\pi \sigma_T^2}} dk_2
\end{eqnarray}
In the first line we used the definition of the expectation value. In the second line we expressed the first two terms on the left in the integrand via their Fourier transforms and we integrated over the $x$-values. In the third line we integrated over $k_1$ . In the fourth line we substituted $k_2 \leftarrow k_2 (T \pi)/2$ and in the fifith line we simplified the result and we defined $\sigma_T := \pi T/2$. Now when $T$ goes to zero the Gaussian in the integral in the fifth line tends to Dirac delta function and therefore from  we have:
\begin{equation}
\lim\limits_{T \rightarrow 0} \mu_T(2) = {\bf C}_{1,2} \sqrt{\frac{2}{\pi}}
\end{equation}
Likewise when $T$ goes to infinity the quantity in square brackets in the third line tends to a Dirac delta function divided by $T$  and therefore we have:
\begin{equation}
\lim\limits_{T \rightarrow \infty} T\mu_T(2) = {\bf C}_{1,2} 
\end{equation}
as it should be.
We also note that the quantity in question  has following series expansions.
Small-$T$ expansion:
Define:
\begin{equation}
a_m := \left. \frac{1}{m!} \frac{d^m}{d k^m} \frac{k}{\sinh(k)} \right|_{k=0}
\end{equation}
for $m=0,1,\cdots$.
Then we have
\begin{equation}
\mu_T(2) = {\bf C}_{1,2} \sqrt{\frac{2}{\pi}}\sum\limits_{m=0}^\infty a_{2 m} \sigma_T^{2 m} 2^m \frac{(m-1/2)!}{(-1/2)!}
\end{equation}
where $|\sigma_T| < 1$.
Large-$T$ expansion:
\begin{equation}
\mu_T(2) = {\bf C}_{1,2} \frac{4}{\pi} \sigma_T\sum\limits_{m=1}^\infty (-1)^{m+1} \frac{(2m-1)!!}{\sigma_T^{2 m}}\left(1-\frac{1}{2^{2m}}\right) \zeta(2 m)
\end{equation}
Now having said all this my question is what is the result for arbitrary values of $n$?
 A: This is not a full answer to this question however I want to post it because I see it as a necessary milestone in order to achieve the full answer. Let us  consider a slightly different expectation value. We define:
\begin{equation}
{\bar \mu}_T(n) := E\left[ \prod\limits_{\xi=1}^n \theta_T(x_\xi) \right]
\end{equation}
Now we take $n=2$ and compute the quantity above in a similar manner as we did in the original question. We have:
\begin{eqnarray}
{\bar \mu}_T(n) &=& \int\limits_{{\mathbb R}^2} \prod\limits_{\xi=1}^2 \theta_T(x_\xi) \cdot \frac{\exp(-\frac{1}{2} \vec{x}^T \cdot {\bf C}^{-1} \cdot \vec{x})}{\sqrt{(2\pi)^2 \det({\bf C})}} \prod\limits_{\xi=1}^2 dx_\xi \\
&=& \int\limits_{{\mathbb R}^2} \prod\limits_{\xi=1}^2 (-\frac{\imath T}{2} \frac{1}{\sinh(k_\xi \sigma_T)}) \cdot \exp\left(-\frac{1}{2} \vec{k}^T \cdot {\bf C} \cdot \vec{k} \right) \cdot \prod\limits_{\xi=1}^2 dk_\xi \\
&=& T^2 \int\limits_{{\mathbb R}_+^2}\exp\left[-\frac{1}{2} \left( k_1^2+k_2^2 \right)\right] \cdot \frac{\sinh(k_1 k_2 {\bf C}_{1,2})}{\sinh(k_1\sigma_T) \sinh(k_2 \sigma_T)} dk_1 dk_2 \\
&=& 2 \pi T^2 \int\limits_{{\mathbb R}_+^2} \prod\limits_{\xi=1}^2 \left[\frac{\exp(-\frac{1}{2} \frac{k_\xi^2}{\sigma_T^2})}{\sqrt{2\pi \sigma_T^2}}\right] \cdot \frac{\sinh\left( k_1 k_2 \frac{{\bf C}_{1,2}}{\sigma_T^2}\right)}{\sinh(k_1) \sinh(k_2)} dk_1 dk_2 \\
&=& \frac{1}{2} \pi T^2 \sum\limits_{m_1=0}^\infty \sum\limits_{m_2=0}^\infty c_{2 m_1,2 m_2}(T) \cdot \left(2 \sigma_T^2\right)^{m_1+m_2} \frac{(m_1-1/2)!}{(-1/2)!} \frac{(m_2-1/2)!}{(-1/2)!}
\end{eqnarray}
The first line is straightforward. In the second line we expressed the "soft-sign" functions through their Fourier transforms and we integrated over the $x$-values. Note that the integral on the right hand side is singular and as such it is understood as its principal value. In the third line we reduced the integral over the whole $k$-space to an integral over the first quadrant. Note that by doing so the singularity at the origin has disappeared as it should be. In the fourth line we substituted $k_\xi \rightarrow k_\xi \sigma_T$ for $\xi=1,2$.  In the last line we expanded  the ratio of hyperbolic sine functions in a Taylor series about the origin and we integrated term by term. Here we defined:
\begin{eqnarray}
c_{2m_1,2m_2}(T) &:=& \left.\frac{1}{(2m_1)! (2m_2)!} \frac{d}{d k_1^{2m_1}} \frac{d}{d k_2^{2m_2}} \left[
\frac{\sinh\left( k_1 k_2 \frac{{\bf C}_{1,2}}{\sigma_T^2}\right)}{\sinh(k_1) \sinh(k_2)}\right] \right|_{k_1=0,k_2=0}\\
&=& \sum\limits_{m_0=0}^{m_1 \wedge m_2} \frac{(\frac{{\bf C}_{1,2}}{\sigma_T^2})^{2 m_0+1}}{(2 m_0+1)!} {\mathcal a}_{2 (m_1-m_0)} {\mathcal a}_{2 (m_2-m_0)}
\end{eqnarray}
where the coefficients $\left\{ {\mathcal a}_{2 m}\right\}_{m\ge 0} = \left\{1,-\frac{1}{6},\frac{7}{360},-\frac{31}{15120},\cdots \right\}$ are solutions to the following system of equations:
\begin{equation}
\sum\limits_{m_1=0}^m {\mathcal a}_{2 m_1} \cdot \frac{1}{(2(m-m_1)+1)!} = \delta_{m,0}
\end{equation}
We have checked numerically that the series converge if $T < \sqrt{2}/\pi$.
A: Here we provide an answer for $n=4$. We will be performing exactly the same steps as we did in the main body  of the question. We have:
\begin{eqnarray}
&&\mu_T(4)=
\int\limits_{{\mathbb R}^4} x_1 \cdot \prod\limits_{\xi=2}^4 \theta_T(x_\xi) \cdot \frac{\exp\left[-\frac{1}{2} \vec{x}^T \cdot {\bf C}^{-1} \cdot \vec{x} \right]}{\sqrt{(2 \pi)^4 \det({\bf C})}} \prod\limits_{\xi=1}^4 dx_\xi \\
&&= \int\limits_{{\mathbb R}^4} \left[\frac{\partial \delta(k_1)}{\partial (-\imath k_1)}\right] \cdot \left[ \prod\limits_{\xi=2}^4 (\frac{-\imath T}{2}) \frac{1}{\sinh(k_\xi \sigma_T)}\right] \cdot \exp\left[-\frac{1}{2} \vec{k}^T \cdot {\bf C} \cdot \vec{k}\right] \cdot \prod\limits_{\xi=1}^4 dk_\xi\\
&&=\frac{(-1)^3}{\pi^3 \sigma_T} \int\limits_{{\mathbb R}^3} \left( \sum\limits_{j=2}^4 {\bf C}_{1,j} k_j\right) \cdot \prod\limits_{\xi=2}^4 \frac{1}{\sinh(k_\xi)} \cdot \exp\left[-\frac{1}{2} (k_2,k_3,k_4) {\bf C}^{(1,1)}\left( \begin{array}{r} k_2 \\ k_3 \\ k_4\end{array}\right) \right] \cdot \prod\limits_{\xi=2}^4 dk_\xi\\
&&= \sqrt{\frac{2^3}{\pi^3}} 8 \sigma_T^2 \int\limits_{{\mathbb R}_+^3} \left[\prod\limits_{\xi=2}^4 \frac{e^{-\frac{k_\xi^2}{2 \sigma_T^2}}}{\sqrt{2 \pi \sigma_T^2}}\right] \cdot \\
&&\mbox{$\frac{\sum\limits_{j=2}^4 k_j {\bf C}_{1,j} \cdot\left[\sinh(k_{j\oplus 1} k_{j\oplus 2} \frac{{\bf C}_{j\oplus 1,j\oplus 2}}{\sigma_T^2}) \cosh(k_j k_{j\oplus1}\frac{{\bf C}_{j,j\oplus 1}}{\sigma_T^2}) \cosh(k_j k_{j\oplus2}\frac{{\bf C}_{j,j\oplus 2}}{\sigma_T^2}) - \sinh \leftrightarrow \cosh\right]}{\prod\limits_{\xi=2}^4 \sinh(k_\xi)} \cdot \prod\limits_{\xi=2}^4 dk_\xi$}
\end{eqnarray}
The first line is just the definition of the expectation value. In the second line we "went to the $k$-space", i.e. expressed the main parts of the integrand through Fourier transforms and integrated over $x$-values. In the third line we integrated over $k_1$ using properties of the Dirac delta function and finally in the fourth line we reduced the integral over the whole space to that over the first octant. Here we used a following notation:
\begin{equation}
j\oplus n := \left\{ 
\begin{array}{r} (j+n) & \mbox{$if \quad j+n \le 4$}\\ 
(j+n)\%4 +1 & \mbox{otherwise}\end{array}\right.
\end{equation}
In order to bring this calculation to completion it woul be nice to provide some sort of small-$T$ expansion in the same way as we did this in the previous answers. Since this task is tedious even though pretty straightforward we leave this aside for the time being. Here we only note the limit when $T$ goes to zero. In this case each of the Gaussians in the square bracket in the last line tends to a delta function and since the integration is over the first octant one eighth of the value at the origin is picked out. Therefore we have:
\begin{equation}
\lim\limits_{T \rightarrow 0} \mu_T(4) = \sqrt{\frac{2^3}{\pi^3}} \cdot \sum\limits_{j=2}^4 {\bf C}_{1,j} {\bf C}_{j\oplus1,j\oplus 2}
\end{equation}
The right hand side is just the Isserlis'-theorem value multiplied by a constant factor $(\sqrt{2/\pi})^{n-1}$.
