# Butterworth filter's gain formula does not agree with R's signal package

I'm trying to calculate the Butterworth filter gain. If I use the formula mentioned on Wikipedia: $$G^2(\omega) = \frac{G_0^2} {1+\left(\frac{j\omega}{j\omega_c}\right)^{2n}}$$

I don't get a matching result from calculating the gain directly from the filter's magnitude using R's signal package.

library(signal)

# Butterworth filter

# Gain formula from wikipedia
Butterworth_gain <- function(freq, cutoff_frequency, n = 1) {
1/(1+(freq/cutoff_frequency)^(2*n))

}

bf <- butter(n = 2, W = .6, type = "low")
bfr <- freqz(bf)

plot(bfr$$f, abs(bfr$$h)^2, ty ='l')
lines(bfr$$f, Butterworth_gain(bfr$$f,  pi*.6, n = 2), col = 'red')


The function butter() computes the coefficients of a discrete-time ("digital") Butterworth filter, whereas the gain formula you used is valid for a continuous-time ("analog") Butterworth filter. According to the R documentation you can use butter() to compute an analog filter using the plane argument.
The gain of a discrete-time Butterworth filter is obtained by the bilinear transform, which substitutes the "analog" frequency variable by $$\tan(\omega/2)$$, where $$\omega$$ is the frequency in radians, normalized by the sampling frequency. For a discrete-time unit gain lowpass filter of order $$n$$, the squared gain is given by
$$G^2(\omega)=\frac{1}{1+\left(\frac{\tan(\omega/2)}{\tan(\omega_c/2)}\right)^{2n}},\qquad |\omega|\le\pi\tag{1}$$