You always need to multiply by the cosine term indeed (that's part of the rendering equation). Though when you do indirect diffuse using ray-tracing and thus monte-carol integration (which is the most common technique in this case), you have to divide the contribution of each sample by your PDF. This is well exampled here.
Note also that in the mentioned reference, if the PDF has terms that you also find in the rendering equations then you can optimise the code if you wish by cancelling out these terms.
Don't forget that the BRDF of a diffuse surface is ρ/π where ρ stands for the surface albedo. So we need to divide the result by π. Though in the case of the indirect diffuse component, don't forget that we should have divided the result of castRay by the PDF of the random variable, which as we showed earlier in this chapter is 1/(2π). Dividing indirectDiffuseby 1/(2π) mis the same as multiplying this value by 2π. And since the albedo is also divided by π we can simplify the code...
You have a similar situation. If you look at the PDF for the cosine sampling, then you will realise that terms can be cancelled out. Which doesn't mean they are 'not' strictly necessary. They are, they just cancel each other out which allows to optimize the code slightly (and avoid a few division, multiplication, etc.). You are more in the micro-optimisation here... which can be confusing if you try to learn the theory by just looking at optimised code (which is often not properly commented).
$ \dfrac{(cos(\theta) ... )}{PDF} = \dfrac{(cos(\theta) ... )}{\dfrac{cos(\theta)}{\pi}} = ... $