Summing indices during matrix multiplication

Question

Lets say I have the following matrix of probabilities that I see either or both a cat or a dog on any given day. The probability of seeing a cat or a dog are not independent but generated via some bivariate distribution.

For example where rows are indexed [did not see cat, saw cat] and columns are indexed [did not see dog, saw dog]

$$ \begin{matrix} 0.1 & 0.2 \\ 0.3 & 0.4 \\ \end{matrix} $$

(values themselves obviously are unimportant so long as they sum to 1).

Let's say I then sample from two (independent) days. If I wanted to know whether or not in total I had seen cats or dogs I would simply square the matrix. However, what I want instead is the matrix of possible combinations over two days, i.e.:

$$ \begin{matrix} (no animals) & (cat one day) & (cat both days) \\ (dog one day) & (cat one day, dog one day) & (dog one day, cat both days) \\ (dogs both days) & (dog both days, cat one day) & (both animals) \\ \end{matrix} $$

for which the manually calculated probabilities are:

$$ \begin{matrix} 0.01 & 0.04 & 0.04 \\ 0.06 & 0.20 & 0.16 \\ 0.09 & 0.24 & 0.16 \\ \end{matrix} $$

as calculated using this probably quite janky quick bit of R

mat1 <- matrix(c(0.1, 0.2, 0.3, 0.4), nrow = 2, byrow = TRUE)
mat2 <- mat1

# assume they are always square
size_mat1 <- sqrt(length(mat1))
size_mat2 <- sqrt(length(mat2))

size_mat_result <- (size_mat1 + size_mat2 - 1)
new_mat <- matrix(data = 0, nrow = size_mat_result, ncol = size_mat_result)

for(i in seq(dim(mat1)[1])){
  for(j in seq(dim(mat1)[2])){
    for(m in seq(dim(mat2)[1])){
      for(n in seq(dim(mat2)[2])){
        im_index = (i-1) + (m-1) + 1
        jn_index = (j-1) + (n-1) + 1
        
        prob_multiplication = mat1[i,j] * mat2[m,n]
        new_mat[im_index, jn_index] = new_mat[im_index, jn_index] + prob_multiplication
      }
    }
  }
}

print(new_mat)

     [,1] [,2] [,3]
[1,] 0.01 0.04 0.04
[2,] 0.06 0.20 0.16
[3,] 0.09 0.24 0.16

is there a better (more concise/ more efficient/ more general) method of calculating such a product? Googling just ends up with normal matrix multiplication.

Ideally would like to be able to iterate this some large $n$ times to (in the example above) get the matrix of total days on which cats and dogs have been seen over a month or so. Also ideally matrices for different days would not be identical and may even be different shapes (saw two cats for example). I think those extensions shouldn't prove problematic but unsure.