I've found in other threads in the mailing list that the units of current is for example (unit of charge)/(hbar/unit of energy) (https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg01100.html). Also, the local density of states has units of energy/volume (https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg00169.html​).

The units of local density of states is rather "per energy per volume" (this is what is written in the linked message) not "energy per volume". Though

My question is, what is the units of the output of the density operator? Is it energy/volume as well?

It is "per energy per volume", as is the local density of states.

If you sum the output of a kwant.operator.Density for all scattering states at a given energy and divide the result by 2pi, it will be identical (up to numerical precision) to the output of kwant.ldos. I've attached a script that illustrates this.

I ask this because intuitively I view it as the square of the wavefunction, but it gives me values larger than 1 for each site when there is only 1 mode involved (see attached picture) ​so it is not just the probability of findinge the electron at that site because this should be maximum 1. I have also noticed that the values I get in the colorbar depend on the value of my hopping (e.g. case of graphene), but overall I'm not so sure of the units.

The scattering wavefunctions are not normalized over the scattering region, so if you sum the absolute square of the wavefunction you will not obtain 1.  The lead modes are normalized such that they carry unit current, and the scattering wavefunctions are thus normalized in a way that is commensurate with this normalization of the lead modes.

In the attached script I also show that the norm of the scattering wavefunction over the scattering region is not 1.

Happy Kwanting,