Confusion about centraliser of $D_5$ in $GL_4(\mathbb{Z})$ I am trying to follow a derivation on a very old paper. My knowledge of group theory is limited, I have the basis but not much experience with advanced concepts.
We are working in 4 dimensions, so the paper quotes the 4D representations of the generators of $D_5$, $r$ (a rotation by $\pi/5$) and $p$ (the reflection $x \rightarrow x, y \rightarrow y$) as:
$$
\mathcal{R}(r) = \left (\begin{array}{ccc}
0 & 0 & 0 & -1  \\
1 & 0 & 0 & -1  \\
0 & 1 & 0 & -1  \\
0 & 0 & 1 & -1  \\
\end{array} \right ), 
$$
$$
\mathcal{R}(p) = \left (\begin{array}{ccc}
0 & 0 & 0 & 1  \\
0 & 0 & 1 & 0  \\
0 & 1 & 0 & 0  \\
1 & 0 & 0 & 0  \\
\end{array} \right ).
$$
Question 1: I thought all dihedral groups had irreducible representations of at most 2 dimensions, hence all >2 dimensional representations should be constructed from direct sums of these (so block diagonals)? Where do they get these expressions from?  
Then we calculate the centraliser of $D_5$ in $GL_4(\mathbb{Z})$, i.e. the largest subgroup of $GL_4(\mathbb{Z})$ which commutes with $D_5$. Unfortunately they skip directly to the answer because apparently "computing these groups in straightforward".  The generators of the centraliser $C$ are given by:
$$
\mathcal{R}(\delta) = \left (\begin{array}{ccc}
-1 & 1 & 0 & -1  \\
0 & 0 & 1 & -1  \\
-1 & 1 & 0 & 0  \\
-1 & 0 & 1 & -1  \\
\end{array} \right ), 
$$
$$
\mathcal{R}(\tau) = \left (\begin{array}{ccc}
-1 & 0 & 0 & 0  \\
0 & -1 & 0 & 0  \\
0 & 0 & -1 & 0  \\
0 & 0 & 0 & -1  \\
\end{array} \right ).
$$
Question 2: where do these come from? Is it there something very trivial that I'm not getting? Where would I even start?
My attempts:  
1)I started this by building a 4D representation of $D_5$ just by direct summing its 2D irreps. I am using $r_{2D}$ as the 2D rotation matrix and $p_{2D} = \left (\begin{array}{cc}
1 & 0   \\
0 & -1   \\
\end{array} \right ), 
$, to then make $
\mathcal{R}(r_{4D}) = \left (\begin{array}{cc}
r_{2D} & 0  \\
0 & r_{2D}  \\
\end{array} \right ), 
$
$
\mathcal{R}(p_{4D}) = \left (\begin{array}{ccc}
1 & 0 & 0 & 0  \\
0 & -1 & 0 & 0  \\
0 & 0 & 1 & 0  \\
0 & 0 & 0 & -1  \\
\end{array} \right ),
$
etc.
2) I used the generic element of $GL_4(\mathbb{Z})$ as $
\beta = \left (\begin{array}{ccc}
a & b & c & d  \\
e & f & g & h  \\
i & l & m & n  \\
q & s & t & u  \\
\end{array} \right ). 
$ 
Then I brute-forced computed the relationships between the entries by requiring $h$ to commute with every element $d$ of $D_5$: $h^{-1}dh = d$.  I found that the matrix above is reduced to $\beta = \left (\begin{array}{ccc}
a & 0 & c & 0  \\
0 & a & 0 & c  \\
s & 0 & u & 0  \\
0 & s & 0 & u  \\
\end{array} \right ). $ However this is where I stop, for I don't know how to compute the generator of the group of these matrices.
 A: $\def\glz{\operatorname{GL}_4(\mathbb Z)}$
