Full-Vectorial Light Propagation Simulation of Optimized Beams in Scattering Media

Volumetric scattering prevents imaging modalities in biomedical optics from imaging deep inside tissue. The optimization of the incident wavefront has the potential to improve these imaging modalities. To investigate the optimization and light propagation of such beams inside scattering media rigorously, full-vectorial simulations based on solutions of Maxwell’s equations are necessary. In this publication, we present a versatile two-step beam synthesis method to efficiently simulate the scanning and phase optimization of a focused beam inside a static scattering medium. We present four different approaches to the phase optimization of the energy density and the absolute value of the Poynting vector. We find that these quantities have two regions with different, almost exponential decays over depth for a non-optimized beam. Optimization by conjugating the phase of the projected electric field in various directions at the focus shows an improvement below a certain penetration depth. Seeking global solutions to the optimization problems reveals an even better enhancement in the energy density and the absolute value of the Poynting vector in the focus. For Poynting vector optimization, the differences between the presented optimization approaches are more significant than for the energy density. With the presented method, it is possible to efficiently simulate different imaging methods improved by wavefront shaping to investigate their possible penetration depths.


Introduction
The scattering of light is one of the greatest obstacles preventing optical imaging modalities from imaging deep inside scattering media such as biological tissue. To reduce the amount of unwanted scattered light in the resulting image, different methods have been developed, such as confocal scanning microscopy and optical coherence tomography. Most of these techniques try to effectively reject the multiply scattered light. Therefore, this has led to the sole use of ballistic light for imaging results in the limitation of the penetration depth due to the Beer-Lambert law. Another way to reduce the scattering of light is optical clearing, which modifies the optical properties of the specimen to an almost homogeneous refractive index distribution [1]. This approach is impractical in many applications. Common approaches, which also use scattered light for imaging in the optical spectrum, such as diffuse tomography [2], have a reduced resolution compared to common microscopy techniques [3].
With the advent of spatial light modulators, the manipulation of the incident light to regain a focus inside or behind a scattering medium became possible [4,5]. Spatial light modulators are high-resolution modulators that use liquid crystal display technology, deformable mirrors, or microelectromechanical systems [6][7][8]. All of these technologies have successfully been applied to correct aberrations or focus behind a scattering medium [9][10][11]. Wavefront shaping methods have since been shown to be able to improve imaging in multiple scattering media with imaging resolutions similar to state-of-the-art microscopy techniques [12]. Frequently used methods to optimize the incident wavefront are via feedback [13], phase conjugation [14], and using the transmission matrix [15,16]. To find the optimal wavefront, one needs access to the electromagnetic field at the point of interest, which is feasible outside the scattering medium but inside the medium requires, for example, a guide star [17]. If the refractive index distribution of a specimen is known, it is theoretically possible to calculate the scattering of the incident light and optimize it via simulation.
To investigate how the electromagnetic field behaves when optimizing the incident wavefront, full-vectorial electromagnetic field simulations are an invaluable tool. Common methods to model light propagation inside diffuse media based on the radiative transfer equation are not applicable in this case because they do not describe interference effects [18]. Therefore, methods to investigate light optimization by manipulating the incident wavefront must describe interference effects. Often, these types of simulations are carried out in a waveguide geometry [19] or in two-dimensional scattering media, where only out-of-plane polarization is considered [20,21]. Other simulation approaches for wavefront shaping scenarios have used approximations of wave propagation, such as the beam propagation method [22,23].
To cover all aspects of the electromagnetic field, such as polarization and phase, fullvectorial electromagnetic field simulations based on Maxwell's equations are needed, which is a computationally demanding task for a large medium. In this publication, we describe in Section 2 how to efficiently simulate the full-vectorial light propagation of focused beams inside static scattering media by a two-step beam synthesis method. This approach makes it possible to optimize the incident beam to enhance the energy density or the absolute value of the Poynting vector in the focus. Additionally, average values of the mentioned quantities can be obtained with our proposed method. Results for a focused beam scanned into a scattering medium and optimized via phase modulation are presented in Section 3. We conclude this study with Section 4, where we also discuss the main results.

