You know, looking at this again my estimate of 0.585 is too high.
Using the errata, one makes a stress die roll, and in addition to accepting result of the experimentation, can add or subtract the risk modifier again. So, we know that the +3 is the best one to take since we can use add any number from 0 to 6 to the stress die roll in order to get a discover.
So, it goes like this:
roll a 4 with a probability of 0.11. Result is 7, "no benefit" and can add +3 again for discovery.
roll a 5 with a probability of 0.1 Result is 8 "complete failure" which according to my interpretation, is literally a complete failure and you do not get a discovery.
roll a 6 with a probability of 0.11. Result is 9, "something else or story" and can add +1 again for discovery.
roll a 7 with a probability of 0.1 Result of 10 is a discovery.
roll a 8 with a probability of 0.111 Result of 11 is modified effect. Can add -1 again for discovery.
roll a 9 with a probability of 0.1 result is 12 - reroll twice, can add -2 again for discovery.
roll a 10 with a probability of 0.01 result is 13 - can add -3 again for discovery.
Now, it seems that our base P(Discovery) = 0.541
Now on this result of reroll twice, we must accept the results, which means on a result of 12 or 13 we can still get a botch of failure that will destroy the discovery. I am explicitly assuming a botch or failure nullifies any other result. I am assuming that further "discoveries" are from the normal table of discovery results, not an original discovery used for pushing the bounds of hermetic magic. I am also assuming that a result of a failure or botch stops the endless rolling.
Assuming 4 botch dice due to doing the research in no aura at all.
First we need to look at this in terms of a single roll per the two rerolls:
P( reroll on discovery chart with +3 modifier) = P(roll a 9+ on a stress die) = which is: P(9)+ P(1,6 to 10) + P(1,1,(3 to 10))+P(1,1,1 (anything)) = 0.1+0.05+0.008+0.001=0.159 = P(R)
P(failure or botch)=0.134 = P(BF)
P(not a failure, botch or reroll) = 1-0.134-0.159 = 0.707 = P(something) = P(S)
So, per reroll we are concerned with a botch or fail since that will invalidate our original discovery.
P(Botch or Fail including all rerolls) = P(BF)+P(S)*P(BF)+P(S)P(R)(P(BF)+P(S)*P(BF)+P(S)P(R)(P(BF)+P(S)*P(BF)+P(S)P(R)(P(BF)+P(S)*P(BF)+P(S)P(R)....
=P(BF)(1+P(S))/(1-P(S)P(R))=(0.134(1+0.707))/(1-0.707*0.159)=0.2577
So, we have D= 0.541 - 0.2577*0.11=0.5127.
Using my previous formula of D/(1-(nr+r)) = P(Discovery) given , n =0.2660, r=0.059 then the net P(Discovery) = 0.554
YR7 after looking at what you did, I think it sort of works, but underestimates the variability of the process. i.e. A random number of boxes with a random number of objects inside is a stochastic probability (I think) and can be trickier than at first thought.
I'll be playing around with this on a spreadsheet and see what I can get.