The other answers have already covered how to generate the formatted checkpoint file (.fchk). So in my answer I will attempt to explain how to use Multiwfn to get the electron densities from the .fchk file.
First, download the software from http://sobereva.com/multiwfn/ .
Then open the MultiWfn executable, and either type in the name of the .fchk file and press enter, or press enter and then open the file through the GUI file selection interface.
If the file is loaded successfully, it should give you the message Loaded <filename> successfully!.
Then from the Main function menu, selection option 5 and press enter:
This option is for calculating parameters that have to be calculated on a 3D grid.
Then select option 1 for electron density. Then the program will ask you for the size of the grid required.
You can choose one of the inbuilt defaults (low, medium, or high quality grids). Note that grids with more points will take a longer calculation time. The paper you mentioned in the question has used a $256\times 256\times 256$ grid with a step size of $0.1$ Bohr.
If I have understood everything correctly (and I am not sure), then this is difficult to do in Multiwfn. Because in Multiwfn, you can either set the total number of grid points in each direction, or you can set the step size. This is because the program usually uses an extension distance ($6.0$ Bohr by default) to generate the grid. If you want to change both the number of grid points and the step size, then you might have to do some math to determine what would be the extension distance. You could try the options 5 to 10 and see if they work.
After setting the grid size, press enter, and the program will calculate the electron density of the file.
You can press -1 to see the isosurface graph. Select option 2 to save the electron density as a cube file in the folder (by default it saves the file as density.cub. And that's the cube file you need.
However, in the paper, they use cubman utility of Gaussian to get the difference in the electron densities between two cube files (for the same molecule calculated with different methods). Multiwfn can also do this.
It doesn't need to be said that the two cube files must have the same molecule and the exact same grid parameters, otherwise the calculation wouldn't run.
So, you need two cube files in the folder, like density1.cub and density2.cub.
Open Multiwfn and open one cube file. Select option 13 i.e. "Process grid data".
Then select 11 for "Grid data calculation".
Then you will be to see options for subtracting a grid file from the current grid file. (I believe the subtraction option is the one which is required, but I am not entirely sure about how Gaussian's cubman works, and how cube files store data, so please check if it's a simple subtract or a square subtract)
Select 4 to subtract a grid file. The program will ask you for the other grid file. (Note that which is subtracted from which depends on which file you loaded first)
Here, you have to type in the name of the other cube file, say density2.cub. Press enter, and if the calculation was successful, then you will see the message Done, grid data has been updated. If the grid sizes are different it will throw an error.
Then press 0 to save the current grid data to a cube file. Enter the name of the new file, and the difference in electron densities will be saved to that. In the paper the researchers did some further calculations from this point using their own C++ programs.
The advantage of this method is that this isn't restricted to just .fchk files from Gaussian. Multiwfn can open output files from many other programs like GAMESS, Orca etc. So, you can use the output of those programs to do this calculation. Note that the program needs the full SCF solution and the basis set data to calculate the density so those need to be in the file that you are opening.
