Compensation of Beer-Lambert attenuation using non-diffracting Bessel beams

We report on a versatile method to compensate the linear attenuation in a medium, independently of its microscopic origin. The method exploits diffraction-limited Bessel beams and tailored on-axis intensity profiles which are generated using a phase-only spatial light modulator. This technique for compensating one of the most fundamental limiting processes in linear optics is shown to be efficient for a wide range of experimental conditions (modifying the refractive index and the attenuation coefficient). Finally, we explain how this method can be advantageously exploited in applications ranging from bio-imaging light sheet microscopy to quantum memories for future quantum communication networks.


Introduction
Imaging through diffusive or absorptive media is an extremely difficult task for an optician. It spans topics from imaging through atmospheric clouds, to in-vivo microscopy, to imaging in dense atomic gases. Several techniques are currently intensively investigated including speckle correlations [1] and transmission matrix reconstruction [2]. In bio-imaging, light-sheet or selected plane illumination microscopy allows selective illumination of tissues and fast 3D imaging of live organisms at the cellular scale which are much more precise than conventional confocal or multi-photon imaging [3,4]. However, the main limitation for the field-of-view of light-sheet microscopy is the penetration depth of the illumination through the tissue [5]. This fundamental limitation makes light-sheet imaging of deep tissues in living animals a challenging task, as it is complicated to illuminate deep structures effectively. In this paper we propose and implement a general method to resolve this strong limitation. We design a non-diffracting Bessel beam with the intensity in the central spot exponentially rising as function of the propagation, in order to exactly compensate for the losses in the tissue (or in any lossy medium) for illumination application.
Actually, the exponential attenuation of the beam due to scattering is common in optics. As soon as a beam propagates through a medium, scattering will inevitably degrade the signal. Our technique is of major interest not only for imaging biological samples but it could also be used to study disordered media or light propagation in dilute atomic clouds and non-linear crystals. In all these systems the exponential attenuation of a beam due to linear losses is a fundamental limitation [6][7][8][9].
Introduced by Durnin et al in 1987 [10], zero-order Bessel beams are one of the most representative non-diffracting solution of the Helmholtz equation together with Airy (or parabolic) beams [11]. The radial intensity profile of zero-order Bessel beams is described by the zero-order Bessel function of the first kind; a high intensity central peak is surrounded by an infinity of concentric rings of decreasing intensity. These optical fields can be considered as the superposition of an infinite number of plane waves with wave-vectors forming a cone around the z axis, the so-called Bessel cone. Although perfect Bessel beams are only mathematical objects as they should carry an infinite energy, spatially limited quasi-Bessel beams can be realized experimentally. Those beams have found various applications -in optical trapping [12,13], laser machining [14], nonlinear optics [15][16][17] and imaging [18,19] for example -as their central cores stay collimated on a distance that is orders of magnitude longer than the Rayleigh length.
The modification of the intensity profile in the propagation direction has been recently studied experimentally in the context of counterbalancing the intensity decay induced by light absorption in weakly absorbing dye solutions [20]. The general idea is to tailor the on-axis intensity of a Bessel beam. The use of an exicon (exponential intensity axicon) [21,22] and the generation of attenuation-resistant beams with computer generated holograms [20,23] have been reported. This last approach is known as frozen waves and results from the superposition of equal frequency Bessel beams produced by modulating only the amplitude of an incident plane-wave with a Spatial Light Modulator (SLM) [24][25][26][27][28].
In this paper, we report on a more general and versatile method based on both phase and amplitude shaping of an incident Gaussian beam.
It allows for compensating any absorption coefficients up to α = 200 m −1 , independently of the refractive index and of loss mechanism by using real space shaping with a reflective phase-only SLM.
We experimentally demonstrate the accuracy of this method in two very different media: a scattering sample based on a dilute solution of milk in water (index n = 1.33) and an absorptive sample of near resonance atomic vapor (n 1).
In the first section, we present the theoretical background of the method including a general approach to compensate for the refractive index of the medium. A detailed description of our experimental setup follows, as well as the measurement of the on-axis intensity profile of the central peak in air. In the next part, we present our results on the compensated absorption for two different media and show that our procedure is efficient independently of the refractive index and for a wide range of loss coefficients. Finally, we describe the potential improvement of 2 orders of magnitude on the field of view for light-sheet microscopy using this approach.