Two-Step Beam Synthesis Method
We use the Debye-Wolf theory [24,25] for an aplanatic optical system to simulate the scanning and optimization of a monochromatic focused beam inside a scattering static medium. A sketch of the focusing system is depicted in Figure 1. As we are interested in the simulation of a spatial light modulator in the back focal plane of the imaging system, it is convenient to model the incident beams with their angular spectrum representation.
The incident geometrical ray with an electric field E inc (s) propagating parallel to the optical axis (z-axis) refracts at the Gaussian reference sphere, where s = (s x , s y , s z ) denotes the normalized wavevector k/k, which is fully characterized by the two coordinates s x and s y in k-space. The wavevector is denoted by k and the absolute value of the wavevector by k = 2πn m /λ with the vacuum wavelength λ and the refractive index of the surrounding medium n m . For a monochromatic beam, the third component of the normalized wavevector can be calculated by s z = (1 − s 2 x − s 2 y ) 1/2 . The refracted ray has an electric field E ∞ (s), which is the electrical field vector of a plane wave propagating in the direction s to the focus region. To obtain the electric field distribution E(r) in the focal region, we must integrate over the whole numerical aperture, each point in k-space. This results in the representation of the focusing field with the so-called Debye-Wolf integral [24] where r denotes the position vector and f is the focal length of the imaging system. To obtain E ∞ (s) from E inc (s), we follow the approach presented in [26]. The incident electric field for a specific ray is decomposed into a component that is perpendicular to the meridional plane E inc (s) · n φ and a component that is parallel to the meridional plane E inc (s) · n ρ . Figure 1. Illustration of the focusing system with the optical axis along the z direction. The incident light is characterized in k-space by the normalized coordinates s x and s y with its angular spectrum E inc (s). The in-and out-of-plane polarizations of E inc (s) with respect to the meridional plane (plane rotated by φ around the optical axis) can be calculated with the normal vectors n φ and n ρ . The aplanatic lens refracts the incident light at the Gaussian reference sphere. The resulting electric field vector E ∞ (s) for the plane wave propagating in direction s is calculated with the normal vectors n φ and n θ . The refracted light then propagates with angle θ relative to the optical axis to the focus region.
The unit vector n φ = (− sin φ, cos φ, 0) is perpendicular to the meridional plane and n ρ = (cos φ, sin φ, 0), the corresponding parallel unit vector. The angle φ describes the rotation of the meridional plane around the z-axis. The electric field components are also perpendicular to the propagation direction that is parallel to the optical axis. After the geometrical ray is refracted from the Gaussian reference sphere, the propagation direction changes to s. The perpendicular part of the electric field vector remains the same after the Gaussian reference sphere but the parallel part has to change as the electric field has to be perpendicular to s. Therefore, the electric field vector for a plane wave in the angular spectrum representation can be calculated by [26] The unit vector n θ = (cos θ cos φ, cos θ sin φ, − sin θ) is perpendicular to the propagation direction s and the factor √ s z is due to energy conservation. The angle θ describes the angle between the refracted ray direction s and the z-axis, as shown in Figure 1. We further assume that the aplanatic optical system has a polarization-independent transmission of unity. Therefore, we can write the electric field distribution at the focus with Equation (1) in terms of the angular spectrum of the incident field generated by a fully illuminated aperture as where we omit the constant factors outside the integral. For a given incident field distribution in k-space, Equation (3) can be calculated numerically to obtain the field distribution around the focus in free space. The incident field to the focus region can be used as the input for a numerical Maxwell's equations solver, such as the FDTD method [27], to calculate the scattering of a beam incident on a turbid medium. Scanning the incident beam to many positions relative to the scattering medium is computationally expensive as Maxwell's equations have to be solved for each position. In this study, we use a two-step beam synthesis method, which we already have applied for beam propagation in two-dimensional scattering media [21]. We first discretize the double integral in Equation (3) with a sum over equidistant points in k-space where ∆s x and ∆s y denote the spacing between the equidistant points in k-space. Each term in the sum corresponds to a plane wave propagating in the direction s i,j to the focus region with a polarization vector calculated by Equation (2). For a given scattering medium, Maxwell's equations are solved for each plane wave with s i,j , which results in a set of N near-field solutions. In this publication, we are interested in monochromatic beams that are scattered by a dielectric medium. Therefore, we solve Maxwell's equations on a grid with a modified Born series approach similar to the one described in [28]. Because we only examine phase modulation in this study, we restrict ourselves to incident fields with the polarization E inc = (E 0 , 0, 0) and constant amplitude E 0 in the back focal plane. The calculation of the set of N near fields is the most time-consuming part, where the bottleneck is a single plane wave simulation. Discretization of the double integral in Equation (3) leads to the repetition of the fields in the x-and y-planes. To avoid aliasing, we have to consider the Whittaker-Shannon sampling theorem [29]. Therefore, the spacing between neighboring k-vectors has to be where L i is the lateral height and width of the simulation region.
In the second step of the two-step beam synthesis method, all N near fields are added according to the angular spectrum in k-space to numerically approximate the integral in Equation (3). Thus, the plane wave part in Equation (4) is replaced by the field distribution over the simulation grid E i,j (r) to obtain the scattered electric field distribution inside the medium The summation is also done for the magnetic field distribution and results in the corresponding total magnetic field H(r) due to the linearity of Maxwell's equations. The magnetic field has to be computed in order to calculate the energy density distribution w(r) and the complex Poynting vector S(r) inside the medium.
Scanning the incident beam by ∆r only requires a phase shift for each plane wave in the angular spectrum of exp(− iks · ∆r). Therefore, the term exp − iks i,j · ∆r needs to be added in Equation (6) under the double sum. In addition, the angular spectrum can also be altered to optimize the focused incident beam to increase the energy density or the Poynting vector at a specific position. In this publication, we consider only phase modulation in the angular spectrum, although polarization and amplitude modulation are also possible with this method. For each optimization channel, a phase factor of exp(− iφ opti,i,k ) can be applied to spatially modulate the phase of the incident beam in the back focal plane. There are different ways to optimize a certain observable, such as the energy density or the Poynting vector, with the knowledge of the electromagnetic field distribution from the plane wave set but, in general, this is a global optimization problem.

