PLL phase and frequency characteristics
Before I address your questions, you should understand:

a. the integral branch of the loop filter maintains a average phase increment in units of radians/sample.  It is not a frequency value, though in the right context it can be converted to a frequency value, using your sample rate as a conversion factor.

b. the total output of the loop filter is an instantaneous phase increment which is used to advance the numerically controlled oscillator's phase.

So to your questions:


  
  Why is the input phase error of the Costas loop limited in the [1,-1] region?
  


It looks like an arbitrary limit to keep the instantaneous phase increment output of the loop filter to some sane number of radians/sample.  (Correcting more than $\pm\pi$ radians/sample is pointless.)

The average phase increment maintained by the integral branch is limited to $\pm 1.0$ radians/sample as well.  In a no-signal, noise processing situation, the integrator output can grow without bound, so the limiters are needed.  The $\pm 1.0$ radians/sample is just an arbitrary limit to keep the loop near sane operating values for when a signal does come in.  A Costas loop that is applying a correction of more than 1 radian/sample is probably not anywhere near converged anyway.


  
  How do you choose the frequency region [max_freq and min_freq] of the integrator part of the loop? why do we need to limit the frequency range?
  


In your application context, which you really haven't provided, you have to decide what are reasonable limits for the average phase increment coming out of the integrator.  You must provide limits, because the integrator is unstable and its output can grow without bound, which would make the PLL useless.

For a Costas loop, the spacing in radians of your constellation points and the number of samples/symbols will dictate what are sensible limits for the maximum &amp; minimum increment of the NCO in radians for each sample processed.  Anything outside of the interval $(-\pi, \pi)$ radians/sample would be nonsensical due to phase wrapping.

For a PLL that is trying to lock on to a tone in a particular frequency range $[f_{low}, f_{high}]$ in Hz, you can use your sample rate, $F_s$, to convert those to limits in units of radians/sample, i.e. $\left[\dfrac{\pi f_{low}}{F_s/2},\dfrac{\pi f_{high}}{F_s/2}\right]$ radians/sample.


  
  Why is the output phase Φ_o[n] (in both Costas and FLL) wrapped between -2PI and 2PI?
  


Because computation of $e^{j\phi_o}$, used for setting the NCO output, is generally numerically better with smaller angle arguments in most libraries.  Phase is cyclic so it's OK to just wrap down to small magnitude angles.