Shaping Bessel beams on-axis intensity
At a given position z on the optical axis (we assume z = 0 in the following), the electric field in the transverse plane of a laser beam E(x, y, z = 0) can be expressed as a function of its spatial spectrum S(k x , k y , z = 0) For a radially symmetric laser beam, this can be written as [29]: where J 0 is the zero order Bessel function of the first kind, r and k ⊥ stand respectively for the transverse radial coordinate and the associated spatial angular frequency. The spatial spectrum S(k ⊥ , z = 0) is also known as the Hankel transform of the electric field amplitude E(r, z = 0). As underscored in [29], this shows that a radially symmetrical field can be considered as a superposition of zero order Bessel fields. Each of these Bessel components then propagates without diffraction [10] as J 0 (rk ⊥ )exp(ik z z), where k z = k 2 0 − k 2 ⊥ is the longitudinal spatial frequency of the considered Bessel mode and k 0 = 2π/λ is the laser wave-vector (λ its wavelength). As a result, the on-axis electric field at point z, E(r = 0, z), can be obtained from (2) as: The angular spectrum of a perfect zero order Bessel beam is a ring in k space, with an angle θ 0 given by k z = k 0 cos(θ 0 ). The cone angle θ 0 sets the spot size i.e. the full width at halfmaximum (FWHM) of the central peak in the transverse intensity profile, which is equal to 2.27/(k 0 sin(θ 0 )) for a perfect zero order Bessel beam. Each of the Bessel modes coming in the spectral decomposition (2) propagates in free-space with a slightly different longitudinal wave-vector k z , as can be seen in (3). They thus merge with different cone angles at distinct positions along the optical axis. The on-axis electric field results then from the interference arising between the individual modes. If one wants to design a Bessel beam with a given on-axis intensity profile I(z) = |E(r = 0, z)| 2 along the optical axis, the spatial spectrum S must be engineered according to the following formula : The spectrum S is centered around the axial wave-vector of the target Bessel beam k z0 = k 0 cos(θ 0 ). This formula gives a physical insight about the engineering process that will be used to overcome attenuation. The initial electric field E(r, z = 0) that will produce a Bessel beam with a given cone angle θ 0 and an on-axis intensity profile I(z) can be evaluated using (2) and (4). In the following, we briefly describe how to generate the target beam by real-space shaping of an incident Gaussian beam on a Spatial Light Modulator (SLM). Technical details can be found in the technical details section at the end of this manuscript. Fourier space shaping may also be considered [29]. However, as the intensity distribution of a Bessel beam in Fourier space is a thin ring, the small overlap between the latter and the incident Gaussian profile would filter out most of the incident energy. Much higher efficiency can then be obtained using real space shaping presented here.
We define z = 0 to be the SLM plane position along the optical axis. Discretizing the electric field according to the SLM matrix (N x × N y ), the target electric field E(i, j, z = 0 + ) right after the SLM can be decomposed in amplitude A(i, j) and phase Φ(i, j) (where 0 ≤ i ≤ N x and 0 ≤ j ≤ N y stand for the pixel coordinates). As suggested by Davis et al. [30], locally reducing the phase wrapping contrast allows for a modulation of the light scattered in the first diffraction order, using a single hologram. We apply this technique with a phase-only SLM. The expression of the SLM phase mask Ψ is given by [31,32]: The function F contains the phase information of the target electric field and Φ g stands for the grating linear phase ramp, used to separate the different diffraction orders in Fourier space. The modulo operation provides the phase wrapping. The diffraction efficiency is locally tuned by the modulation function M (0 ≤ M(i, j) ≤ 1). The complex amplitude of the field diffracted in the first order can be expressed as follow [31,32]: where A in is the amplitude of the incident laser beam on the SLM. By identifying E 1 with the target electric field, one can obtain the functions F and M (see the technical details section A at the end of this manuscript).
In principle, arbitrary on-axis intensity profiles can be generated using the method described above. In the following section, we introduce the target profile I(z) we use to maintain the central peak intensity constant along the propagation in a uniform and linear lossy medium. Our approach is independent of the loss origin; we demonstrate its validity for both absorbing and scattering type of losses. Let L and α stand respectively for the propagation length and the linear attenuation coefficient of the medium. According to the Beer-Lambert's law, the transmittance t of the medium decays exponentially with the propagation distance: T = exp(−αz). Therefore, to compensate for these losses, the on-axis intensity should exponentially increase along the propagation such as I(z) ∼ exp(αz). We ramp the on-axis intensity up (from 0 to I(z 1 ) = I 0 ), until the entrance plane position z 1 , before exponentially increasing it over the length L. We then make it go back to 0 smoothly. The full on-axis target profile we designed is described as: For all the measurements we performed, we set z 1 × G 2 = 1.5 cm, z 2 = z 1 + L G 2 (with L = 7.5 cm) and z 4 = 3 z 1 + L G 2 , where G = 0.5 stands for the telescope demagnification factor which optically conjugates the SLM and the z = 0 planes. C 1,2 and z 3 are constants chosen in order to make the profile continuous and differentiable; they are defined in the technical details section at the end of this manuscript. In the following, the profile has been normalized to 1 dividing I(z) by the maximum intensity I max = I 0 (1 + exp(αL)). The spatial spectrum associated to this on-axis profile is analytically derived in the appendix. We obtain the target electric field by computing the inverse Hankel transform using (2).
When the linear refractive index n of the medium is not equal to 1 (as implicitly assumed before), the target Bessel beam will undergo refraction at both the entrance and the output plane of the medium. Applying the Snell's law at the entrance (resp. the output) plane, we get: sin(θ i ) = n sin(θ r ), where θ i and θ r stand respectively for the incident and refractive cone angle of a given Bessel mode. Using the transverse spatial angular frequency k ⊥ = nk 0 sin(θ), we find that k (i) ⊥ = k (r) ⊥ . Therefore, according to (2), the transverse shape of the target Bessel beam is not modified by successive refractions [33]. Nevertheless, the cone angle is modified when the Bessel beam enters the medium.
As n>1, the inner cone angle θ r is smaller than the external one and the Bessel beam will cover a distance longer than in air. This stretching of the beam inside the medium will necessarily reduce the compensation coefficient by a factor n. To counteract this effect, we constrict the exponentially rising part of the target on-axis profile beforehand by a factor n (as suggested in [34]). In other words, we replace L by L/n and α by αn in the second line of Eq. (7). By doing so, the stretching of the beam will compensate exactly for the exponential attenuation in the medium, as shown in Fig. 1. This compensation procedure generalizes easily to more complex situations when several layers of materials (with different attenuation coefficients and refractive indices) are involved. When the medium refractive index n is larger than 1, the calculated profile must be stretched by n to compensate for the refraction. Blue, red and grey circles are obtained by solving numerically the evolution of the Bessel beam. Blue dots are calculated for n = 1 in air and red dots correponds to the strectched case for n = 1.33 (without losses). To verify that we are able to compensate losses in a medium with n>1, we plot in gray the propagation in a lossy medium of refractive index 1.33. For red and greys data, this can be deduced from the vacuum case by stretching the z axis by a factor n between z 1 and z 2 . When losses are present, the on-axis intensity remains constant along the propagation (black dots) as expected.