$\def\inv{^{-1}}$
The question is about filling in some details about assertions made in this paper. The answer got quite long as I added more details. Guide: the part of the answer through the first edit shows that the commutant of the copy of $D_5$ in $\operatorname{GL}_4(\mathbb C)$ is two dimensional, abelian, and generated as a unital algebra by a certain matrix $M = R + R\inv$.  However, the copy of $D_5$ and the matrix $M$ both lie in $\glz$, and the subsequent part of the argument in the 4th postlude shows that $M$ and $- 1$ (minus the identity matrix) generate the centralizer of $D_5$  in $\glz$ as a group. In between, I added more detail about the first part, and also explained the construction of the 4 dimensional representation of $D_5$ given in the paper. 
I'm going to write $R$ and $J$ for your two matrix generators.  $R$ is diagonalizable with eigenvectors the 4 primitive 5th roots of 1, call them $\omega, \omega^{-1},  \omega^2, \omega^{-2}$.  In the basis of eigenvectors, matrices commuting with your $D_5$ in particular commute with $R$, so are diagonal.  Also in the basis of eigenvectors,  $J$, up to some scaling which is probably irrelevant, permutes the eigenvectors in pairs (belonging to inverse pairs of eigenvalues).  Thus your diagonal commuting matrices have only two distinct diagonal entries.  Thus the commutant is commutative, generated by two commuting idempotents of rank 2.  
Now, still in the basis of eigenvectors, $M = R + R^{-1}$ is a diagonal matrix with just two eigenvalues $\omega + \omega^{-1}$  and $\omega^2 + \omega^{-2}$  corresponding to the range of the two desired idempotents. So in fact, the whole commutant is generated by polynomials in $M$.  
Note that this is where $\mathcal R(\delta)$ came from, it is $R + R^{-1}$, in the original basis, not the basis of eigenvectors.  
I've considered the problem from the point of view of the commuting algebra, but the centralizer subgroup just consists of the invertible elements in the commuting algebra.  I've also worked over the complex numbers.   Anyway, I think this is close to what you need.
Edit:  Looking at this from more of a group representations point of view, we can get to the same place more conceptually, perhaps.  Clearly $R + R^{-1}$  is in the commutant for any representation.  So $R + R^{-1}$ acts as a scalar in each irreducible representation.  In a 2 dimensional irreducible representation, the value of that scalar is the character of $R$, which is $2 \cos(2 \pi/k)$, where $k$ can be taken as a label for the representation.  In either of the two one dimensional representations, $R + R^{-1}$  acts as the scalar $2$.   
Now in our particular representation, by examining the spectrum of $R + R^{-1}$ or of $R$, you can see that both 2 dimensional irreducible representations occur with multiplicity 1, so it follows that $R + R^{-1}$ actually generates the commutant (as an algebra).
2nd edit:  You asked how to use polynomials in $R + R^{-1}$ to generate the whole commutant.  In fact linear polynomials will do, because of the special structure here, as I will next explain.  
Mini-mini-course in representation theory:   Suppose I have a group representation $\rho : G \to \operatorname{GL}(V)$, and it decomposes as a direct sum of  mutually inequivalent irreducible subrepresentations, $V = W_1 \oplus W_2 \oplus \cdots \oplus W_n$.  Let $P_j$ be the projection operator which is $1$ on $W_j$ and zero on $W_i$  for $i \ne j$.  Thus $P_j^2 = P_j$, $P_i P_j = 0 $ for $i \ne j$ and $\sum_j P_j = 1$,
where I write $1$ for the identity operator.  Then the entire commutant $\mathcal C$ of $\rho(G)$ consists of linear combinations of $P_1,\dots, P_n$.  
Mini-mini-course in linear algebra  Now suppose I have a collection of projection operators as above $P_1,\dots, P_n$ with
$P_j^2 = P_j$, $P_i P_j = 0 $ for $i \ne j$ and $\sum_j P_j = 1$.  No matter if they came from a group representation or from somewhere else.
Let $W_j = P_j V$, so $V$ is the direct sum of the subspaces $W_j$.  Let $\mathcal C$ be the algebra generated by the projection operators $P_j$;  in fact, it consists of linear combinations of the operators $P_j$.  Suppose someone hands me a particular linear combination $T= \sum_i \alpha_i P_i$ with all the $\alpha_i$ distinct.  Then I claim that I can recover the $P_j$ and hence all of $\mathcal C$ from polynomials in $T$, and if $n = 2$ I can do it with linear polynomials.
This goes by the name interpolation formulas.  
Fix $i$ and consider
$$\prod_{j \ne i} \frac{ T - \alpha_j 1}{\alpha_i - \alpha_j},$$
which is a polynomial of degree $n-1$ in $T$.  You can check that applied to a vector in $W_k$ for $k \ne i$, it gives zero, but applied to a vector in $W_i$ it gives back the vector.  Thus the polynomial in $T$ is actually equal to $P_i$.
3rd postlude I see from your discussion/chat that you are still puzzled about the derivation of the 4d representation matrices.  As "explained" in the physics paper, the representation was derived from the 5d permutation representation exactly as Josh B. suggested.  First of all, you have to know that given a representation $V$ and a subrepresentation $W$, you get a quotient representation on $V/W$.    Start with the permutation representation of $D_5$ on $F^5$,  $F$ your favorite field,  given by
$$
R: e_0 \to e_1 \to e_2 \to e_3 \to e_4 \to e_0
$$
where $e_i$ are the standard basis elements of $F^5$  (with shifted labels)
and
$$
J:  e_0 \to e_0,\quad  e_1 \to e_4 \to e_1, \quad e_2 \to e_3 \to e_2.
$$
Now $w = \sum_{j = 0}^4 e_j$ is fixed by both $R$ and $J$. In fact, it's fixed by all permutation matrices.  Let $W = F w$ and $V = F^5/W$.
Give $V$ the basis $\overline e_1 = e_1 + W, \dots, \overline e_4 = e_4 + W$.  Note $e_0 + W = -(\overline e_1 + \overline e_2 + \overline e_3 + \overline e_4)$,  since  $e_0 = -\sum_{j = 1}^4 e_j  + w$.  It follows that the matrices of $R$ and $J$ in $V = F^5/W$ with respect to the basis
$\{\overline e_i\}_{1 \le i \le 4}$ are exactly those that you quoted from the paper.
4th postlude:  I finished it off.
Demonstration that the centralizer of $D_5$ in  $\glz$ is generated (as a group) by 
$M = R + R\inv$ and $- 1$.   Recall that $\glz$ is the set of integer matrices with determinant $\pm 1$.
We already know that any matrix that commutes with $D_5$ pointwise is a linear combination of $M$ and the identity matrix.  We need to find out when $a M + b 1$
has integer entries  and  determinant $\pm 1$.   By inspection,
 $$
 a M + b 1 =\left(
\begin{array}{cccc}
 b-a & a & 0 & -a \\
 0 & b & a & -a \\
 -a & a & b & 0 \\
 -a & 0 & a & b-a \\
\end{array}
\right)
 $$
 is integer valued if and only if $a$ and $b$ are integers.  By computation,
 $$
 \det(a M + b 1) = \left(a^2+a b-b^2\right)^2,
 $$
 and in particular $\det M = 1$.    Moreover, $\det(a M + b 1) \ge 0$, so $$a M + b 1 \in \glz  \implies \det(a M + b 1) = 1.$$ 
