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.