Solving a Modified Biharmonic Equation on a Square

Question

I am seeking to solve the differential equation

\begin{equation} \left[\partial_{\overline{x}}^{4}+2\left(1+\delta\right)\partial_{\overline{x}}^{2}\partial_{\overline{y}}^{2}+\partial_{\overline{y}}^{4}\right]\psi=w^{2}\gamma\left(\partial_{\overline{x}}^{2}+\partial_{\overline{y}}^{2}\right)\psi \end{equation}

on a square geometry $(\overline{x},\overline{y})\in[-1/2,1/2]^2$. The boundary conditions on the square edges are expressed in the minimal (non-)working example of Mathematica code below:

reg = Rectangle[{-.5, -.5}, {.5, .5}];

dd[x_, a_] := UnitBox[x/a];
f[x_, a_, d_] = dd[x - d, a] - dd[x + d, a];

bi2[w_, γ_, δ_] := {
   Derivative[4, 0][ψ][x, y] + 
    2*(1 + δ)*Derivative[2, 2][ψ][x, y] + 
    Derivative[0, 4][ψ][x, y] - w^2 * γ * (Derivative[2, 0][ψ][x, y] + Derivative[0, 2][ψ][x, y])
   };

BCs[a_, d_] := {
   Derivative[1, 0][ψ][x, -.5] == f[x, a, d],
   Derivative[1, 0][ψ][x, .5] == f[x, a, d],
   Derivative[0, 1][ψ][-.5, y] == f[y, a, d],
   Derivative[0, 1][ψ][.5, y] == f[y, a, d]
   };
 
ψSol = ParametricNDSolveValue[
  {
   bi2[w2, γ2, δ2] == 0,
   BCs[a2, d2]
   },
  ψ, {x, y} ∈ reg, {a2, d2, w2, γ2, δ2},
  Method -> {"PDEDiscretization" -> {"MethodOfLines", 
      "SpatialDiscretization" -> "TensorProductGrid"}}
  ]

ψSol[.01, .1, 1., .4, .2][.1, .1]

If I try NDSolve with the method set to "FiniteElement", Mathematica complains about the degree of the PDE being greater than 2, but if I use another method (as in the above code), Mathematica complains about the boundary conditions: "Boundary values may only be specified for one independent variable."

I would greatly appreciate any advice about how to tackle this problem with Mathematica. If it helps, I don't need the full solution in the square, just values of second derivatives at the center of the square, e.g. $\psi^{\left(2,0\right)}\left(0,0\right)$, $\psi^{\left(1,1\right)}\left(0,0\right)$.

EDIT:

I've attempted a slightly different approach by rewriting the differential equation as $$ \left(\nabla^{2}\right)^{2}\psi-w^{2}\gamma\nabla^{2}\psi+2\delta\partial_{\overline{x}}^{2}\partial_{\overline{y}}^{2}\psi=0 $$ and following the advice given in the following comment: https://mathematica.stackexchange.com/a/185530/53559. Namely, by introducing auxiliary functions, we can make the PDE "look" second order. We can also use the same trick to make the boundary conditions no longer contain derivatives. However, despite these tricks my difficulties remain unresolved. Here again is a minimal (non-)working example:

reg = Rectangle[{-.5, -.5}, {.5, .5}];

dd[x_, a_] := (Pi * a)^(-1) * 1/(1 + (x/a)^2);
f[x_, a_, d_] = -dd[x + d, a] + dd[x - d, a];

bi2[w_, γ_, δ_] := {
   Derivative[2, 0][ψ][x, y] + Derivative[0, 2][ψ][x, y] - ψ2[x, y],
   Derivative[1, 1][ψ][x, y] - ψ3[x, y],
   Derivative[2, 0][ψ2][x, y] + Derivative[0, 2][ψ2][x, y] 
       + 2 * δ * Derivative[1, 1][ψ3][x, y] - w^2 * γ * ψ2[x, y],
   D[ψ[x, y], x] - ψx[x, y],
   D[ψ[x, y], y] - ψy[x, y]
};

BCs[a_, d_] := {
   ψx[x, -.5] == f[x, a, d],
   ψx[x, .5] == f[x, a, d],
   ψy[-.5, y] == f[y, a, d],
   ψy[.5, y] == f[y, a, d]
};
 
ψSol = ParametricNDSolveValue[
   {
      bi2[w2, γ2, δ2] == {0, 0, 0, 0, 0},
      BCs[a2, d2]
   },
   {ψ, ψ2, ψ3, ψx, ψy}, {x, y} ∈ reg, {a2, d2, w2, γ2, δ2}
];

ψSol[.01, .1, 1., .4, .2][[1]][.2, .1]

Note that I've also changed the boundary condition functions from box functions to Lorentzians, in response to comments to my original question.

The above code runs, and I get no complaints about the order of the PDE or the boundary conditions, but now the returned function $\psi$ is identically zero over the square. I feel like this is progress toward a working solution, but I would greatly appreciate any insight on how this might be modified to work.