Experimental setup
Our experimental setup is shown in Fig. 2. A continuous-wave laser beam, produced by a tapered amplifier laser system, gets 4 times magnified by a 4-f telescope system (lenses L 1 and L 2 ) and is spatially filtered in the Fourier space of L 1 . The resulting 6.6 mm diameter radially symmetric Gaussian beam reaches the center of the SLM chip with a normal incidence. The SLM used for the experiment is a liquid crystal on silicon (LCOS) phase-only modulator, with an effective area of 1272 × 1024 pixels and a pitch of 12.5 µm. After shaping, a 50 : 50 non-polarizing beam splitter separates the diffracted beam from the incoming one. Another 4-f telescope (L 3 and L 4 , NA = 0.017) optically conjugates the SLM and the z = 0 planes with a demagnification factor G = 0.5. The choice of the telescope lenses L 3 and L 4 , as well as the demagnification factor G, is conditioned by the length of the lossy medium we are dealing with. For biological applications, G should be divided by 10 at least, as pointed out in section 2. The first order diffracted beam is then selected by masking the zero and the higher order ones in the Fourier plane of L 3 . The Bessel beam finally starts forming from the focal plane of L 4 (at z = 0) and propagates through a lossy medium. The output plane of the medium is imaged by a third 4-f arrangement (lenses L 5 and L 6 , NA = 0.042) onto a microscope objective which is set up on a computer controlled translation stage. By moving the objective along the optical axis, we can monitor the Bessel beam evolution along z. The last lens (L 7 ) images the plane we look at on the CMOS camera. The magnification factor G of the whole imaging system is 13.6 ± 0.1.