The minimal polynomial for $M$ is 
 $x^2 + x -1$,  as follows for example  from the computation of the eigenvalues of  $M$.  Thus we have
 $$
 M\inv = (M + 1) \quad \text{and} \quad M^2 = (1 - M).
 $$
Observation 1.  Suppose $a M + b  1 \in \glz$ then $a$ and $b$ are relatively prime.
Proof.   If $a$ and $b$ have a common prime factor $p$, then $\det(a M + b 1)$ is divisible by $p^4$.  qed
Suppose that $T  = a M + b 1 \in \glz$.  We want to show that $T$ is in the subgroup of $\glz$ generated by $M$ and $ - 1$.  If one of the coefficients  $a, b$ is zero, then $T = \pm M$ or $T = \pm 1$, and we are done,  If both of the coefficients  $a, b$ are $\pm 1$,  then $T = \pm M\inv$ or $T = \pm M^2$,  by the discussion of the minimal polynomial, so again we are done.   Now we can proceed by induction on $|a| + |b|$. 
Let us consider the case that $a b > 0$.  By extracting a factor of $ -1$  (i.e. the matrix $-1$), we can assume $a, b > 0$, and relatively prime.  Moreover, we are only interested in the case that $ a b > 1$. In particular neither of $a, b$ is divisible by the other. 
Observation 2.    $ b > a > 0$. 
Proof.  If $a >b$ then $$a ^2 + ab - b^2 > a^2 - b^2 = (a -b)(a + b) > a + b > a > 1.$$   Hence, $\det(T) >1$, a contradiction.  qed
Now we take our element $T = a M + b 1 \in \glz$ and factor out $M\inv$:
 $$
 T =  M\inv (a M^2 + b M) = M\inv ( a(1 - M) + b  M) = M\inv ( (b - a) M + a 1).
