Let's take a data-analytical approach first. We can try running the function with smaller Ns and look for a pattern.
import numpy as np
global call_cnt
def cc(arr,i,j,N):
global call_cnt
call_cnt += 1
#print(call_cnt)
if i==N or j==N or i<0 or j<0:
return
if arr[i][j]==1 or arr[i][j]=='X':
return
arr[i][j]='X'
cc(arr,i,j+1,N)
cc(arr,i,j-1,N)
cc(arr,i+1,j,N)
cc(arr,i-1,j,N)
return(arr)
n_On = []
for i in range(1,30):
global call_cnt
call_cnt = 0
cc(np.zeros((i,i)).tolist(), 0, 0, i)
print(call_cnt)
n_On.append(call_cnt)
import matplotlib.pyplot as plt
plt.plot(n_On)
plt.show()
We have initialized the arrays with non-1 values so that we don't meet the base case arr[i][j] == 1, so the program will run for the worst-case scenario.
After 29, we get a maximum recursion depth, so we just chose 1-29 for our N.
After plotting the number of total calls for each N, we get this plot.
5
17
37
65
101
145
197
257
325
401
485
577
677
785
901
1025
1157
1297
1445
1601
1765
1937
2117
2305
2501
2705
2917
3137
3365

Now, we can fit these points using many regression lines, one thing is for sure it's not exponential.
After a little bit of tuning with numbers, you can see the function 4xN^2 is very close to it, in fact, the values perfectly fit 4xN^2 + 1 function.
Let's try to plot this function on top of the number of calls.
import matplotlib.pyplot as plt
plt.plot(n_On)
plt.plot([4*(i**2) + 1 for i in range(1,30)], 'k+')
plt.xlabel('N')
plt.ylabel('number of calls')
plt.show()

It's a perfect fit!
Now, why? In your recursive function, you have mainly 2 variables that are changing i and j. As i and j are only changing in range N (due to your check i, j > 0 and i, j < N) (ixj -> NxN -> N^2) and you also did use memoization (not visiting already visited cells), so here's N^2 possibility in total.
Finally, there are 4 calls for each case. Hence, this complexity.
So, it's O(N^2).