Just like voltage, current, power and so forth, that we can plot/probe it. And it can be used in optimization.  
And how can we use the harmonics in FFT analysis just as the other variables.

If you're using an optimizer around your simulation (for example, using ADE-GXL in IC610, or NeoCircuit in IC5141, or even the ADE optimizer in IC5141), you should be able to create a calculator expression which uses the dft() function and picks the value at a particular harmonic - and use that as a goal in your optimizer.

If you're hoping to do that within the simulator, that's going to be harder - although it's probably possible with spectreMDL using mvarsearch (which is effectively a simple optimizer). I've just never tried it myself.

Regards,

Andrew.