## 3 Results and Discussion  
For the sake of clarity, in the following we will focus on the most meaningful simulations (simulations E and F) and report the main conclusions for the simulations A, B, C and D. For simulations from A to D, we provide external references to FigShare links, which give access to the relative figures that are not reported in this paper, so as not to overload the discussion [30; 31; 32; 33; 34; 35]. Simulation A has provided indications as for the best voxel size to use during the reconstruction algorithm showing how the 1 mm and 2 mm voxel size are totally equivalent, while the 5 mm voxel size slightly worsens the localization errors and spatial resolution in favor, though, of a significant reduction of calculation time [30; 32; 35]. We used the 2 mm voxel size for the rest of simulations. Simulation B showed no significant changes with the Time Channel duration up to 48 ps. A default 12 ps time channel duration was chosen. The number of temporal steps was fixed to be 600, thus leading to a \(\sim\)7 ns long temporal scale, which fully includes the dynamic range of the simulated curves [33]. Simulation C suggested that the optimal temporal width of the software gating windows should be less than 500 ps. A wider time window introduces errors especially on the absorption coefficient estimation and reduces the penetration depth. The time window duration is given by the duration of the time scale (i.e., the number of time channels multiplied by the time channel duration) divided by the number of time windows. As we fixed the time channel duration to 12 ps and the number of temporal steps to 600, the minimum number of temporal windows is 14: A rounded value of 20 was chosen [31]. Simulation D showed how hardware gating of the detector improves the penetration depth. The relative volumetric error on the absorption coefficient, without gating, is up to 30% at 2 cm depth, while it decreases as close as 10% when several gating windows are used. We identified three as the optimal number of hardware gates: A greater number possibly reduces the photon counts in each window, which may increase error and worsen localization [34]. The previous simulations helped us to define the best variables set to provide a lower absorption quantification error and a better localization of the perturbation. Simulation **E** is of great interest as it studies the effects of the signal level, which is often a concern in optical breast imaging due to high absorption (e.g., at long wavelengths due to water concentration) and breast thickness variability. In the following we will discuss in detail the figures of merit. Each panel of Figure 3 shows the error of the reconstructed \(x\) coordinate of the inclusion with respect to the simulated one, i.e., the localization on \(x\), expressed in millimeters. The horizontal axis of each panel corresponds to the simulated \(x\) coordinate of the inclusion (from \(-\)20 mm to 20 mm at 5 mm step), while the Y axis is the depth of the inclusion, i.e., the \(z\) coordinate (from 5 mm to 25 mm in 5 mm steps). Different rows correspond to different count rates (\(10^{5}\), \(10^{6}\), \(10^{7}\), \(10^{8}\) from top to bottom) while columns represent different values of the regularization parameter (\(10^{-4}\), \(10^{-3}\), \(10^{-2}\), \(10^{-1}\), \(10^{-0}\) from left to right). The error is color-coded: Red colors correspond to positive errors, while blue colors are negative errors; small errors are identified by greenish colors. The thick black line is the error-free contour line, while a label indicates the error value on the other contour lines. Figure 3 suggests a general high localization of the \(x\) coordinate in most of the conditions as the green color (indicating low error, close to 0 mm) is spread over most of the panels except for low regularization (\(\tau=10^{-4}\)). This is expected, as a very small value of the regularization parameter leaves excessive room to the effects of the reconstruction noise.