What is the first non-vanishing multipole moment of this configuration? Imagine that you have a triangle where each side has the length $a$ and a charge $q$ sitting at every vertex. Additionally, we have a charge $-3q$ sitting in the center of the triangle. What is the first non-vanishing multipole moment?
I thought that there cannot be a net charge or dipole moment, but I was uncertain about the quadrupole moment?
 A: Spherical Multipole Expansion
The first nonvanishing moment is the quadrupole moment. Assuming that you are looking for an exterior multipole expansion, the potential from the charge distribution $\rho$ can be expanded as
$$V(\mathbf{r})=\frac{1}{4\pi\epsilon_0}\sum_{L=0}^\infty\sum_{m=-L}^LI_{Lm}(\mathbf{r})\langle R_{Lm},\rho\rangle$$
where $I_{Lm},R_{Lm}$ are the real solid irregular and regular harmonics, respectively (alternately, you can use the complex solid harmonics, but since $\rho$ is real, there is no need to), and $\langle,\rangle$ denotes inner product.
From this we obtain the first several moments (in Mathematica):
SphericalHarmonicYr[L_, m_, \[Theta]_, \[Phi]_] := 
  Piecewise[{{(I*(SphericalHarmonicY[L, m, \[Theta], \[Phi]] - (-1)^m*
           SphericalHarmonicY[L, -m, \[Theta], \[Phi]]))/Sqrt[2], 
     m < 0}, {SphericalHarmonicY[L, 0, \[Theta], \[Phi]], 
     m == 0}, {(SphericalHarmonicY[L, -m, \[Theta], \[Phi]] + (-1)^m*
         SphericalHarmonicY[L, m, \[Theta], \[Phi]])/Sqrt[2], m > 0}}];
SolidHarmonicRr[l_, m_, r_, \[Theta]_, \[Phi]_] := 
  Limit[Sqrt[(4*Pi)/(2*l + 1)]*R^l*
    SphericalHarmonicYr[l, m, \[Theta], \[Phi]], R -> r];
position = {{0, 0, 0}, {a/Sqrt[3], Pi/2, 0}, {a/Sqrt[3], 
    Pi/2, (2*Pi)/3}, {a/Sqrt[3], Pi/2, (4*Pi)/3}};
charge = {-3*q, q, q, q};
Q[L_, m_] := 
  Sum[charge[[i]]*SolidHarmonicRr[L, m, ##] & @@ position[[i]], {i, 
    4}];
QQ = Simplify[Table[Q[L, m], {L, 0, 2}, {m, -L, L}]]


{{0}, {0, 0, 0}, {0, 0, -((a^2 q)/2), 0, 0}}

This means that the potential due to the configuration can be well-approximated by
$$\begin{array}{||}\hline
V(\mathbf{r})\approx -\frac{1}{4\pi\epsilon_0}\frac{qa^2}{2}I_{2,0}(\mathbf{r})=-\frac{qa^2 (3 \cos (2 \theta )+1)}{32 \pi  r^3 \epsilon _0}\\\hline
\end{array}$$
which numerically agrees with the exact potential to within $\approx20\%$ at a radius of $|\mathbf{r}|=10a$, and agrees to within $\approx1.5\%$ at $|\mathbf{r}|>100a$.
For reference, the regular real solid harmonics are given by
$$R_{Lm}(\mathbf{r})=\sqrt{\frac{4\pi}{2L+1}}r^L\begin{cases}
\sqrt{2}(-1)^m\text{Im}\left[Y_L^{|m|}(\theta,\phi)\right]&m<0\\
Y_L^0(\theta,\phi)&m=0\\
\sqrt{2}(-1)^m\text{Re}\left[Y_L^{|m|}(\theta,\phi)\right]&m>0
\end{cases}\\
=(-1)^mr^L\sqrt{\frac{(L-|m|)!}{(L+|m|)!}}P_L^{|m|}(\cos(\theta))\times\begin{cases}
\sqrt{2}\sin(|m|\phi)&m<0\\
1&m=0\\
\sqrt{2}\cos(|m|\phi)&m>0
\end{cases}$$
and the irregular harmonics $I_{Lm}$ are the same, but with $r^L$ replaced by $r^{1-L}$.
Cartesian Multipole Expansion
As vesofilev answered, the traceless Cartesian quadrupole tensor can be obtained via
$$\mathbf{Q}=\langle3\mathbf{r}\otimes\mathbf{r}-|\mathbf{r}|^2\mathbf{I},\rho\rangle$$
which in Mathematica becomes
position = 
  Prepend[Table[
    a/Sqrt[3] {Cos[2 \[Pi] k/3], Sin[2 \[Pi] k/3], 0}, {k, 3}], {0, 0,
     0}];
charge = {-3 q, q, q, q};
f[r_] := 3 r\[TensorProduct]r - (r.r) IdentityMatrix[3];
Q = Sum[charge[[k]] f[position[[k]]], {k, 4}] // MatrixForm // TeXForm

$$\mathbf{Q}=\left(
\begin{array}{ccc}
 \frac{a^2 q}{2} & 0 & 0 \\
 0 & \frac{a^2 q}{2} & 0 \\
 0 & 0 & -a^2 q \\
\end{array}
\right).$$
The potential can then be approximated as
$$V(\mathbf{r})\approx \frac{1}{4\pi\epsilon_0}\frac{\mathbf{r}^\mathsf{T}\mathbf{Q}\mathbf{r}}{2|\mathbf{r}|^5}.$$
A: It has a non-zero quadrupole moment. The easiest way to see this (Other then computing it) is to consider the expansion of the potential near a point perpendicular to the plane of the triangle and passing through the centre. If you place the origin of the coordinate system also at the origin then the radius vectors of all charges will be perpendicular to the radius vector of the point that you consider. In general you have:
$\begin{equation}V(\vec R)=\sum\limits_i\frac{q_i}{4\pi\varepsilon_0}\frac{1}{\left|\,\vec R -\vec r_i\,\right|}=\frac{\sum_iq_i}{4\pi\varepsilon_0\,R}+\frac{\sum_i q_i \vec r_i .\vec R}{4\pi\varepsilon_0\,R^3}+\frac{\sum_iq_i\left[3(\vec R.\vec r_i)^2\,-\,\vec r_i^2\vec R^2\right]}{8\pi\varepsilon_0\,R^5}+O(R^{-4})\end{equation}$
So you can easily check that in your case $\sum_iq_i=0$ and $\sum_i q_i \vec r_i .\vec R =0$. For the third therm the easiest way to check that it is non zero is to calculate it at a point when $\vec R\perp r_i$ (as I explained in the beginning). Indeed, in this case the third term is simply: $-\frac{\sum_i q_i \vec r_i^2}{8\pi\varepsilon_0\,R^3}\neq 0$ so you do have a quadrupole moment.
I hope that helps.
