Here is my integer programming formulation. Let binary variables $r_{ij}, s_{ij}, p_{ij}$ indicate whether the tile $(i,j)$ is reachable/a sprinkler/a plant accordingly.

The objective function is $\alpha \sum_{ij} s_{ij} -\sum_{ij} p_{ij}$. I set $\alpha$ sufficently small ($1/(mn)$), and when proving the upper bound of plant count, I set $\alpha=0$.

Since p,s,e are exclusive: $$ r_{ij}+p_{ij}+s_{ij} \leq 1 $$ it is $\leq$ because a tile could be empty but not reachable.

P needs R in neighborhood: $$ p_{ij} \leq \sum_{i',j' \in N(i,j)} r_{i'j'} $$ where $N(i,j)$ means the orthogonally neighborhood of $(i,j)$.

P needs S in 5x5 neighborhood: $$ p_{ij} \leq \sum_{i',j' \in F(i,j)} s_{i'j'} $$ where $F(i,j)$ means the 5x5 neighborhood of $(i,j)$.

With only constraints above, any empty tile could be R. Let binary variables $re_{ijd}$ ($d=U,D,L,R$ representing up, down, left, right) indicates whether tile $(i,j)$ "fetches" reachablility from its $d$ neighbor. So it is necessary for that neighbor to be R $$ re_{ijd} \leq r_{i'j'}, \quad (i',j') = step(i,j,d) $$

Every R except one needs reachability from others: $$ r_{ij} = \sum_d re_{ijd} \quad (i,j)\neq (i_0, j_0) $$

From a perspective of graph theory, it is actually a directed graph. The tiles are nodes and the neighboring relationship forms edges. We need to select a subset of nodes (with $r$) and a subset of edges (with $re$) which forms a tree, whose root is $(i_0, j_0)$, and every child node points to parent node. The above formulation only describes the graph structure, so there might be cycles in it (e.g. two adjacent tiles point to each other). Then it becomes subtour elimination problem in TSP in my opinion. The only difference is that we ought to forbid all cycles but in TSP we reserve the big cycles.

I forbid cycles with lazy constraints, whenever the solver yields an incumbent solution, it is easy to check and find cycles in it. Say I find a cycle $C=\{(i,j,d)\}$, I forbid it with $$ \sum_C re_{ijd} \leq |C|-1 $$ (I also generate the reversed cycle and forbid it too.)

That is all constraints of my model. It is able to return the same solution given by others in 2 minutes, and to prove the upper bound of plant count within 4 minutes (with $\alpha=0$).

Here is my integer programming formulation. Let binary variables $r_{ij}, s_{ij}, p_{ij}$ indicate whether the tile $(i,j)$ is reachable/a sprinkler/a plant accordingly.

The objective function is $\alpha \sum_{ij} s_{ij} -\sum_{ij} p_{ij}$. I set $\alpha$ sufficently small ($1/(mn)$), and when proving the upper bound of plant count, I set $\alpha=0$.

Since p,s,e are exclusive: $$ r_{ij}+p_{ij}+s_{ij} \leq 1 $$ it is $\leq$ because a tile could be empty but not reachable.

P needs R in neighborhood: $$ p_{ij} \leq \sum_{i',j' \in N(i,j)} r_{i'j'} $$ where $N(i,j)$ means the orthogonally neighborhood of $(i,j)$.

P needs S in 5x5 neighborhood: $$ p_{ij} \leq \sum_{i',j' \in F(i,j)} s_{i'j'} $$ where $F(i,j)$ means the 5x5 neighborhood of $(i,j)$.

With only constraints above, any empty tile could be R. Let binary variables $re_{ijd}$ ($d=U,D,L,R$ representing up, down, left, right) indicates whether tile $(i,j)$ "fetches" reachablility from its $d$ neighbor. So it is necessary for that neighbor to be R $$ re_{ijd} \leq r_{i'j'}, \quad (i',j') = step(i,j,d) $$

Every R except one needs reachability from others: $$ r_{ij} = \sum_d re_{ijd} \quad (i,j)\neq (i_0, j_0) $$

From a perspective of graph theory, it is actually a directed graph. The tiles are nodes and the neighboring relationship forms edges (arcs, to be accurate). We need to select a subset of nodes (with $r$) and a subset of edges (with $re$) which forms a tree, whose root is $(i_0, j_0)$, and every child node points to parent node. The above formulation only describes the graph structure, so there might be cycles in it (e.g. two adjacent tiles point to each other). Then it becomes subtour elimination problem in TSP in my opinion. The only difference is that we ought to forbid all cycles but in TSP we reserve the big cycles.

I forbid cycles with lazy constraints, whenever the solver yields an incumbent solution, it is easy to check and find cycles in it. Say I find a cycle $C=\{(i,j,d)\}$, I forbid it with $$ \sum_C re_{ijd} \leq |C|-1 $$ (I also generate the reversed cycle and forbid it too.)

That is all constraints of my model. It is able to return the same solution given by others in 2 minutes, and to prove the upper bound of plant count within 4 minutes (with $\alpha=0$).