When is separation of variables an acceptable assumption to solve a PDE?

We know that one of the classical methods for solving some PDEs is the method of separation of variables. It works for known types of PDEs and many examples of physical phenomena are successfully represented in PDE systems where an assumption that the functions are separable in variables seems to work just fine, and we get correct solutions.

My question is how do we arrive at the conclusion that the assumption that our function is 'separable in variables' is valid/correct for a particular equation in the first place? Also, I suspect that it will rely on the given scenario of boundary conditions too (not only the type of equation), such as how symmetric the boundaries are to the coordinates, whether they are curved or separable in variables themselves, etc.

Here is an example of what I mean: suppose I have the Helmholtz equation, say $\nabla^{2}f(x,y)+k^{2}f(x,y)=0$ and the boundary values are defined for some lines, like $f=0$ for $x=0$ and $x=x_{0}$ and for $y=0$ and $y=y_{0}$ (like a rectangular domain). Suppose that separation of variables worked here and we got the solution. Can we solve the same equation using the same method (i.e. separation of variables) for a new domain of boundaries, say that $f=0$ over the line segments $|y|=ax$ for $x\in[0,x_{0}]$ with $a$ being a scalar, or over a curved boundary like $y=ax^{2}$, or over the above rectangular domain of boundaries but with a cut/trim at one of its corners?

The point is, if the boundary conditions are not seperable in the given coordinate variables, wouldn't the method of separation of varaibles fail? This begs the question whether its utility is also subject to the nature of BVs, and not only to the form of the PDEs.

Thanks for any help.

I think you are asking about the possibility of satisfying the desired boundary conditions by each solution of the Helmholtz equation in the product form $X(x) Y(y)$, where $X(x)$ an $Y(y)$ were obtained by separation of variables. This can only happen if the mode functions $X(x)$ and $Y(y)$ have zero sets that are compatible with the chosen boundary. Or more specifically, if the coordinates $x$ and $y$ are compatible with a chosen boundary. For instance, using separation of variables in Cartesian coordinates to solve a boundary value problem with a circular boundary is hopeless, because neither of the $x$ or $y$ coordinates is constant even along any piece of the boundary. So, as you remarked, to use this method, one needs to choose a coordinate system that is compatible both with the differential equation and with the boundary.
You can see the complete list of 11 coordinate systems in which the Helmholtz equation is separable here on MathWorld. I think your example $|y|=ax$ with $a\ne \pm 1$ is not appropriate for any of these coordinate system. But the $y=ax^2$ example may be appropriate for parabolic coordinates.