Results and discussion
The first step to verify the compensation of the beam attenuation is to compare the transverse and longitudinal intensity profiles of the experimentally measured beam with the target ones from simulations, in air. We design the Bessel beam to overcome 96% attenuation over a lossy, 7.5 cm long medium. The 2D map in Fig. 3 and profiles in Fig. 4 are obtained by scanning slowly (v = 2 mm.s −1 ) the microscope objective along the z axis. Both the transverse, Figs. 4(a) and (b), and the longitudinal, Fig. 4(c), measured intensity distributions of the tailored beam (blue circle and line) are in excellent agreement with the simulation (dashed black lines). The target profiles (dashed black lines) are obtained by numerically solving the evolution of the transverse electric field from z = 0 to L with the second order split-step method. We take as initial condition a field with the SLM imprinted phase Ψ and the radially symmetric Gaussian envelope of the SLM input beam.
To accurately determine the central peak intensity along z as presented in Fig. 4(b), we fit with a Gaussian profile the region delimited by the two white dashed lines on both sides of the central peak (Fig. 3) as illustrated in Fig. 4(a) (red line). The width of the peak along the propagation is found to be constant (±5%), as shown in Fig. 4(b). This measurement demonstrates that we are therefore able to control the longitudinal intensity profile without altering the non-diffracting behavior of the Bessel beam. Along the propagation axis, we have also measured a transversal cut of the intensity (at r = 0) and plotted it as function of z. We observed that our experimental data are increasing exponentially as expected and follows very well the calculated profile. Nevertheless, small intensity oscillations can be observed at the beginning of the measured on-axis profile  Fig. 4(c). They are due to high longitudinal frequency truncation, as k z is upper bounded by the laser wave-vector k 0 . The oscillation amplitude can be reduced further by increasing the Bessel cone angle θ 0 . In our setup, we are limited by the mirror size as the beam will start to clip and will loose its radial symmetry.
The lossy medium is then positioned on the beam path. By fitting the on-axis intensity profile with the function of (7), we find the position z 2 where the medium output plane should set. We then move the lossy medium cell along the optical axis until imaging this plane on the camera. The 1 mm depth-of-field of the imaging system and the standard deviation on the fit parameters translate into an uncertainty of ±2 mm on the medium output plane position.
Three different media (contained in three different glass cells) have been used to check our ability to compensate for the attenuation of the Bessel peak intensity along propagation. Two cells are filled with isotopically pure Rubidium vapor (the first (7.5 cm long) with 87 Rb only and the second (2.5 cm long) with 85 Rb only); the third one (2.5 cm long) contains a diffusive water-milk mixture. Rubidium cells are heated up to 140 • C. At this temperature, the atomic density is large (n a 2 − 5 × 10 13 atoms/cm 3 ). By tuning the laser frequency ν 0 over the D 2 Rubidium absorption lines, we can change the transmission over several orders of magnitude, without affecting significantly the refractive index n Rb . The latter is estimated theoretically taking the Rubidium hyperfine structure and the Doppler broadening into account [35]: n Rb (ν 0 ) 1.00±0.02 (scanning ν 0 over the whole absorption spectrum). The transmission of the water-milk mixture can be tuned changing the milk concentration. Remaining under highly diluted condition, the medium refractive index stays close to the water one n w 1.33. As explained above, we should balance in this case the change of refracting index stretching the Bessel beam along the optical axis, replacing beforehand in the target profile L and α with L/n and α n respectively.
We design the target profile to overcome attenuation over 7.5 cm long materials, whatever the length of the cell we use. The overall Bessel power is reduced to keep the input peak intensity lower than the Rubidium saturation intensity (I sat 2.5 mW/cm 2 for linearly polarized laser beam). We finally measure both the peak intensity in the entrance plane (without cell) and in the output plane (with cell) to evaluate the transmission trough the medium. We perform the fitting procedure on five different images of the central spot with a 2D-Gaussian function as detailed before. The measured transmission is shown in Fig. 5. The experimental data obtained with the 87 Rb vapor cell, the 85 Rb one and the water-milk mixture are respectively plotted in blue stars, orange circles and grey diamonds. The 4% reflectively of the cell windows has been taken into account. The black dashed line represents a perfect on-axis compensation (T = 1); most of the experimental points lie right under it. The small discrepancy comes from the input plane intensity measurements rather than from the output plane ones. The oscillations of the intensity at the medium input plane, visible on Fig. 4(c), adding to a positioning uncertainty of 4 mm of the input plane, lead to the errorbar uncertainty reported in Fig. 5. The red lines represent the transmission of a non-saturating collimated Gaussian beam with respect to the attenuation coefficient α for a 2.5 cm and a 7.5 cm long lossy medium.

