## 2 Materials and Methods
### The DBIM-FISTA Algorithm
We have compared results from these two formulations for all our imaging scenarios and concluded that the scalar approximation produced results of similar quality. More details on this important contribution of this work are presented in Section 3. By substituting the [PERSON]'s function in (3), the ill-posed linear system is built as,  
\[\mathbf{b}=\mathbf{A}\delta_{\varepsilon}, \tag{15}\]  
where \(\mathbf{A}=\{\mathbf{A}(p,q)\}\) is an \(M\times N\) matrix (\(M\ll N\)), with \(M\) transmit-receive pairs and \(N\) voxels of the reconstruction region \(V\),  
\[\mathbf{A}(p,q)=k_{0}^{2}\Delta h^{3}\bar{\mathbf{G}}_{\mathbf{b}}(\mathbf{r_{n}},q)\bar{\mathbf{E}}_{\mathbf{b}}(q,\mathbf{r_{m}}), \tag{16}\]  
\(p\) is the \(p\)-th TR pair \((\mathbf{r_{n}},\mathbf{r_{m}})_{p}\), \(q\) is the \(q\)-th voxel inside region \(V\), and \(\mathbf{b}=\{E_{s}(p)\}\) is a \(M\times 1\) vector of the scattered fields. The computational complexity of (16) using the dyadic Green's function is \(O(M\times 9N)\), while the computational complexity further reduced to \(O(M\times 3N)\) and \(O(M\times N)\) respectively when using (13) and (14). At each iteration \(K\), the forward solver (FDTD for our implementation) is run first to obtain the background data and then used to build (15), which is solved by the inverse solver FISTA. Finally, the background profile is updated by  
\[\varepsilon_{b}^{K+1}(\mathbf{r})=\varepsilon_{b}^{K}(\mathbf{r})+\delta_{\varepsilon }(\mathbf{r}) \tag{17}\]  
and the DBIM continues to next iteration \(K+1\). We use a single-pole Debye model to simulate the frequency-dependent materials in FDTD as,  
\[\varepsilon_{r}(\omega)=\varepsilon_{\infty}+\frac{\Delta\varepsilon}{1+j \omega\tau}+\frac{\sigma_{s}}{j\omega\varepsilon_{0}}, \tag{18}\]  
where \(\varepsilon_{\infty}\), \(\Delta\varepsilon\), \(\tau\) and \(\sigma_{s}\) are the parameters of the Debye model. Thus (15) is constructed as,  
\[(\mathbf{A}^{\Re}+j\mathbf{A}^{\Im})(\delta_{\varepsilon_{\infty}}+\frac{ \delta_{\Delta\varepsilon}}{1+j\omega\tau}+\frac{\delta_{\sigma_{s}}}{j \omega\varepsilon_{0}})=\mathbf{b}^{\Re}+j\mathbf{b}^{\Im}, \tag{19}\]  
where \(\Re\) and \(\Im\) represents the real and imaginary part respectively and  
\[E_{z}^{s}=k_{0}^{2}\Delta h^{3}\sum_{r\in V}[\bar{\mathbf{G}}_{\mathbf{b}}(\mathbf{r }_{n},\mathbf{r})\bar{\mathbf{E}}_{\mathbf{b}}(\mathbf{r},\mathbf{r}_{m})\delta_{\varepsilon} (\mathbf{r})]\cdot\hat{z}, \tag{11}\]0As we only consider real valued computations, the following linear system is built as,  
\[ \begin{bmatrix}\mathbf{b}^{\Re}(1)\\ \mathbf{b}^{3}(1)\\ \mathbf{b}^{\Re}(2)\\ \cdots\\ \mathbf{b}^{\Re}(M)\\ \mathbf{b}^{\Re}(M)\end{bmatrix} = \begin{bmatrix}\bar{\mathbf{A}}_{\infty}^{ \Re}(1),\bar{\mathbf{A}}_{\infty}^{\Re}(1),\bar{\mathbf{A}}_{\infty}^{\Re}(1) \\ \bar{\mathbf{A}}_{\infty}^{\Im}(1),\bar{\mathbf{A}}_{\infty}^{\Im}(1),\bar{ \mathbf{A}}_{\infty}^{\Im}(1)\\ \bar{\mathbf{A}}_{\infty}^{\Re}(2),\bar{\mathbf{A}}_{\infty}^{\Re}(2),\bar{ \mathbf{A}}_{\infty}^{\Re}(2)\\ \bar{\mathbf{A}}_{\infty}^{\Im}(2),\bar{\mathbf{A}}_{\infty}^{\Im}(2),\bar{ \mathbf{A}}_{\infty}^{\Im}(2)\\ \cdots\\ \cdots\\ \mathbf{A}_{\infty}^{\Re}(M),\bar{\mathbf{A}}_{\infty}^{\Re}(M),\bar{ \mathbf{A}}_{\infty}^{\Re}(M)\\ \bar{\mathbf{A}}_{\infty}^{\Im}(M),\bar{\mathbf{A}}_{\infty}^{\Im}(M),\bar{ \mathbf{A}}_{\infty}^{\Im}(M)\\ \cdots\\ \delta_{\sigma_{1}}(N)\\ \delta_{\sigma_{2}}(1)\\ \delta_{\sigma_{3}}(N)\end{bmatrix}\begin{bmatrix}\delta_{\varepsilon_{\infty} }(1)\\ \delta_{\varepsilon_{\infty}}(2)\\ \ldots\\ \delta_{\varepsilon_{\infty}}(N)\\ \delta_{\Delta\varepsilon}(1)\\ \delta_{\Delta\varepsilon}(2)\\ \ldots\\ \delta_{\Delta\varepsilon}(N)\\ \delta_{\sigma_{1}}(1)\\ \delta_{\sigma_{2}}(2)\\ \ldots\\ \delta_{\sigma_{3}}(N)\end{bmatrix}