Phase Optimization Techniques
The optimization of the incident electromagnetic field to obtain, for example, a focus inside a scattering medium requires either a feedback signal [13] or information about the electromagnetic field [30] at the region of interest. The total number of optimization channels N ch is usually given by the experiment setup. Often, a spatial light modulator is used, where each pixel or a group of pixels equals one optimization channel. In a fullvectorial light simulation, it is possible to access the electromagnetic field anywhere in the simulation region. Therefore, the information for optimization is easily accessible. In general, the incident light can be optimized for any objective function. Here, we focus on the energy density w(r) and the absolute value of the Poynting vector | (S(r))|. () denotes the real part. The energy density is calculated as w(r) = 1 4 ( (r)E(r) · E * (r) + µ 0 H(r) · H * (r)) and the complex Poynting vector as S(r) = 1 2 E(r) × H * (r). Each channel has an electromagnetic field denoted byẼ n (r),H n (r), which can be calculated by summing up the electromagnetic fields E i,j (r) and H i,j (r) over the corresponding region of the n-th channel in k-space.
One of the most common attempts to optimize the energy density or intensity of the electromagnetic field inside or behind a scattering medium is to use phase conjugation. This technique results in the maximum constructive interference of the electric fields, but only if light can be modeled as a scalar or all incident channels have the same polarization state in the focus position. For phase conjugation, looking only at the x component of the electric fieldẼ n,x (r foc ) at the optimization position r foc , we apply a phase change of exp(− i arg{Ẽ n,x (r foc )}) to each channel. This phase shift is also applied to the magnetic fields. The energy density and the Poynting vector are then calculated with these fields and are further denoted by "x-pc".
Scattering randomizes the phase and polarization state of light for each channel. Thus, it can be expected that phase conjugation with respect to a random direction n rand inside strongly scattering media can also enhance the energy density or absolute value of the Poynting vector. Phase conjugation along a random direction is obtained by projecting the electric field at r foc of each optimization channel along n rand and applying this phase shift of exp − i arg{Ẽ n (r foc ) · n rand } to the electromagnetic field of each channel. We further denote the optimized quantities obtained by this method as "n-rand-pc".
It is known that the electric field vector can only oscillate in a plane [31], even inside a strongly scattering medium. Therefore, the optimized electric field vector oscillates also in a plane. One option to optimize the incident field is to find the best direction n opti (φ, θ) = (cos φ sin θ, sin φ sin θ, cos θ) in which to project the electric fields of each channel to find the optimal phase shift. To find the optimal direction, we must solve the following optimization problem for the energy density w foc = w t (r foc ) at the focus r foc Here, it should be noted that we are using the electric field to project along the optimal direction n opti to find the optimal angles for φ and θ. This could also be done with the magnetic field. In the same way, we can write the optimization problem for the average value of the Poynting vector (S(r foc )) as Solving Equations (7) and (8) involves only the variables φ and θ. We further denote this optimization scheme as "n-opti".
The maximal number of variables that can be optimized is the number of channels N ch if each channel can be phase modulated. Therefore, the optimization problem for the energy density in the focus point can be written as Furthermore, the optimization problem for the average value of the Poynting vector in the focus position is In further sections, we denote quantities optimized with all the possible N ch channels as "g-opti". To find optimal solutions for Equations (7)-(10), we use the Julia packages Optim.jl [32] and Blackboxoptim.jl [33].