$$
The factor $ (b - a) M + a 1$ is necessarily in $\glz$, and has positive coefficients whose sum $b$ is less than the sum of the coefficients of $T$.  So our conclusion follows from the induction hypothesis.
Next consider the case that $a b < 0$.  By extracting a factor of $-1$ if necessary, we can assume
$a > 0 > b$, and $|a b| > 1$. 
Observation 3.    $ a > -b$. 
Proof.  if $-b > a$,  then 
 $$
 b^2 - ab  - a^2 > b^2 - a^2 = (b-a)(b+a) > 1.
 $$
 Hence $\det(T) > 1$, a contradiction.  qed
We take our element $T = a M + b 1 \in \glz$ and factor out $M$:
 $$
 T = M( a 1 + b M\inv)  = M ( a 1 + b (M + 1)) = M (  b M + (a + b) 1).
 $$
 The factor $b M + (a + b) 1$ is in $\glz$ and has coefficients  $a + b > 0$ and $b <0$.  The sum of absolute values of these coefficients is $(a + b) - b = a$, which is less than the corresponding sum for the coefficients of $T$, namely  $a - b$.  Again our conclusion follows from the induction hypothesis.
A: $\def\glz{\operatorname{GL}_4(\mathbb Z)}$
$\def\inv{^{-1}}$
I'm posting some further observations, as if my previous gigantic post wasn't enough.
We have specific matrices $R  = \mathcal R(r)$ and $J = \mathcal R(p)$ in $\glz$ providing a representation of the dihedral group $D_5$.   It is required to find the centralizer of $D_5$ in $\glz$
We also have the matrix $M = R  + R\inv$, denoted $\mathcal R(\delta)$ in the original post, which is an element of $\glz$ commuting with $D_5$.
The first step is to find the commutant of $D_5$, that is the set of all matrices commuting with all elements of $D_5$.
Your original brute force method for computing the commutant of $R$ works fine.  I asked Mathematica to do it, and was told the solutions are of the form
$$\left(
\begin{array}{cccc}
 a & b & 0 & -b \\
 0 & a+b & b & -b \\
 -b & b & a+b & 0 \\
 -b & 0 & b & a \\
\end{array}
\right)$$
which by examination is 
$b M + (a + b) 1$.    So this gives back the solution for the commutant, that it is the set of linear combinations of $M$ and the identity, without any of the "fancy" arguments from my previous answer.
The determinant of $ a M + b 1$  is  $D(a, b) = (a^2 + ab - b^2)^2$.   So we have a Diophantine  problem to determine the integral solutions of $D(a, b) = 1$.  
I'm not going to provide a simpler solution to this problem than already appeared in my other answer (4th postlude), but I want to point out that  Fibonacci numbers pervade this problem.  
In fact, we can first observe that the sequence of powers of $M$ is a matrix Fibonacci sequence in the following sense.   Using the minimal polynomial of $M$, namely $x^2 + x - 1$ repeatedly, we find that, with $(f_k)_{k \ge 0}$ denoting the Fibonacci sequence,
$$
(-M)^k =  -f_k M + f_{k-1}  \quad \text{and} \quad M^{-k} = f_k M + f_{k+1} 1,
$$
for $k \ge 1$.
So what we have to show is that an integer pair $(a, b)$ satisfies $D(a, b) = 1$ if and only if 
$$(a, b) =  \pm (-f_k, f_{k-1})  \quad \text{or} \quad (a, b) = \pm(f_k, f_{k+1}),
$$
for some $k$.   (I am extending the Fibonacci sequence by taking $f_0 = 0$.)
This requires more or less the same sort of work as was contained in my previous 4th postlude, and I am going to skip the details.  (This is known: Jones, James P.
Diophantine representation of the Fibonacci numbers. 
Fibonacci Quart. 13 (1975), 84–88. )
In any case the OP is interested in generalizing these results to $D_n$, $n >5$, and all of these Fibonacci phenomena are very much tied to the case of $D_5$. 
Edit:  I am adding this in response to some questions posed to me by the OP.
$\def\glz{\operatorname{GL}_4(\mathbb Z)}
\def\inv{^{-1}}
\def\boldG{\mathbf G}
\def\boldn{\mathbf n}
\def\Z{\mathbb Z}
\def\R{\mathbb R}
\def\boldq{\mathbf q}
$
Let me start over with a few things. 
The starting point in the paper is the set of points $\boldG_j  = (\cos(2 \pi  j/ 5), \sin(2 \pi  j/5))$,  $0 \le j \le 4$ in the real vector space $\mathbb R^2$.  (It is also convenient to take the indices as representing residue classes mod 5.)   These $\boldG_j$ satisfy $\sum_{j = 0}^4 \boldG_j = 0$, but
$\{\boldG_j :  1 \le j \le 4\}$ is linearly independent over $\mathbb Z$, while of course linearly dependent over $\mathbb R$.  So the set $P$ of integer linear combinations of $\{\boldG_j :  1 \le j \le 4\}$ is a complicated configuration of points in the plane,  but is also   a free $\mathbb Z$--module of rank $4$.  The dihedral group $D_5$, generated by the geometric rotation $\tau$  by angle $2 \pi/5$ in the plane and the reflection through the $x$-axis $\sigma$ leaves invariant the set of 
$\{\boldG_j :  0 \le j \le 4\}$.  The group $D_5$ acts on $\mathbb R^2$ by real linear automorphisms,  and hence acts on $P$ by $\mathbb Z$--linear automorphisms.  The matrices of the generators  $\tau$, $\sigma$  with respect to the $\Z$--basis $\{\boldG_j :  1 \le j \le 4\}$ of $P$  are the matrices $M(\tau)$ and $M(\sigma)$ given in the paper. 
We have a $\Z$--module isomorphism  $\boldq : \Z^4 \to P$, $\boldq(\boldn) = \sum_{j = 1}^4 n_j \boldG_j$. This gives us a representation of $D_5$ on $\Z^4$ with the same representation matrices.
We can regard this representation as a representation on $\mathbb R^4$ or on $\mathbb C^4$, or whatever, by using the same representation matrices.
There is also an element $\delta = \tau + \tau\inv$. (It lives in the the group ring consisting of formal linear combinations of the group elements of $D^5$.).  It acts on $P$ or on $\Z^4$ by the matrix
$M(\delta) = M(\tau) + M(\tau)\inv$.  
All of these transformations $\tau, \sigma, \delta$ will act via the same matrices, either on $P$, with respect to the basis $\{\boldG_j :  1 \le j \le 4\}$ or on $\Z^4$ with respect to the standard basis.
