{'fraction_non_alphanumeric': 0.04293346902429716, 'fraction_numerical': 0.011925963617860323, 'mean_word_length': 4.619192280200143, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Elastodynamic single-sided homogeneous Green's function representation: Theory and numerical examples

The homogeneous Green's function is the difference between an impulse response and its time-reversal. According to existing representation theorems, the homogeneous Green's function associated with source-receiver pairs inside a medium can be computed from measurements at a boundary enclosing the medium. However, in many applications such as seismic imaging, time-lapse monitoring, medical imaging, non-destructive testing, etc., media are only accessible from one side. A recent development of wave theory has provided a representation of the homogeneous Green's function in an elastic medium in terms of wavefield recordings at a single (open) boundary. Despite its single-sidedness, the elastodynamic homogeneous Green's function representation accounts for all orders of scattering inside the medium. We present the theory of the elastodynamic single-sided homogeneous Green's function representation and illustrate it with numerical examples for 2D laterally-invariant media. For propagating waves, the resulting homogeneous Green's functions match the exact ones within numerical precision, demonstrating the accuracy of the theory. In addition, we analyse the accuracy of the single-sided representation of the homogeneous Green's function for evanescent wave tunnelling.


Introduction
The homogeneous Green's function is the difference between an impulse response and its time-reversal. In the absence of losses, the wave equation is symmetric in time. Therefore, an impulse response to a source and its time-reversal satisfy the same wave equation. By subtracting the wave equations for these two responses from each other, we obtain a wave equation with a source term equal to zero, and a solution: the homogeneous Green's function.
In optics, Porter [1] used a closed-boundary representation of the homogeneous Green's function to retrieve the wavefield inside a medium. This representation has been the basis for inverse source problems [2] as well as inversescattering methods [3]. Unfortunately, in many practical situations, there is only single-sided access to the medium. When measurements are absent at a substantial part of the closed boundary, the retrieved homogeneous Green's function will suffer from significant artefacts. In particular in the presence of strong internal multiple scattering, these artefacts become more severe.
The closed-boundary representation can be modified to become an integral representation over the top and bottom boundaries of the medium if the medium has infinite horizontal extent [e.g. 4,5]. Further, a recent progress of wave theory has demonstrated that, after appropriate modification of the homogeneous Green's function representation, the integral contribution from the bottom boundary vanishes [6]. The result is a single-sided homogeneous Green's function representation. This representation correctly describes the wavefield inside the medium, including all orders of scattering, but excluding evanescent waves. The form of the single-sided representation is similar to the closed-boundary representation. However, the single-sided representation uses a so-called focusing function instead of the time-reversed Green's function [6]. For acoustic waves, the focusing function can be retrieved from a single-sided reflection response and an estimate of the direct arrival, using the Marchenko method [e.g. [7][8][9]. In the elastodynamic case, the approximate focusing function can be retrieved in a similar way [10,11]. However, an exact retrieval of the elastodynamic focusing function requires additional information about the medium [11]. In this paper, we assume that the elastodynamic focusing function is available (obtained either approximately by the Marchenko method or by direct modelling when the medium is known). The single-sided representation theorem provides the mathematical framework to place virtual sources and/or receivers inside the medium. Imaging techniques, e.g. for medical or geophysical applications with limited access to the medium, could benefit from this. Furthermore, virtual receivers inside a medium could be used for time-lapse monitoring purposes, i.e. to observe changes in a medium over time. Other potential applications could be forecasting the medium response to induced sources, or the localisation of passive sources inside a medium such as an earthquake [12].
We outline the theory of the single-sided homogeneous Green's function representation for elastodynamic waves in lossless media. Further, we evaluate the accuracy of the elastodynamic single-sided homogeneous Green's function representation numerically for 2D layered media. For the evaluation we use a directly modelled homogeneous Green's function as a reference. Eventually, we present an example in which the wavefield has wavenumber-frequency components that are evanescent only inside a thin layer between the recording boundary and the target depth. We demonstrate that these evanescent wavenumber-frequency components of the elastodynamic homogeneous Green's function are accounted for by the single-sided representation, except for small numerical errors.

Theory
This section consists of three parts. In the first part we show the definition of the decomposed matrix-vector wave equation. In the second part we define the elastodynamic homogeneous Green's function, and in the third part, we introduce the elastodynamic single-sided homogeneous Green's function representation. A detailed derivation of this representation can be found in Appendix A.

Matrix-vector wave equation for decomposed wavefields
We represent the elastodynamic wavefield using powerflux-normalised P-and S-wavefield potentials. Besides, we choose the depth direction x 3 as a preferential direction of propagation. For this reason, we decompose the wavefield in downward and upward propagating waves [13][14][15], and we define a 6 ⇥ 1 wave vector p containing decomposed wavefields and a 6 ⇥ 1 source vector s containing sources for these decomposed wavefields, The superscript + refers to downgoing waves, the superscript to upgoing waves. The wavefield potentials ± , ± , ⌥ ± represent the P, S1 and S2 waves, respectively (in cylindrical coordinates in a laterally-invariant medium, S1 and S2 waves are SV and SH waves). The decomposed source terms s ± are defined analogous to the decomposed wavefields p ± . After applying a forward Fourier transform from the space-time to the space-frequency domain, the matrix-vector wave equation for decomposed wavefields can be written as, Here, i is an imaginary unit (i 2 = 1), @ 3 denotes a spatial derivative in the x 3 direction and the operator B = B(x, !) accounts for the propagation and the mutual coupling of the decomposed fields [an expression of B can be found in 6].
The spatial coordinates are denoted by x = (x 1 , x 2 , x 3 ) T , the time is denoted by t and the frequency is denoted by !.

Homogeneous Green's function
Consider the decomposed matrix-vector wave equation with a delta source term, s = I (x x s ), where I is an identity matrix of appropriate size. The solution of this equation, is the Green's matrix (x, x s , !) containing decomposed wavefields,

Elastodynamic single-sided homogeneous Green's function representation
Consider a medium which is bounded by a reflection-free boundary @D 0 at the top (x 3 = x 3,0 ) as depicted in Fig. 1(a). Let R(x, x 0 , !) be the reflection response associated with a source for downgoing waves located at x 0 just above @D 0 , and a receiver for upgoing waves located at x on @D 0 . We define the depth level of x 0 as @D 0 0 . In the notation of decomposed wavefields, we can state that R(x, x 0 , !) = G ,+ (x, x 0 , !). The direct (downgoing) wave from x 0 to x on @D 0 is part of the decomposed Green's function G +,+ (x, x 0 , !) = I (x H x 0 H ). Here, the subscript H refers to the horizontal components, i.e. x H = (x 1 , x 2 ) T . Since the medium is reflection-free above @D 0 , the decomposed Green's functions associated with sources for upgoing waves at x 0 and receivers at x on @D 0 are zero, G ±, (x, x 0 , !) = O. According to the matrix notation in Eq. (12), we can write the homogeneous Green's matrix h (x, x 0 , !) for x at @D 0 in terms of the reflection response R(x, x 0 , !) and an identity matrix I of appropriate size, Next, we introduce the focusing function F(x, x s , !). The focusing function is defined in a so-called truncated medium which is identical to the true medium for x 3,0  x 3 < x 3,s and homogeneous for x 3 x 3,s , where x 3,s is the depth of the focal point at x s (see Fig. 1(b)). We assume that x s is located below @D 0 . The decomposed focusing matrix consists of a down-and an upgoing part, with, The superscripts ± and the first subscript describe the wavefield direction and the wavefield potential at the receiver location x, respectively. The second subscript describes the wavefield potential at the focusing position x s . The downgoing focusing function F + (x, x s , !) for x at @D 0 is the inverse of a transmission response of the truncated medium between @D 0 and the depth level and the complete focusing function F obeys the focusing condition, Here, we introduced the matrix I 1 = (I, O) T . The upgoing focusing function F (x, x s , !) is the reflection response of the downgoing focusing function in the truncated medium. In a physical interpretation, the focusing function is the Fourier transform of a wavefield injected at the surface @D 0 , which focuses at time zero at the position x s (see Fig. 1(b)). Note, that the first event of the focusing function is injected at negative times.
Consider the homogeneous Green's function h (x, x s , !) related to a source at x s inside the medium and receivers on the surface @D 0 as depicted in Fig. 2 from which the homogeneous Green's function h (x, x s , !) can be represented as, Evaluating the matrix products in Eqs.  is the reflection response of a medium which is reflection-free above @D 0 . The source is located at x 0 on @D 0 0 (just above @D 0 ), the receivers are located at x on @D 0 . (b) The decomposed focusing function F(x, x s , !) is defined in a truncated medium which is identical to the true medium for x 3 < x 3,s , and reflection-free for x 3 x 3,s . from which the homogeneous Green's function h (x r , x s , !) can be constructed, The quantity F(x, x r , !) is the focusing function related to a focal point at x r . Further, we introduced the matrix Eqs. (20)-(23) together form an elastodynamic single-sided homogeneous Green's function representation: It expresses the homogeneous Green's function between x s and x r inside the medium in terms of the single-sided data h (x, x 0 , !) at the upper boundary @D 0 [ @D 0 0 . An illustration of the homogeneous Green's function h (x r , x s , !) is displayed in Fig. 2 By evaluating the matrix products in Eqs.

Numerical example
In this section, we show a numerical example of the elastodynamic single-sided homogeneous Green's function representation for a 2D laterally-invariant medium. Further, we investigate the accuracy of the single-sided representation for wavenumber-frequency components of the elastodynamic homogeneous Green's function that are evanescent only inside a thin layer between the recording boundary and the virtual receiver depth. From here on, we consider a 2D medium, i.e. in all expressions of Section 2 the spatial coordinate x simplifies to x = (x 1 , x 3 ) T and the horizontal coordinate x H simplifies to x 1 . Besides, we only consider P and SV waves, indicated by the subscripts and , respectively. In the provided numerical examples, we use modelled focusing functions. This allows us to analyse the properties of the single-sided representation, independent of approximations of the focusing function retrieval via the Marchenko method.

Wavenumber-frequency domain expressions
Since we consider a laterally-invariant model the required data can be modelled efficiently by wavefield extrapolation in the wavenumber-frequency domain [21,22]. However, in the theory section the single-sided homogeneous Green's function representation is formulated in the space-frequency domain. To transform the expressions to the wavenumberfrequency domain, we use the fact that all the presented expressions have a similar form, namely a product of two 3 ) and f(x 00 In laterally-invariant media, integrals of the above form can be rewritten as a spatial convolution, which corresponds to a multiplication in the wavenumber domain [23], Here, we introduced the horizontal wavenumber k 1 . Note, when we say wavenumber domain, we refer to the horizontalwavenumber domain (k 1 , x 3 ) on the receiver side (first coordinate). Wavefields in the wavenumber domain are written with a tilde on top. Expressions in the space and wavenumber domain are mutually related via the Fourier transform, We model the required input data, the reflection , in the wavenumber-frequency domain by wavefield extrapolation. Since we model all fields for a source with a horizontal space coordinate x 1 = 0, we will omit this coordinate in the following expressions. Next, we transform Eqs.
Here, we replaced the horizontal receiver coordinate x 1,r in the exponent by the horizontal offset between the receiver and the source, x 1,r x 1,s , to account for the actual horizontal position of the source x 1,s .

Results
We present a numerical example of the elastodynamic single-sided homogeneous Green's function representation. The result is compared to the exact homogeneous Green's function, which is computed by wavefield extrapolation. For a clear illustration, we choose a simple model as depicted in Fig. 3. Results for a more complex model can be found in Appendix B. Note that we use superscripts (i) to refer to the ith layer of the medium.
According to Appendix A, the single-sided homogeneous Green's function representation ignores evanescent waves, which are associated with imaginary-valued vertical-wavenumbers . Thus, all wavenumberfrequency components that satisfy the relation, should be excluded from the analysis. Here, c max is the maximum propagation velocity of the truncated medium. In the following, before displaying a wavefield in the space-time domain, we mute evanescent waves according to Eq. (29) using the maximum P-wave velocity of the truncated medium c max , apply an inverse Fourier transform according to Eq. (28) and convolve the resulting wavefield with a 30 Hz Ricker wavelet [defined by Eq. 7 in 24]. Note that, we create a virtual source as well as a virtual receiver, meaning that there are two truncated media, bounded at the bottom by x 3,s and x 3,r , respectively. The above mentioned maximum P-wave velocity c max is the overall maximum P-wave velocity of both truncated media.
The reflection responseR and the focusing functionsF ± are modelled by wavefield extrapolation using the modelling parameters shown in Table 1.  The modelled reflection responseR(k 1 , x 3,0 , x 0 3,0 , !) contains P-wave and S-wave recordings. The components associated with a P-wave source are superimposed and shown in the space-time domain in Fig. 4(a). For clearer illustration, we only show source-receiver offsets between 2000 m and +2000 m and times between 0 s and 4 s. The reflection response contains primary reflections, converted reflections and internal multiples. Mode conversions between P-and S-waves do not occur at zero-incidence. Therefore, events with a gap at x 1 = 0 m are easily identified as converted waves. Next, we define a virtual source inside the medium at x s = (0 m, 1500 m) T . Thus, we model a focusing functioñ F(k 1 , x 3,0 , x 3,s , !) with focusing depth equal to 1500 m. The up-and downgoing P-and S-wave components of the focusing function, that are associated with a P-wave focus, are superimposed and displayed in the space-time domain in Fig. 4 Since the focusing function contains both causal and acausal events it is displayed for times between 2 s and +2 s. Subsequently, we compute the single-sided representation of the homogeneous Green's function˜ h (k 1 , x 3,r , x 3,s , !) associated with virtual receivers at x 3,r = 1700 m according to Eqs. (22)- (23). We superpose the downgoing components of the homogeneous Green's function,G , and display the absolute value of the result in Fig. 5. Here, we display the four elastic components separately. Due to the symmetry of the homogeneous Green's (12)), the sum of its upgoing components, produces an identical result. Further, since the medium is horizontally-layered we only show positive wavenumbers k 1 . In Fig. 5, the amplitude of the , and components of the quantity,G , becomes unstable. This instability could be due to either the behaviour of the elastodynamic single-sided homogeneous Green's function representation for evanescent wavenumber-frequency components, or numerical instabilities, or both. Nevertheless, the analytic elastodynamic homogeneous Green's function is characterised by an exponential amplitude decay for evanescent wavenumber-frequency components. The amplitude of the component of the quantity,G , the transition from propagating to evanescent waves is defined by the P-wave velocity of the second layer c (2) p , instead of the P-wave velocity of the third layer c p . This is expected because its component only requires an S-wave focus in the third layer of the medium. Creating an S-wave focus in the third layer only allows for P-wave conversions above layer (3). As such, the highest P-wave velocity associated with an S-wave focus in layer (3) is c (2) p . In other words, the unstable behaviour of the component for |k 1 | > ! c (2) p (to the right of the dotted red line) is caused by waves that are S-waves at the source, convert to P-waves in the second layer and convert back to S-waves before reaching the receivers. For unconverted S-waves the highest propagation velocity is c  of the single-sided representation of the elastodynamic homogeneous Green's function according to, where the absolute value is evaluated element-wise. The resulting relative errorẼ(k 1 , x 3,r , x 3,s , !) is shown in Fig. 6.
, the single-sided representation of the quantity,G , is accurate within numerical precision up to wavenumbers defined by the P-wave velocity of the second layer of the medium, |k 1 The resulting elastodynamic homogeneous Green's function˜ h (k 1 , x 3,r , x 3,s , !) is a decomposed wavefield, described by up-and downgoing P-and S-waves. To obtain the full homogeneous Green's functionG h (k 1 , x 3,r , x 3,s , !), that is associated with measurable wavefield quantities, we apply wavefield composition, whereL(k 1 , x 3,r/s , !) is the composition operator. Further details about wavefield composition can be found in Wapenaar et al. [6]. The full homogeneous Green's functionG h (k 1 , x 3,r , x 3,s , !) relates force sources f and deformation sources h to traction recordings ⌧ and particle velocity recordings v,  = 1200 m s 1 , we find that the linear artefact at the focusing position at time zero has an inclination angle ↵ = 28.7 (see Fig. 7). An inspection of the time slices, e.g. for t = 0.15 s and t = 0.45 s, demonstrates that the single-sided representation of the elastodynamic homogeneous Green's function contains primary waves, converted waves, and multiply scattered waves.
We evaluate the accuracy of the elastodynamic single-sided homogeneous Green's function representation by comparing it to the exact elastodynamic homogeneous Green's function. To that end, we modelled the elastodynamic homogeneous Green's function for an actual source at x s = (0 m, 1500 m) T . The modelled elastodynamic homogeneous Green's function is subtracted from the elastodynamic homogeneous Green's function obtained from the single-sided representation. To exclude evanescent waves, the result ˜ h (k 1 , x 3,r , x 3,s , !) is element-wise multiplied by a k 1 -! filter M determined by the maximum propagation velocity of the medium. Next, we compute the normalised Frobenius norm consists of four nt ⇥ nr block-matrices. In Fig. 8, we show the resulting residual norm as a function of virtual receiver depth x 3,r . The difference plot demonstrates that, for all modelled wavenumber-frequency components of the propagating wavefield, the elastodynamic single-sided homogeneous Green's function representation is accurate within numerical precision. For the evaluation, we used the residual of the elastodynamic decomposed homogeneous Green's function ˜ h (k 1 , x 3,r , x 3,s , !) instead of its composed equivalent to exclude effects of the wavefield composition (see Eq. (31)) from the analysis. In this case, however, the wavefield composition also performs within numerical precision.

Evanescent wave tunnelling
In this section, we investigate tunnelling of evanescent waves using a new model (see Fig. 9). The new model is nearly identical to the model in Fig. 3. However, the layer between , we would expect a rapid increase of the relative errorẼ(k 1 , x 3,r , x 3,s , !) due to instabilities, similar to Fig. 6. However, Fig. 10 shows that for , the relative errorẼ(k 1 , x 3,r , x 3,s , !) ranges from about 10 15 to about 10 6 . The relative error is not within numerical precision but still very small. We interpret this effect as evanescent wave tunnelling. In the theory section, we excluded evanescent waves. However, this restriction is only required on the boundaries on which the reciprocity theorem of the correlation-type is evaluated (see Appendix A). Therefore, the theory does not exclude evanescent wave tunnelling. Besides, the theory does not make any assumption about the thickness of the tunnelled layer. In the presented numerical example, we observed that the relative error of the tunnelling experiment depends on the thickness x 3 of the tunnelled layer. This observation is expected because the amplitude of the tunnelled evanescent waves is reduced by the factor, The single-sided representation compensates for this amplitude decay. For increasing thickness of the tunnelled layer, i.e. stronger amplitude decay, the amplitude ratio of tunnelled and propagating waves becomes smaller, and the numerical errors become larger. To estimate the maximum amplitude decay in the presented example, we maximise . Thus, we choose c = c (3) p , maximise the frequency ! = 1025 ⇥ 0.511 s 1 (see Table 1) and maximise the horizontal wavenumber k 1 = !   Green's function representation (see Fig. 7) and its modelled equivalent as a function of virtual receiver depth x 3,r . shows that, for tunnelling through the high-velocity layer of 80.0 m thickness, the smallest amplitude ratio of tunnelled and propagating waves is already in the order of 10 6 .

Discussion
We presented numerical examples of the single-sided homogeneous Green's function representation in elastic laterallyinvariant media. Such 1D models are relatively simple, but they allow us to model the required fields with very high numerical accuracy via wavefield extrapolation. Nonetheless, the theory is valid for laterally varying media, and has already been tested numerically for acoustic waves [25]. In the future, we will extend the here presented numerical example to elastic laterally varying media.
Here, we modelled the focusing functions. In practice, for the acoustic situation the focusing functions can be retrieved from the reflection response combined with a smooth velocity model via the Marchenko method, which is described by Broggini et al. [26], Van der Neut et al. [27] and others. For example, Wapenaar et al. [25] presented an acoustic single-sided homogeneous Green's function representation that uses focusing functions retrieved via the Marchenko method. Da Costa Filho et al. [10] and Wapenaar [11] extended the Marchenko method to elastodynamic waves, however, it still requires more prior knowledge of the medium than in the acoustic case. In the here presented numerical example, we used numerically modelled elastodynamic focusing functions to analyse the properties of the elastodynamic singlesided homogeneous Green's function representation, independent of approximations of the focusing functions. Currently, we are developing the (elastodynamic) Marchenko method further to minimise the required amount of prior knowledge of the medium.  Fig. 9 with the thin high-velocity layer. Note that also for tunnelled waves (the region between the dashed and dotted red lines) the relative errors are very small.
According to the theory, the single-sided representation of the elastodynamic homogeneous Green's function is accurate for waves that are non-evanescent at the virtual source and receiver depth levels. In our first numerical experiment, we only consider propagating waves and confirm the theory within numerical precision. In our second experiment, we investigate wavefield components that are evanescent inside a thin layer, but propagating at the virtual source and receiver depth levels. In this case, the numerical result for the elastodynamic homogeneous Green's function using the single-sided representation performs well also for evanescent waves, except for a small relative error. The latter experiment is a tunnelling experiment, which performs better if the tunnelled layer is thin with respect to the inverse of the absolute vertical wavenumber. Although evanescent wave tunnelling appears to be possible in theory, in practice it will suffer from noise and it will not be possible to retrieve the evanescent components of the focusing function via the Marchenko method.

Conclusion
We presented a numerical example demonstrating that single-sided access to a medium suffices to correctly retrieve the non-evanescent components of an elastodynamic homogeneous Green's function with virtual sources and virtual receivers inside the medium. Despite the single-sided access to the medium, all events of the homogeneous Green's function including primaries, internal multiples and converted waves are represented correctly. Waves that are evanescent, either on the recording boundary, or at the virtual source/receiver depth level, are neglected by the single-sided representation of the elastodynamic homogeneous Green's function, and hence, can lead to artefacts. Nevertheless, the evanescent components can be filtered to suppress these artefacts. The remaining, i.e. propagating, components can be used e.g. for imaging.
In addition, we observed in a numerical experiment that the elastodynamic single-sided homogeneous Green's function representation can tunnel evanescent components through a thin layer. In theory, this tunnelling of evanescent waves is independent of the thickness of the tunnelled layer. In practice, the single-sided representation of tunnelled waves has limited numerical accuracy because it has to compensate for the exponential amplitude decay of evanescent waves, which is stronger for thicker layers. Hence, the tunnelled layer has to be sufficiently thin with respect to the inverse of the absolute vertical wavenumber of the evanescent wave. Note that we refer to tunnelling because the single-sided representation requires that at the virtual source and receiver depth levels the elastodynamic homogeneous Green's function is propagating.
To conclude, we provided a mathematical framework to create virtual wavefield measurements inside a medium that is only accessible from a single side. We foresee potential applications in the fields of imaging, time-lapse monitoring, forecasting of the medium response as well as localisation of passive sources.

A.2. Focusing function
The downgoing part of the focusing function F + (x, x r , !) is the inverse of a transmission response related to sources at the boundary @D 0 and a receiver at x r , Z @D 0 and it obeys the focusing condition, with F containing F + and F , see Eq. (16). The upgoing focusing function F (x, x r , !) is the reflection response of the downgoing focusing function in the truncated medium. In a physical interpretation the focusing function, when transformed to the time domain, is a wavefield injected at @D 0 , which focuses at time zero at the location x r . Fig. 1(b) depicts the focusing function in a cartoon. The solid arrows represent the downgoing focusing function F + . When the downgoing focusing function is sent into the truncated medium it is partially reflected, leading to the upgoing focusing function, indicated by dashed arrows in Fig. 1(b). In the absence of a coda, the upgoing focusing function would be reflected downward again. Consequently, the focusing function would not focus at x r . This scenario is prevented by sending a coda of the downgoing focusing function into the medium to cancel the downward reflections of the upgoing focusing function. The coda is also shown in Fig. 1(b). For a 3D acoustic medium, the focusing function can be retrieved from the reflection response of the medium combined with a smooth velocity model via the Marchenko method [e.g. [7][8][9]. For an elastic medium, the focusing function retrieval still requires additional information about the medium [11].

A.3. Derivation
We insert states A and state B (Eqs. and, H(x 3 ) is the Heaviside function. We multiply Eq. (A.8) by K from the right and substitute the identities J = NK as well as JK = N, We eliminate the term with the Heaviside function by subtracting Eq. (A.9) from Eq. (A.7). The resulting expression can be written in terms of the homogeneous Green's function (using Eq. (11)), The multiplication by I T 1 N from the left in Eq. (A.10) deletes the upper submatrices of the homogeneous Green's function h (x r , x s , !). We retrieve the complete matrix h (x r , x s , !) by multiplying Eq. (A.10) by I 2 from the left, and by using Eqs. (11) and (12), representation in a similar way. First, in Eq. (A.10), we substitute x by x 0 on @D 0 0 (just above @D 0 ), x s by x (on @D 0 ) and x r by x s ,   From the reflection response and the focusing function we compute the single-sided representation of the elastodynamic homogeneous Green's function˜ h (k 1 , x 3,r , x 3,s , !) and apply a k 1 -! filter (determined by the P-wave velocity as a function of the virtual receiver depth x 3,r ). We compose the result according to Eq. (31) and obtain the full elastodynamic homogeneous Green's functionG h (k 1 , x 3,r , x 3,s , !). Next, we apply a transformation to the space-time domain and a convolution with a 30 Hz Ricker wavelet. To analyse the accuracy of the single-sided representation, we model the elastodynamic homogeneous Green's function for an actual source at x s = (0 m, 1500 m) T . We compute the difference between the modelled and the single-sided representation of the elastodynamic homogeneous Green's function. To exclude the evanescent wavefield, we element-wise multiply the residual ˜ h (k 1 , x 3,r , x 3,s , !) by a k 1 -! filterM, which is determined by the maximum propagation velocity of the medium. Subsequently, we evaluate the normalised Frobenius norm N F = 1 p 4nt4nr kM ˜ h (k 1 , x 3,r , x 3,s , !)k 2 and show the result as a function of virtual receiver depth x 3,r in Fig. B.13. The error plot demonstrates that, also in case of the 20 layer model, the single-sided homogeneous Green's function representation is accurate for propagating waves within numerical precision.