## 4 Discussion  
Our methodology requires knowledge of the intercept times \(\tau_{Ud}^{p}\) and \(\tau_{Ud}^{S}\), whereas the focal depth \(z_{i}\) may be unknown. Hence, we can effectively retrieve [PERSON]'s functions at a desired (intercept) time, even in the absence of velocity information [2]. A velocity model is only required to convert intercept times into depths, akin to acoustic Marchenko imaging [45]. An important observation in this context is that the construction of a virtual \(P\)-wave source is intrinsically decoupled from the construction of a virtual \(S\)-wave source in our formalism and both tasks may even be processed independently. Consequently, we may choose \(\tau_{Ud}^{p}\) and \(\tau_{Ud}^{S}\) at mutually different focal depths without affecting the accuracy of our results. Hence, we may conclude that neither \(c^{p}(z)\) and \(c^{S}(z)\), nor the ratio \(\frac{c^{p}}{c^{S}}(z)\), is intrinsically required for the application of our methodology.  
For our numerical simulations, we have designed a medium such that the arrival times of all waveforms coincide with exact time samples (see Appendix A). In this way, we could avoid problems related to discretization. In practical applications, data are recorded within a finite frequency band only, posing limitations to our resolution, especially in the presence of thin layers [2]. It has been shown previously that some of these limitations can be overcome by enforcing energy conservation and minimum-phase conditions in the single-sided Marchenko equation [16; 42; 43]. Similar strategies may be applied to the system of equations that we presented in this paper.  
Although our methodology has been derived for a layered lossless medium with homogeneous halfspaces above \(z_{U}\) and below \(z_{L}\), it could potentially be applied to a broader range of problems. Mild lateral variations of the medium's properties may be tolerable, akin to elastodynamic Marchenko imaging of single-sided data [12; 18]. The effects of dissipation might be incorporated by computing all correlation-based reflection and transmission operators in an effectual medium, akin to the equivalent two-sided acoustic problem [33; 34]. Heterogeneities above \(z_{U}\) might be accounted for by convolving our representations with areal sources that take interactions with this part of the medium into account, akin to Rayleigh-Marchenko redatuming [46; 47]. A similar strategy may allow us to account for heterogeneities below \(z_{L}\).  
Figure 12: \(SS\)-component of the retrieved Green’s functions (organized as in Figure 9).  
The invertibility of matrix **A** in our formalism is likely to depend on the medium's properties, the available bandwidth and the wavefield components that can be emitted and recorded in practice. For transcranial applications, we may modify the theory [48] or apply redatuming [49] to account for spherical arrays. Moreover, it might be necessary to place our transducers in a water layer, which is impenetrable for \(S\)-waves. In terms of matrix algebra, such a configuration induces a projection of our multi-component wavefields to a reduced domain of \(P\)-waves only [50]. It seems plausible that such a projection would affect the invertibility of **A**, but this remains to be investigated.