Scattering Medium and System Specifications
In Figure 2, the turbid medium that is used in this study is shown. The scattering medium is built up by cubic scatterers randomly positioned in a host medium with refractive index n m = 1. We choose to use cubic scatterers with the same orientation because it is easy to randomly position them without overlapping with a large concentration. Each cubic scatterer has a side length of 0.7λ and a real refractive index of 1.5. Hence, the medium is non-absorbing. The scatterers are randomly distributed inside a rectangular volume with side lengths L x = 41.7λ, L y = 41.7λ, and L z = 20λ. The volume fraction occupied inside the rectangular volume is f V = 0.3. The concentration and refractive index of the scatterers were determined by test simulations in such a way that the medium scatters relatively strongly. This medium is similar, for example, to a thin powder layer of silica particles of similar size and concentration. The refractive index distribution of the randomly positioned cubic scatterer was transferred into a simulation grid that had a resolution of λ/6. To quantify the scattering mean free path of the turbid medium, we performed a simulation, where a plane wave was scattered by the turbid medium. The propagation of the plane wave was in the z-direction and we applied periodic boundary conditions in the x-and y-direction to avoid artifacts from the finite size of the medium. Furthermore, the polarization of the incident plane wave was in the x-direction. The coherent intensity, which is often defined as E x 2 x,y (z), can be used as an estimator for [35]. The averaging of E x (r) was performed over the x-and y-direction of the simulation volume. The normalized coherent intensity E x 2 x,y (z)/ E 0 2 is shown in Figure 2b, where E 0 2 is the absolute square of the amplitude of the incident plane wave. A different approach to obtaining an estimate of is to mimic a collimated transmission experiment by performing the plane wave simulation, as described before, for different widths L z of the medium. Calculating the two-dimensional Fourier transform of E x at a transverse plane behind the medium gives the far field. The normalized intensity from the collimated transmission simulations is shown with orange markers in Figure 2b. Both approaches give similar estimates of the coherent (unscattered) intensity. A fit of an exponential function gives a scattering mean free path of = 1.05λ. The scattering mean free path is often estimated with the scattering efficiency of a single scatterer and the concentration of the scatterers [36]. If we assume a medium with the spherical scatterers but the same volume as a single scatterer, as in this study, we find that this approach would give an estimation of = 0.62λ, which is a ≈40 % relative difference from the numerical approaches. A concentration of f V = 0.1 would give a relative difference of ≈4.6 % between the different approaches. It is known that this linear dependency between the inverse of the scattering mean free path and scatterer concentration is only valid for small concentrations. Therefore, we used the decay of the coherent intensity and the collimated transmission simulations to estimate the scattering mean free path. The simulated incident focused beam is defined by the imaging system described in the previous section. In this study, we choose a numerical aperture of NA = 0.45, which results in a maximum angle of θ max = 26.74°relative to the optical axis for the allowed k-vectors. The numerical aperture limits the area in k-space as depicted in Figure 1. In order to use the two-step beam synthesis method, we have to simulate each k-vector in a separate simulation; therefore, we have to sample the aperture. We choose an equidistant spacing of ∆s x = ∆s y = 0.025 according to Equation (5) to avoid artifacts coming from the repetition of the incident beams in the lateral direction. The sampling of the aperture in normalized k-space is shown in Figure 3. In total, we simulated N = 1009 plane waves for the uniform illumination of the back focal plane and polarization of E inc = (E 0 , 0, 0) .

Results
The described vectorial two-step beam synthesis method is capable of simulating the scanning and optimization of arbitrary electromagnetic beams inside turbid media. In this section, we present results obtained with this simulation approach to scan inside a scattering medium and visualize the energy density distribution and the absolute value of the Poynting vector. Afterward, we investigate the impact of phase optimization for different numbers of optimization channels for these quantities.