Applications
Compared to previous experiments, in which compensation for 30 % and 10 % attenuation have been achieved using exponential intensity axicon [34] and attenuation-resistant frozen waves [36] respectively, we manage to maintain the Bessel beam on-axis intensity quasi-constant along its propagation by compensating for an attenuation of 200 (equivalent to a transmission of 5 × 10 −3 ). This is a crucial advantage for biological applications based on scattered light as the diffusive coefficients of biological tissues observed in light-sheet microscopy range typically' from 50 to 200 cm −1 [37,38]. For comparison, let's assume that the sample transmission is 6 × 10 −3 under Gaussian illumination (as for the last orange circle of Fig. 5) for a diffusive coefficient of 200 cm −1 . The length of the sample is then equal to − ln (6 × 10 −3 )/200 250 µm. By increasing the demagnification of the 4-f telescope {L 3 , L 4 } by 10, we could easily decrease the length of the Bessel zone by a factor 100 (as the length varies with G 2 ) and, using the same target profile, uniformly illuminate biological tissues over hundred of microns up to diffusive coefficient α of 200 cm −1 . For such attenuation coefficient, the field of view would then be more than 100 µm longer than the best one obtained in the literature so far [38] (or even more if partial attenuation-compensation is considered). Moreover, we extend the approach of [38] with Airy beams to Bessel beams.
For biological applications maintaining the transverse size through diffusive medium is important. To verify that our technique allows this, we measured the transverse profile after the diffusive medium (water+milk). In this situation we did not observe any significant broadening larger than in air (Fig. 4). This can be understood simply because we continuously add energy in the spot center to compensate absorption (from the outter rings) and not the tail of the spot. The center spot is therefore continuously reshaped towards the initial size.
On the other hand, any light-matter interface that relies on maximising interactions over a long distance along the propagation axis, such as EIT-based quantum memory or gradient echo memory [39], will benefit from this technique to improve significantly the effective optical depth. As it is known that the quantum memory efficiency (for Raman schemes) is proportional to the Rabi frequency of the coupling field [40], compensating for losses will immediately increase the storage efficiency especially in the case of an ultra-high optical depth medium [41]. This of course implies that the region used for the memory is limited to a cylinder with the diameter of the Bessel beam and assumes that atoms are slow enough to stay within the interaction region during the memory time. Finally, attenuation compensated Bessel beams are needed for the generation of stationary non-diffracting potentials in fluid of light experiments in the propagation configuration [42], where superfluidity has recently been observed [43,44]. So far, exponential attenuation of the defect was the main limitation of these experiments in hot atomic media. With the approach described in this paper, future experiments can be envisioned where a Bessel beam pumps a vapor close to resonance and modifies the medium refractive index locally in the transverse plane and uniformly along the beam axis.

Conclusions
We have reported the shaping of the longitudinal intensity profile of Bessel beams using a phase-only SLM. We have shown that this method can be used to compensate the Beer-Lambert law and generate a constant intensity profile along the propagation direction in a lossy medium. We verified that our approach is robust independently of the loss mechanism at play in the medium and for a wide range of conditions (various refraction indices and attenuation coefficients). The results are in agreement with numerical simulations and the most crucial limitations are clearly identified. We identified two applications where this can be advantageously used: in bio-imaging and in quantum optics. Finally, this method can be easily generalized to tailor any kind of on-axis intensity profile.

