Why entropy decreases in this coupled quantum systems? I have calculated the exact time evolution of a simple 1-D qubit lattice (2008 paper) and this is what I've found for $\rho(t)$ containing one excitation of 2 qubit site $(|1\rangle,|2\rangle)$ + 1 sink $(|3\rangle)$ + a vacuum state $(|0\rangle)$:

We can observe that in beginning the excitation starts in site $1$, hops into $2$ and ends in the sink (with stable population $\approx0.7$); some are dissipated (transferred to vacuum). This corresponds to the evolution of von-Neumann entropy ($k_B\equiv1$) like
$\;\;\;\;$
This is a complete system $(\mathrm{Tr}\rho(t)=1,\;t>0)$ and it is said, i.e in this answer, that the equilibrium is achieved when entropy is maximum. For my result this is clearly not the case since the entropy bumps and reaches stability at $t\rightarrow\infty$. I know the is partially because the population $\rho_{11}$ goes from $1$ to $0$, but what is really going on here? What makes this result seemingly inconsistent with the statement "entropy always increases"?
Edit:
The Hamiltonian is
\begin{equation}
H=\sum_{k=1}^N \omega_k\sigma^+_k\sigma^-_k + \sum_{k<l}\nu_{kl}(\sigma^+_k\sigma^-_l + \sigma^-_k\sigma^-_l)
\end{equation}
and the system follows Lindbladian evolution with
\begin{equation}
\mathcal{L}_{\mathrm{dissipation}}(\rho)=\sum_{k=1}^N \Gamma_k [-\{\sigma_k^+\sigma_k^-,\rho\}+2\sigma^-_k\rho\sigma^+_k], \\
\mathcal{L}_{\mathrm{dephasing}}(\rho)=\sum_{k=1}^N \gamma_k [-\{\sigma_k^+\sigma_k^-,\rho\}+2\sigma_k^+\sigma_k^-\rho\sigma_k^+\sigma_k^-], \\
\end{equation}
here site $2$ is connected to the sink, where the population cannot escape,
\begin{equation}
\mathcal{L}_{\mathrm{sink}}(\rho)=\sum_{k=1}^N \Gamma_{N+1} [-\{\sigma_2^+\sigma_{N+1}^-\sigma_{N+1}^+\sigma_{2}^-,\rho\}+2\sigma_{N+1}^+\sigma_{2}^-\rho\sigma_{k}^+\sigma_{N+1}^-]. \\
\end{equation}
I took $N=2$ and work in one-exciton manifold (denoted by sites $|1\rangle,|2\rangle$ and the sink $|3\rangle$) plus a vacuum $|0\rangle$ so that
\begin{equation}
\rho=\rho_{11}|1\rangle\langle 1| +\rho_{22}|2\rangle\langle2|+\rho_{33}|3\rangle\langle3|+\rho_{00}|0\rangle\langle 0| + \mathrm{off-diagonals},
\end{equation} 
and also we can write $\sigma^-_n=|0\rangle\langle n|$. The parameters: $\Gamma_{1,2}=0.01,\Gamma_3=0.2,\nu_{12}=0.1,\gamma=0.02$. I got the same analytical result as in the paper.
 A: There is no reason to expect that the entropy of an open system should always increase. Indeed, in many examples it clearly should decrease. Consider, for example, a two-level system undergoing spontaneous decay to the ground state. Since the final state of the two-level system is pure (i.e. with zero entropy), the entropy of the two-level system is strictly non-increasing in this case. However, the decay process is accompanied by the emission of a quantum of energy into the environment. This process typically increases the entropy of the environment to compensate for the decrease of entropy of the open system. 
The moral of the story is that only the total entropy of both system and environment should be non-decreasing if the relaxation process is supposed to describe some sort of equilibration. However, in the open-system formalism, the entropy of the environment is not accessible by construction. One can nevertheless define something like an entropy production for the environment in terms of the heat flow into it:
$$ \frac{{\rm d} S_B}{{\rm d} t} = -T^{-1}\, {\rm Tr}[\hat{H} \mathcal{L}\rho],$$
where $T$ is the temperature of the bath ($B$) described by the dissipator $\mathcal{L}$. You should be able to convince yourself that the right-hand side is proportional to the rate of change of the mean energy of the open system, which is generated by the dissipative coupling to the bath.
Note that this definition relies on the old-school thermodynamic notion of entropy changes $\delta S = \delta Q/T$, which is only correct for transitions between thermal states (assuming we want to talk about the von Neumann entropy). Hence, the entropy production as defined above is only an approximation, since the bath will not truly remain in a thermal state. Nevertheless, this approximation is typically consistent with the other assumptions underlying the master equation, which is usually valid only if the bath is weakly perturbed from equilibrium.