Scattering of a Focused Beam by the Turbid Medium
Scanning a focused beam into a scattering medium results in the redistribution of the incident light. If the scattering mean free path is much smaller than the scanning depth, the focus completely deteriorates, which makes the focused beam unusable for imaging modalities such as confocal microscopy. To illustrate the scattering of a focused beam inside a highly scattering slab, we scanned the focus to a depth of z = 15λ inside the scattering medium. The normalized energy density distribution of the beam in vacuum w(r)/w 0 is shown in Figure 4a. A clear focus at the position r foc = (0, 0, 15λ) with w(r)/w 0 = 1 can be seen. w 0 denotes the energy density at the focus position in vacuum. The full width at half maximum of the focus in vacuum is ≈8.5λ in the axial direction and ≈1.1λ in the lateral direction. Scanning the focused beam to the same position in the presence of the scattering medium results in the energy density distribution shown in Figure 4b. The effect of scattering reduces the energy density inside the medium because a lower amount of light enters the scattering medium and also the energy is distributed over a larger volume. The maximum normalized energy density inside the scattering medium is 0.35 and its location is about a distance of λ from the entrance surface. At the original focus, the energy density drops by around two orders of magnitude to w(r foc )/w 0 = 0.003.  In the same way, the Poynting vector can be investigated. In Figure 5a, the absolute value of the Poynting vector is shown, which is normalized by the peak value in vacuum at the focus position further denoted by | (S 0 )|. The largest value can be found in the focus as expected. The full width at half maximum of the focus for | (S foc (r))|/| (S 0 )| is ≈8.5λ in the axial direction and ≈1.1λ in the lateral direction. These are the same values as for the focus size of the energy density. Introducing the scattering medium leads to the redirection of the flow of light, as can be seen in Figure 5b. Furthermore, a reduction in the intensity of the Poynting vector is additionally caused by the backscattered light, which leads to a counterpropagating flow of light, reducing the overall Poynting vector.