SLM hologram design
One can obtain the function M and F of (6) by solving the following system: The inverse sinc function (sinc −1 ) is defined on [−π, 0]. Computing it for each points of the hologram is usually demanding (N x × N y operations). However, if both the incident and the first order diffracted beams are radially symmetric, we only need to determine the radial profile of the modulation function. For beams centered in the SLM matrix, (8) can be simplified such as: where i is an integer running from 0 to N x /2 (for N x ≥ N y ). Using a circular interpolation, M can be entirely reconstructed from m, computing the inverse sinc function for N x /2 points only instead of N x × N y . In practice, we start with the clean-up of the incident laser beam, filtering out in Fourier space its high k-vector components with a small pinhole aperture. Afterwards, the widths of the input Gaussian beam is radially symmetric in the SLM plane (ω x,y 3.3 ± 0.1 mm). For non-evanescent modes, k z should lie in the interval [k min , k 0 ], where k min is defined by the numerical aperture (NA) of the imaging system as k min = k 0 1 − NA 2 . As explained in Ref. [32], the target intensity profile should not vary on a length scale smaller than ∆z = 4π/(k 0 − k min ) to avoid significant frequency truncations in the associated spectrum, leading to undesirable oscillations in the measured on-axis intensity profile.

On-axis intensity profile and spatial spectrum
The target profile we designed to compensate on the optical axis the attenuation of the central peak intensity is given by (7). The SLM mask is optically conjugated with the z = 0 plane by the 4-f telescope formed with the lenses L 3 and L 4 . For clarity, we have ignored the telescope magnification factor G = 0.5 in the analytical derivations above, but we need to consider it in practice when we compute the target intensity profile. For all the measurements we performed, we set z 1 × G 2 = 1.5 cm, z 2 = z 1 + L G 2 (with L = 7.5 cm) and z 4 = 3 z 1 + L G 2 . We choose C 1 and C 2 in order to make the target profile continuous and differentiable at z 1 and z 2 : the constant C 1 is obtained by solving the equation tan(x) = 2x αL (deriving from the diffentiability condition at z 1 ) and C 2 = sin −1 I 0 /I max exp αL 2 . The maximum intensity I max is obtained for z 3 = z 2 + 2C 2 α tan(C 2 ) . The spatial spectrum associated with the target profile can be derived analytically using (4). As all the parts composing the target profile can either be expressed by an exponential rising function or a sine square function, computing the spatial spectrum associated to the following generic functions is sufficient: The derivation of the associated spectra S where l = z j − z i and δk = k z0 − k z . We finally obtain the spectrum summing the spectral contributions coming from the different parts of the profile: S = S (0,1) sin + S (1,2) exp + S (2,3) sin + S (3,4) sin .

Clearing of the refractive stretching
Let's assume that the target Bessel beam enters at z 1 a material with an attenuation coefficient α and a refractive index n. In order to counteract the refractive stretching of the beam, let's also replace L with L/n and α with α n in the second line of (14). Using (4) and the change of variable z →z = n(z − z 1 ), we can derive the spectrum S 1,2 associated to the exponential rising part of the on-axis profile (between z 1 and z 2 ): The on-axis electric field E(r = 0, z) is related to the spatial spectrum S by the Fourier transform (3). In practice, k min and k 0 set respectively the lower and upper bounds of the integral coming in this equation (in air). Using (14) and (3) and the change of variablek z = (k z − k z0 )/n, we can derive the on-axis electric field E 1,2 (r = 0, z) associated to S 1,2 : As z lies in the interval [z 1 , z 2 ] and z 2 = z 1 + L/n, δz = n(z − z 1 ) varies from 0 to L. The phase Φ l = k z0 (z 1 + δz/n) is the phase accumulated by the Bessel beam along its propagation until z. The medium is supposed to be linear; this phase term is therefore the only expected. The inside brackets (15) should then be real. Let's divide the integral in two parts I 1 and I 2 as follow : From (16) and (17), we derive the real parts of I 1 and I 2 : Re (I 1 ) = 0 and Re (I 2 ) = exp αδz 2 . The on-axis electric field E 1,2 (r = 0, δz) is finally given by : By replacing L with L/n and α with α n in the expression of the target on-axis intensity profile (7), we manage to overcome the refractive stretching of the Bessel beam and compensate for the good attenuation coefficient α.