Here is my understanding, gathered largely from the paper https://www.sciencedirect.com/science/article/abs/pii/0165168494900299
First, some non-rigorous reasoning. PCA produces an orthogonal matrix defined up to a +-1 scaling of the columns, while ICA (when applicable) yields a generic matrix defined up to any scaling and reordering of the columns. If ICA and PCA give the same matrix, then the true matrix $A$ can be obtained from it by column scaling and reordering, and thus has orthogonal columns, and this is the condition of PCA and ICA coincidence. When the sources are Gaussian, ICA simply does not work.
Now for an attempt at a more thorough argument. I am describing this in terms of random variables, not samples from them - the latter will behave similarly, but with an effect of sampling noise. (So, a 'known random vector $x$' below means that the joint distribution of its components is fully known).
ICA
Suppose $x$ and $s$ are column vectors of random variables with zero means, and $x = As$, where A is a full-rank constant matrix (possibly non-square with more rows than columns), and the components of $s$ are independent. ICA is defined as finding the matrix $A$ for a known vector $x$ (and then $s$ can also be found from the equation). It is not possible to find exactly $A$ as any matrix from the set obtained from $A$ by an arbitrary scaling of each column and an arbitrary permutation of the columns will also yield an $s$ whose components are independent. Thus, the best that can be done is finding this set of matrices. And this is possible if the components of true vector $s$ have finite, nonzero variances, and at most one of them is Gaussian (this is the theorem of ICA identifiability, see COROLLARY 13 in the paper).
PCA
Let $x$ be a column vector of random variables with zero mean and a finite covariance matrix. PCA is defined as finding an orthogonal matrix $Q$ such that for $\tilde{s}$ defined by $x = Q \tilde{s}$, the components of $\tilde{s}$ are uncorrelated and their variances form a (weakly) decreasing sequence. Again, there is an indetermination: for any such $Q$, a matrix obtained by multiplying any subset of the columns by -1 will also satisfy the conditions. If the variances of the components of $\tilde{s}$ are all different, this -1 scaling is the only indetermination, otherwise, for any set of components of $s$ with the same variance, $Q$ can be multiplied by an arbitrary rotation matrix in the space of the corresponding columns.
Now let us consider when ICA and PCA can give the same result. Let $x = As$, with $A$ being a square matrix and $s$ satisfying the conditions of ICA identifiability. Let $x = Q \tilde{s}$ be the PCA of $x$, and suppose for now that the variances of the components of $\tilde{s}$ are all different. If $Q$ satisfies the conditions of ICA then it must coincide with the true matrix $A$ up to a scaling and reordering of columns. Thus, $A$ must have orthogonal columns (which is equivalent to $A^T A$ being diagonal). Conversely, if $A$ has orthogonal columns, then by scaling them to unit norm and ordering by strictly decreasing variance of $s_i$, we obtain an orthogonal matrix and uncorrelated sources in the proper order, and thus, one of the possible PCA decompositions.
To sum up, the orthogonality of the columns of $A$ is a necessary and sufficient condition of ICA and PCA 'coincidence' under the above assumptions. The 'coincidence' should be understood as the fact that the sets of matrices obtained by both methods intersect, and each has minimal indetermination.
Concerning the question about an example where the covariance of $x$ is not diagonal, we can take $s_1$ uniform on [-1,1], $s_2$ uniform on [-2,2], and $A$ a rotation in 2D by an angle that is not an integer multiple of $\pi / 2$ (otherwise $cov(x)$ is diagonal). Here is an example for a rotation by $\pi / 6$, with 1000 samples:
The lines are: black - true axes, blue - principal component directions, red - independent component directions. PCA was computed by the MATLAB function pca, while ICA was obtained using the EEGLAB package with a call: "runica(y, 'extended', 1);" The extended algorithm is needed if some components have distributions with negative excess kurtosis (called subgaussian, or platykurtic, distributions), as is the case here for uniform distributions.
If we remove the assumption of different variances of $\tilde{s}_i$, keeping the ICA identifiability and the orthogonality of the columns of $A$, PCA will be unable to select the 'true' basis within any subspace corresponding to the same variance, while ICA, as always, will find this basis up to a reordering and rescaling. Thus, the set of matrices satisfying the PCA conditions will include the set of ICA matrices with normalized columns, but the PCA set will be substantially larger (and thus give less information about the true sources).
Now let us try to remove the assumption of ICA identifiability. As an example, consider the situation when $A$ is square and $s$ has a multivariate Gaussian distribution with a diagonal covariance matrix having different values on the diagonal. In this case, the $\tilde{s}$ obtained by PCA will also be multivariate Gaussian, and its components will be uncorrelated and thus independent. So, the set of matrices produced by PCA will be included in the set provided by ICA. However, the indetermination in ICA will be larger than usual (i.e., than in the case of ICA identifiability). Indeed, if $A$ is the true matrix, and $D$ is the (diagonal) covariance matrix of $s$ then $A D^{1/2} R$ will be an ICA matrix for any rotation R. This can be seen from the fact that $A D^{1/2}$ is an ICA matrix with such a scaling that the corresponding sources have unit covariance (and are multivariate Gaussian). Any rotation of these sources will conserve this property and thus keep them independent. So, under these assumptions, if $A$ does not have orthogonal columns, the PCA matrix set will not include the true $A$, while the ICA set will, but will have a larger than usual indetermination. Conversely, if $A$ has orthogonal columns, PCA will find $A$ with its minimal indetermination, while the ICA matrix set will contain A together with many other matrices (in practice, a random one of them will be output by software).
As a conclusion for applications, if the columns of $A$ are known to be orthogonal and $s_i$ have different variances, go with PCA (regardless of the distribution of $s_i$, but if more than one of them are Gaussian, then, out of the two methods, only PCA will recover $A$ and $s$ under these assumptions). In all other cases, ICA is the way to go.