Phase Optimization to Enhance the Focus Energy Density
By scanning a focused beam into a scattering medium, the focus inevitably deteriorates, as shown in the previous Section 3.1. Phase optimization can lead to an increase in the energy density at the desired focus location.
In Figure 6, two simulations of phase-optimized beams are shown. The non-optimized focused light beams were first scanned to a depth of z = 15λ and afterward optimized with N ch = 137 (Figure 6a) and N ch = 1009 (Figure 6b) equally sized channels. Equally sized channels refer to the channels with the largest rectangular area in k-space that is modulated. There are additional channels with smaller areas at the edge of the circular region in k-space that are modulated, but we do not count them here as full optimization channels. For these examples, an optimal solution according to Equation (9) was sought. The phase distributions applied in k-space are shown in the left plots and the resulting energy density distribution of the optimized incident beams are shown in the right plots. After phase optimization with N ch = 137, a distinct speckle-sized focus is regained at z = 15λ. It can also be noticed that in the volume close to the incident surface, the energy density has a larger value as in the desired focus. Optimization with N ch = 1009 results also in a speckle-sized focus but with a greater value of the energy density in the focus than in any other region of the scattering medium. The reason for this is that as the number of channels increases, the electromagnetic fields constructively interfere only in the focus, and, in other regions, the field results in a random interference pattern. The axial and lateral full width at half maximum of the focus is ≈0.6λ for N ch = 137 and N ch = 1009.
To obtain the quantitative behavior of the energy density in the focus, we scanned the non-optimized and optimized beams into the medium to obtain the energy density in the focus over depth. Different optimization schemes and different numbers of optimization channels were used to obtain an enhancement in the energy density in the focus. Here, we define the enhancement as the quotient of the optimized and the non-optimized quantity. The non-optimized quantity is denoted by "scan". In Figure 7 (top), the averaged energy density w foc (z) over depth normalized by the energy density in the focus in vacuum w 0 for non-optimized ("scan") and optimized ("x-pc", "n-rand-pc", "n-opti", and "g-opti") beams is shown. The results from the optimized beams are shown with different line styles and the colors denote the different numbers of optimization channels. The averaging was performed over 49 lateral positions between the area of −6λ and 6λ in the x and y-direction, where each lateral position was separated by at least a distance of 2λ. By using only a small limited lateral area, we aimed to avoid boundary effects due to the finite size of the medium. The enhancement over depth is shown in Figure 7 (bottom). The shaded areas show the expected enhancement, where the maximum value for a specific N ch is obtained by the scalar theory [4] from the formula η(N ch ) = π/4(N ch − 1) + 1. The minimum value is obtained by η(N ch )/3. This can be justified by the consideration that the polarization of a general electromagnetic field has three spatial degrees of freedom and these orthogonal polarization states cannot interfere. For the scanned non-optimized beam, the normalized energy density in the focus first increases for shallow depths up to a distance 2λ inside the medium. This can be explained by the backscattered light from larger depths enhancing the focus energy density on average. For depths z > 2λ, w foc (z) /w 0 decreases exponentially with a factor ≈1/1.5λ, which is lower than predicted by the Beer-Lambert law. This deviation can be explained due to the presence of scattered light in the focus. The behavior of the energy density changes for depths larger than 10λ to an exponential decay with ≈1/4.8λ. This phenomenon is also present in two-dimensional, strongly scattering media [21]. Behind the medium, w foc (z) /w 0 drops abruptly due to the absence of backscattered light from larger depths.
In Figure 7 (top), we see that applying the conjugate phase of the electric field of each channel at the focus in the x-direction ("x-pc") or a random direction ("n-rand-pc") leads to no enhancement or even a decrease in the energy density in the focus for depths up to 2λ. As the electric field of the incident beam is predominantly polarized in the x-direction, conjugating the phase in this direction leads to almost no enhancement before the scattering medium and at shallow depths because all the electromagnetic fields of the channels are already in phase at the focus. For larger depths, modulating the incident beam with the two different phase conjugation approaches leads also to a significant enhancement. The results of the approaches "n-opti" and "g-opti", which utilize global optimization algorithms to find an optimal solution of the phase patterns, show similar behavior before the scattering medium and at shallow depths with a small or no enhancement compared to the nonoptimized beam. For larger depths, the enhancement is always greater for these types of optimization approaches. Furthermore, optimization according to the "g-opti" approach gives always a greater energy density in the focus than "n-opti", but one has to keep in mind that "g-opti" uses N ch variables in the optimization problems, whereas "n-opti" only uses two. The enhancement with all optimization approaches shows a steep increase until a depth of z = 8λ. This effect can be explained by the fact that the focused incident beam is itself optimized for a vacuum and deteriorates when scanned into the scattering medium. After a depth of z > 8λ, the focus is practically lost, and therefore the enhancements for different numbers of incident channels stagnate or have a decreased slope. For the energy density enhancement, only the results for N ch = 9 reach the theoretically expected enhancement between η(9)/3 = 2.4 and η(9) = 7.3. For more optimization channels, the different electromagnetic fields might not be completely randomized for the scanned non-optimized beam and the consulted scalar theory might not be fully applicable to vectorial light propagation.

Phase Optimization to Enhance the Absolute Value of the Focus Poynting Vector
In addition, it is often of interest to investigate the flow of light described by the Poynting vector and how it can be increased by the phase modulation of the incident beam. As seen in Figure 5, a scattering medium reduces the absolute value of the Poynting vector in the focus compared to the value in vacuum S 0 . Phase optimization according to Equation (10) enhances the absolute value of the Poynting vector at the desired focus. This can be seen for a scanning position at z = 15λ in Figure 8 for two different numbers of optimization channels. We find that in the case of N ch = 137, a distinct focus is regained at the desired focus position, although | (S(r foc ))|/| (S 0 ))| can have greater values in the scattering medium at shallow depths. For N ch = 1009, the maximum absolute value of the Poynting vector is in the desired focus and also larger as for the simulation with N ch = 137. In both cases, the full width at half maximum in the lateral direction is ≈0.7λ and in the axial direction ≈0.3λ.
Examining the average absolute value of the Poynting vector for a non-optimized beam over the scanning depth shows a decrease with penetration depth, as can be seen in Figure 9. Averaging was performed in the same way as for the energy density. Furthermore, the absolute value of the Poynting vector does not show an increase in value at shallow depths compared to the energy density in the previous section. This can be explained by the counter-propagating waves from the scattering medium, which lead to a resulting average decrease in the Poynting vector. We also find two regions within the medium with different exponential decays, at first with an exponential of ≈1/1.4λ until z = 8λ and with ≈1/4.5λ until the end of the medium. Both exponential behaviors have a slightly lower decay rate compared to the Beer-Lambert law. Phase optimization to enhance the absolute value of the Poynting vector has been done for three different numbers of optimization channels N ch ∈ {9, 137, 1009}. The four different optimization approaches described in Section 2.2 were used for optimization. The results are shown in Figure 9 in different colors and line styles. For scanning positions in front of the medium, all optimization techniques besides "g-opti" deteriorate the focus compared to the non-optimized beam and result in a lower value of | (S foc (z))| /| (S 0 )|. For a scanning depth greater than 4λ, all optimization techniques lead to an increase in the absolute value of the Poynting vector in the focus. At scanning positions inside the medium, we find that the differences between the global optimization technique "g-opti" and the other approaches increase with an increase in the number of optimization channels. Behind the medium, the "n-opti" and "g-opti" approaches lead to similar values of the average absolute value of the Poynting vector in the focus. The enhancement with the different techniques increases with the scanning depth for all approaches and is shown in Figure 9 in the bottom plot.

