One might be able to get your second approach, with integrator, to work, however I would think adding feed-forward might be an easier solution. Namely when you reach your reference goal your system has to satisfy $C\,x_r + D\,u_r = r$ and $\dot{x}_r = A\, x_r + B\,u_r = 0$, or when writing it as one system of linear equations
$$
\begin{bmatrix}
A & B \\ C & D
\end{bmatrix}
\begin{bmatrix}
x_r \\ u_r
\end{bmatrix} =
\begin{bmatrix}
0_{12\times1} \\ r
\end{bmatrix}.
$$
Assuming that that matrix is invertible, then solving for $x_r$ and $u_r$ yields
$$
\begin{bmatrix}
x_r \\ u_r
\end{bmatrix} =
\begin{bmatrix}
A & B \\ C & D
\end{bmatrix}^{-1}
\begin{bmatrix}
0_{12\times1} \\ r
\end{bmatrix}.
$$
By defining new state space coordinates and input $\hat{x} = x - x_r$ and $\hat{u} = u - u_r$ respectively, then its dynamics can be written as
$$
\dot{\hat{x}} = \dot{x} - \dot{x}_r = A\, x + B\, u = A (\hat{x} + x_r) + B (\hat{u} + u_r) = A\, \hat{x} + B\, \hat{u} + \underbrace{A\, x_r + B\,u_r}_{0}.
$$
So $x$ can be driven to $x_r$ ($\hat{x} = 0$) by using the control law $\hat{u} = -K\, \hat{x}$ (for all $K$ such that $A-B\,K$ has all its eigenvalues in the open left half plane). This implies that the actual input $u$ should be equal to
$$
u = \hat{u} + u_r = -K (x - x_r) + u_r.
$$
When substituting in the solution for $x_r$ and $u_r$ yields
$$
u = -K\,x +
\begin{bmatrix}
K & I_{4\times4}
\end{bmatrix}
\begin{bmatrix}
A & B \\ C & D
\end{bmatrix}^{-1}
\begin{bmatrix}
0_{12\times4} \\ I_{4\times4}
\end{bmatrix} r.
$$
Since all these matrices are constant, then by defining a new matrix as
$$
F =
\begin{bmatrix}
K & I_{4\times4}
\end{bmatrix}
\begin{bmatrix}
A & B \\ C & D
\end{bmatrix}^{-1}
\begin{bmatrix}
0_{12\times4} \\ I_{4\times4}
\end{bmatrix}
$$
you can see that with your first attempt you where not that far off from this feed-forward method in terms of control law.
It can be noted that this method assumes that $r$ is constant, however when the reference changes over time, then your tracking performance might not be great. And after some playing with this myself it even seemed that choosing a $K$, such that the real parts of the poles of $A-B\,K$ are more negative, negatively affected the tracking performance a high frequencies.
It might also be worth noting that one can define a new state space model which captures the resulting reference-output dynamics as follows:
$$
\left\{
\begin{array}{l}
\dot{x} = \underbrace{(A - B\,K)}_{A_r}\,x + \underbrace{B\,F}_{B_r}\,r \\
y = \underbrace{(C - D\,K)}_{C_r}\,x + \underbrace{D\,F}_{D_r}\,r
\end{array}\right.
$$
I did a quick search online and found this document which basically contains the same derivation of the feed-forward method. It also contains the integral method, however it stops once they defined the extended state space and only suggested a control law $u = K\,x - K_I\,e$, but not how to find such $K$ and $K_I$.
I was determined to also find a solution to the integrator approach myself. But ran into difficulties using only the state space extension of $\dot{e} = r - y$. After some trial and error I decided on an extend state space of
$$
z =
\begin{bmatrix}
x \\ e \\ r \\ u
\end{bmatrix}
$$
By defining a new input to this system as $v = \dot{u}$ then the dynamics can be formulate as
$$
\dot{z} = \underbrace{
\begin{bmatrix}
A & 0 & 0 & B \\
C & 0 & -I & D \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0
\end{bmatrix}}_{\hat{A}} z + \underbrace{
\begin{bmatrix}
0 \\
0 \\
0 \\
I
\end{bmatrix}}_{\hat{B}} v.
$$
As time goes to infinity one would like $y-r$, so $\dot{e}$, to go to zero, however one would also like the system to be internally stable, this would happen if $\dot{x}$ and $\dot{u}$ would go to zero as well. So a cost function for an optimal controller could be defined as
$$
\begin{array}{rl}
J = & \int \left(\dot{e}^\top Q_e\, \dot{e} + \dot{x}^\top Q_x\, \dot{x} + \dot{u}^\top Q_u\, \dot{u}\right) dt \\
= & \int \dot{z}^\top \underbrace{
\begin{bmatrix}
Q_x & 0 & 0 & 0 \\
0 & Q_e & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & Q_u
\end{bmatrix}}_{\hat{Q}} \dot{z}\, dt
\end{array}
$$
where $Q_e$, $Q_x$ and $Q_u$ are positive definite matrices. Substituting in the expression for $\dot{z}$ gives
$$
J = \int \left(z^\top Q\, z + v^\top R\, v + 2\, z^\top N\, v\right) dt,
$$
$$
Q = \hat{A}^\top \hat{Q}\, \hat{A},
$$
$$
R = \hat{B}^\top \hat{Q}\, \hat{B} = Q_u,
$$
$$
N = \hat{A}^\top \hat{Q}\, \hat{B} = 0.
$$
This just in the standard LQR form. However I had some trouble solving this consistently, because often due to numerically errors $Q$ was not positive semi-definite and $(\hat{A},\hat{B})$ has a boundary stable mode that is not controllable. One way I was able to get this to work with some success was by adding an identity matrix times a very small positive value $\epsilon$ to $Q$ to increase all its eigenvalues by $\epsilon$ and by adding very weak damping to $r$ by setting the third matrix on the diagonal of $\hat{A}$ equal to an identity matrix times a very small negative value. If you did manage to find a controller $v = -K\,z$ with
$$
K =
\begin{bmatrix}
-K_x & -K_e & -K_r & -K_u
\end{bmatrix},
$$
then the closed-loop could also be written as
$$
\begin{bmatrix}
\dot{x} \\ \dot{e} \\ \dot{u}
\end{bmatrix} =
\begin{bmatrix}
A & 0 & B \\
C & 0 & D \\
K_x & K_e & K_u
\end{bmatrix}
\begin{bmatrix}
x \\ e \\ u
\end{bmatrix} +
\begin{bmatrix}
0 \\ -I \\ K_r
\end{bmatrix} r,
$$
$$
y =
\begin{bmatrix}
C & 0 & D
\end{bmatrix}
\begin{bmatrix}
x \\ e \\ u
\end{bmatrix}.
$$
When looking at the relevant CARE I believe it will always be the case that $K_e = 0$. If this is true then the closed-loop can be reduced even further, by removing the state $e$. Doing this yields
$$
\begin{bmatrix}
\dot{x} \\ \dot{u}
\end{bmatrix} =
\begin{bmatrix}
A & B \\
K_x & K_u
\end{bmatrix}
\begin{bmatrix}
x \\ u
\end{bmatrix} +
\begin{bmatrix}
0 \\ K_r
\end{bmatrix} r,
$$
$$
y =
\begin{bmatrix}
C & D
\end{bmatrix}
\begin{bmatrix}
x \\ u
\end{bmatrix}.
$$
After some more searching I also came across what is called Linear-Quadratic-Integral. This indeed does defines tries to find an input of the form $u = -K\,x - F\,e$ which tries to minimize the following cost function
$$
J = \int \left(\begin{bmatrix}x \\ e\end{bmatrix}^\top Q\, \begin{bmatrix}x \\ e\end{bmatrix} + u^\top R\, u + 2\,\begin{bmatrix}x \\ e\end{bmatrix}^\top N\,u\right) dt.
$$
According to the function definition of lqi in MATLAB its corresponding ARE it does indeed tries to solve the same as your first ARE
$$
\begin{bmatrix}A & 0 \\ -C & 0\end{bmatrix}^\top P + P
\begin{bmatrix}A & 0 \\ -C & 0\end{bmatrix} - \left(P
\begin{bmatrix}B & -D\end{bmatrix} + N\right) R^{-1} \left(
\begin{bmatrix}B & -D\end{bmatrix}^\top P + N^\top\right) + Q = 0.
$$
After some testing I found that this method seems to be identical to mine when you choose $Q$ for lqi the same as your $\mathcal{Q}$, so placing $Q_x$ and $Q_e$ on its diagonal, $R = Q_u$ and $N = 0$. This was somewhat unexpected since I used those same matrices to minimize $\dot{x}$, $\dot{e}$ and $\dot{u}$ instead of $x$, $e$ and $u$. The advantage of this method over mine would be that it does not seem to run into the same troubles as mine did and the corresponding ARE it tries to solve is of smaller dimension which should also make it a little faster to solve.
Best Answer
See the paper: Kalman, R. E. (1964). When is a linear control system optimal?. Journal of Basic Engineering, 86(1), 51-60.
The answer is positive at least for a class of systems. As far as I remember, the answer is also positive for a general LTI system, but I cannot find a reference at the moment.
UPDATE: Every linear system with nondynamic feedback is optimal with respect to a quadratic performance index that includes a cross-product term between the state and control, see [R1].
If you do not allow for the cross-product term, then several sufficient and necessary conditions are known, see for example [R2] and the references there.
[R1] Kreindler, E., & Jameson, A. (1972). Optimality of linear control systems. IEEE Transactions on Automatic Control, 17(3), 349-351.
[R2] Priess, M. C., Conway, R., Choi, J., Popovich, J. M., & Radcliffe, C. (2015). Solutions to the inverse lqr problem with application to biological systems analysis. IEEE Transactions on control systems technology, 23(2), 770-777.