Discussion and Conclusions
In this study, we extended the two-step beam synthesis approach already used to investigate light propagation and optimization in two-dimensional media [21,37] to be able to model full-vectorial light propagation and optimization in three dimensions. The presented method is split into two main steps. First, a set of plane wave near-field solutions for a static medium is obtained by numerically solving Maxwell's equations, where the incident polarization and the propagation direction of the plane waves are determined through the Debye-Wolf theory. In the second step, the near fields of arbitrary incident beams scattered by the medium can be computed by adding the plane wave near-field solutions according to the angular spectrum of the incident beam. The separation of the calculation of the plane wave solutions and the synthesis of a complex beam makes it possible to reuse the plane wave solutions. Thus, with this method, the scattering of different incident beams can be efficiently calculated for the static medium under consideration. This includes scanning and optimizing the incident beam. Common methods would have to perform the complete numerical calculation again for each change in the incident beam.
We used this approach to simulate the scanning and phase optimization of a focused beam with a numerical aperture of NA = 0.45 inside a strongly scattering medium with a scattering mean free path of = 1.05λ and thickness 20λ. We examined the energy density and absolute value of the Poynting vector at the intended focus position over depth ranging from −2.5λ to 22λ, where 0λ was the location of the front surface of the scattering medium. The optimization of the incident beam was investigated for a different number of optimization channels N ch ∈ {9, 137, 1009} and different optimization schemes. To find an optimal phase pattern to regain a focus, we used the knowledge of the electromagnetic field at the intended focus position. The energy density and the absolute value of the Poynting vector in the focus showed two regions over depth, with first a steep and afterward a flatter, exponential-like decay. This can also be found in two-dimensional scattering media [21] but is not present in thin media or scattering media with a larger . For example, the absence of the two regions can be found in simulations for a scattering medium with = 1.95λ (not presented in this publication) and two-dimensional media [38]. Phase optimization by modulating the phase in k-space and applying the conjugate phase of the electric field in x or a random direction in the focus for each optimization channel can lead to a significant enhancement inside the scattering medium. An increased enhancement can be obtained by seeking a solution to the optimization problem to find the best direction in which to project the electric field to obtain the phase for phase conjugation. This requires two optimization variables. The best performance was always achieved when a solution for the "g-opti" was sought, but with the drawback that the N ch variables have to be optimized for the given objective function. An additional drawback of the presented simulation approach is that one near-field simulation takes around 30 min, and, for a larger volume and a larger numerical aperture, the number of required near-field solutions increases drastically. However, with increasing computing power and a larger number of parallel simulations, this approach can be accelerated.
The versatility of the presented method makes it possible to model arbitrary vectorial light beams. Possible applications of the approach are the investigation of different electromagnetic quantities inside turbid media, as shown in this study for the energy density and the Poynting vector. This can be of importance to quantify the performance of imaging modalities. Another possible use case is the calculation of optimal phase patterns computationally for real turbid media with a measured refractive index distribution to improve imaging modalities [23]. For large turbid media, the numerical solutions of Maxwell's equations are often too time-consuming to obtain. Therefore, it is also possible to use our proposed method but with plane wave solutions obtained with approximations such as the beam propagation method.  Data Availability Statement: Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.