We present a c omputationally efficient full-wave spectral model of OCT -scan formation with the following capabilities/features : (i) the illuminating beam may have arbitrary phase -amplitude profil e with allowance of a sharp diaphragm ; (ii) paraxial approximation that limits the degree of focusing/divergence is not used ; (iii) the broadly used approximation of ballistic (single ) scattering by discrete scatterers is assumed without additional limitat ions on the density of scatter ers, their distribution in space and scattering strengths with possible frequency -dependence ; (iv) besides rigorous accounting for the influence of focusing/divergence of the waves, factors describing the wave decay can r eadil y be introduced by analogy with Monte -Carlo approaches; (v) a rbitrary measurement noises can easily be added to simulate required signal -to-noise ratios. The model also allows one to account for arbitrary (e.g., random or flow/deformation - induced) displace ments of scatterers between subsequently generated scans. Thus, in view of the above -listed features the model can be characterized as comprehensive in the fram ework of ballistic scattering by discrete scatterers . The model is computational ly efficien t due to the use of only rapid summation s and fast Fourier transforms. Main model features are illustrated , including simulati ons of OCT scans for a nearly non -diverging Bessel beam and a focused Gaussian beam , with the possibility to introduce at the tissue bo undary arbitrary aberrations represented via Zernike polynomials often utilized for describing aberrations in ophthalmology. The unprecedented flexibility and high computational efficacy of the model o pen a broad range of possibilities for studying OCT -scan properties and developing new methods of their processing for biomedical diagnostics . Keywords: Optical coherence tomography ; OCT ; OCT -image formation ; Numerical simulation of OCT ; Multimodal OCT 1. Introduction OCT signal simulation suggests unique possibilities for develop ment/optimization of various OCT mo dalities in highly controllable , easily modifi able conditions , which cannot be attained in physical experiments with real biological tissues and even phantom s. Numerous research groups make significant effort in develop ing various tools to enable modeling of OCT signals for general and specific purpos es, e.g.1-12. Main approaches to OCT -scan modeling roughly can be divi ded in such groups as Monte -Carlo methods2,4,5, utilization of the point -spread function concept11,12 and full- wave methods describing in various approximations the forward propagation of the optical wave towards scatterers and then its backward propagation and collection over the receiving aperture6,7,8. In many models the scatterers are approximated by discrete scattering centers , properties of which can be more complex than point scatterers (e.g., accounting for frequency dependence) to better imitate real tissues . The use of discrete scatterer s concerns all three groups of simulation methods (e.g.9,16,17,18). This is quite a meaningful approximation , even if real scatterers are not literally point -like o bjects; nevertheless, th ey still can be treated as discrete sub-resolution particles in OCT scans and be reasonably corresponded to c ells in biological tissues. Also it should be recalled that the very fact of successful realization of OCT imaging is based on utilization of balli stic back -scattering . Even within the single -scattering approximation one can study many important specific features of OCT images (e.g., evolution of speckle patterns in the presence of scatterer motions, statistical properties of OCT speckles, etc.). Mul tiple 2 scattering as a degrading factor in OCT often is of secondary interest, because numerous OCT -scan features can be well captured even in the single -scattering approximation. It is noteworthy that, d epending on the situation , speckles in OCT may be either considered as a degrad ing factor or play a useful role of high diagnostic value1. For instance, motion - induced speckle variability is the basis of OCT -angiography19. Therefore, it is reasonable that the majority of scan - formation models in OCT are based on ballistic approximation. For the illuminating -beam description in such models, usually various additional simplifications are made , which may be mor e or less realistic. In simpler variants like1,16, the illuminating beams may have uniform phase in the cross section (which still may be a reasonable simplification fo r OCT -systems with weakly focused beams without pronounced curvature of wave fronts, whereas the amplitude profile may be variable , e.g. Gaussian ). Even using such simplifi ed beam geometries as in13,16, it is already possible to simulate important features of real OCT -scans formed by beams with weak focusing/divergence and obtain useful results for developing new OCT modalities such as elastography20,21,22,23. For weakly focused beams, t he validity of such sim plified models is also confirmed by similarity with physical experiments and results of more advanced models like17,18 that belong to full-wave approaches and are based on rigorous consideration of Gaussian beams with allowance for pronounced focusing. It should be pointed out that f or simulations of realistic densities of scatterers (imitating the density of cells in biological tissues) in biologically relevant sizes o f visualized regions (~hundreds of microns and greater), sufficiently high computational efficacy of the simulation s is of key importance , especially when large series of OCT -data volumes are required to study effects of moving scatterers. In this context utilization of Gaussian beams in modeling is rather attractive because of the possibility to describe their forward propagation analytically. Collection of the back - propagated waves over the receiving aperture can be made either using numerical integration (like in17), or even analytically for non -diaphragmed beams18. Certainly, the possibility to pass from numerical integration to analytical description of the back -propag ated wave and its collection over the receiving aperture like in18 dramatically improved the computational efficiency. However, this improvement was enabled at the expense of a limitation : the Gaussian beam diameter should be significantly smaller than the illuminating/receiving aperture . For non -Gaussian beam forms (and even for a Gaussian beam with an aperture limited by a diaphragm or somewhat distorted by aberrations ), feasibility of analytical descri ption of the illuminating -beam propagation towards scatterers becomes either impossible or requires a radical reformulation of all basic equations . Numerical realizations (concerning both the forward/ back ward propagation and collection of t he back scattered fields over the receiving aperture ) is very computationally demanding . Consequently, development of realistic and flexible models of OCT -scan formation requiring minimal simplifications/approximations and at the same time enabling high com putationally efficiency remains challenging and highly demanded . In what follows we demonstrate that using the same basic approximations (i.e., the use of discrete scatterers and ballistic scattering) as in17,18 and many other models, e.g.9,14, it is possible to reformulate the OCT -scan formation model , such that the abov e-mentioned limitations can be completely lifted . As a result, the modified model acquires the following previously unavailable features : (i) it may simulate OCT scans using arbitrary illuminating beams (Gaussian, Bessel , etc. ), including the possibility t o introduce sharp diaphragms ; (ii) arbitrary aberrations of the illuminating beam can be introduced (such that forward and backward propagation through the aberrating screen is rigorously accounted for); (iii ) the widely accepted paraxial approximation is not required , so that wide -angle beams can be used; (iv) arbitrary motion s of scattere rs between subsequent scans can be readily accounted in the model; (v ) although the widely used approximation of ballistic scattering at discre te scatterers is used, there are no other limitations on the density of scatterers, their distribution in space and scattering strengths including possible frequency -dependence. The proposed modified model described in the next section uses full-spectral representation of the illuminating /scatterered fields (both in axial directions , like is physically used in spectral -domain OCT, and the description of the transversal profile of the wave is also made in the spectral form instead of the spatial -domain representation ). Due to su ch a modification , the new formulation of the model becomes very straigh tforward and strikingly compact. The possibility of using arbitrary beam profiles is not made at the expense of computational 3 efficiency ; on the contrary, the model computationally is very fast. To the best of ou r knowledge, there were no previously discussed models that could simultaneously suggest such a range possibilities . Then we demonstrate the possibilities of the proposed model, initially by validating it via comparison with th e reference model18 that is based on exact analytical solutions for non-diaphragmed focused Gaussian beams. We demonstrate that in approximation of sufficiently lar ge aperture (required for validity of18) for the same set of scatterers , the proposed new formulation yields the results coinciding with a high accuracy with the exact analytical results18. The universality of the ne w model is then illustrated by obtaining OCT images for a Gaussian beam and a Bessel one (again for the same set of scatterers) . Furthermore, in view of high importance of accounting for aberration s for some applications of OCT (e .g., in ophthalm ology24), the unique flexibility of the new model is demonstrated by generating images with additionally introduced aberrations represented by Zernike functions often used for descriptions of eye -lens aberrations in ophtha lmology . Finally, the model usage in configuration with moving scatterers (like in angiographic /elastographic applications) is demonstrated. 2. Methods of f ull-wave spectral description of OCT -scan formation In the proposed description, each spectra l component of the i lluminating -beam with wavenumber k is characterized by a complex -amplitude distribution ),,( kyxL at the output lens (with allowance for a sharp diaphragm) . The lens is located in the (x,y) plane has the axial coordinate dzz and z -axis is directed to the tissue depth along the axis of the illuminating beam as shown in Fig. 1. Coordinate dz may either directly correspond to the tissue boundary or the boundary of the immersion layer (in which scatterers are absent). The lateral coordinates ),(0 0yx of the beam axis define the position of a single A -scan, so that lateral scanning required for obtaining 3D volume of OCT data corresponds to variation in ),(0 0yx . Each discrete scatterer is characterized by coordinates ),,(sss zyx . Like in real spectral -domain OCT we represent the illuminating beam as a set o f spectral components with wavenumber s k. The forward propagation of the axial spectral components an d their back -propagation from scatterers is often described in the spatial domain17,18. However, in principle the sp ectral approach can be used to describe the optical -wave evolution in axial direction as well . For example, after obtaining 3D images in17,18, the numerical refocusing was made by recalculating the wave fronts considering their spatial spectra in the lateral direction. In what follows we apply such a spectral -domain description (both in the axial and lateral directions) to describe the wave fields for all stages: first, for illuminating -beam propagation towards the scatterers and then for backward propagation of scattered fields and their collection over the ill uminating - receiving aperture. Fig. 1. Schematically shown geometry of the illuminating beam (not necessarily focused like in the figure) with the lateral coordinates ),(0 0yx of the axis directed to the tissue depth along the z -axis of the coordinate system. Coordinates ),,(sss zyx characterize positions of scatterers. 4 In such a rep resentation, the complex -amplitude distribution ),,( kyxL over the input/output aperture is characterized by the spatial spectrum ),,( kvug (for simplicity , below we omit constant factors before the integrals ):   Dyxvyuxidxdy ekyxL kvug ,) (),,( ),,( . (1) Here, for brevity we denote the wave vector projections ukx and v ky . After propagation of the illuminating beam between the input plane dzz and plane szz corresponding to the depth of a scatterer with coordinate sz , the illuminating -wave spectrum can be written as ),,(),,( kvugkvuhs , where the propagator function ),,( kvuhs has the form 2 2 2) (),,(vuk zzi sd sekvuh . (2) Note that a similar method for describing the axial evolution of the transversal spectrum of the illuminating beam using the propagator (2) was recently applied in paper15 to describe degradation of the point -spread function in the presence of refractive -index fluctuations . Since for the scatterer with coordinates ),,(sss zyx , the transversal spatial spectrum of the incident wave is ),,(),,( kvugkvuhs , its Fourier transform over spectral coordina tes ),(vu corresponds to the spatial - domain representation of the field incident on the scatterer : dudv ekvuhkvug kyxU kvuyyvxxui s ss s  2 2 20 0 )] () ([ 0 0 ),,(),,( ),,( , (3) The possibility to describe the field ),,(0 0 kyxUs in the spectral form (3) is a very useful feature, because the transversal spectrum ),,( kvug can be readily found numerically for an arbitr ary beam profile ),,( kyxL only once for all scatterers and does not require additional assumptions (e.g., about a Gaussian profile, absence of a diaphragm, etc. ) which are usually needed to simplify analytical consideration . In can also be noted that b y ana logy with the spectral function ),,( kvug defined by Eq.(1) (i.e., the Fourier transform of the initial beam prof ile ),,( kyxL ), for the propagator ),,( kvuh in the spectral form , its spatial -domain counterpart ),,( kyxHs can be introduced via the inverse Fourier transform over spectral coordinates ),(vu :  dudv ekvuh kyxHvyuxi s s) (),,( ),,( (4) Bearing in mind that according to Eq.(3), the transversal spectrum of ),,(0 0 kyxUs is given by the product of the Fourier transforms ),,( kvug and ),,( kvuhs corresponding to the spatial -domain functions ),,( kyxL and ),,( kyxHs , according to the theorem about the product of spectra, function ),,(0 0 kyxUs can also be expressed as the following convolution of ),,( kyxL and ),,( kyxHs :    Dyxs ss s dxdykyyxxHkyyxxL kyxU ,0 0 0 0 ), , ( ), , ( ),,( (5) 5 In Eq.(5) the integral is found over the aperture D of the lens (see Fig. 1). In fact ),(yxHs is the Green’s function for a localized scatterer located at ),,(sss zyx and Eq. (5) is a spatial -domain counterpart of the spectral representation of the incident wave given by Eq. (3). The scattered wave for a discrete scatterer located at ),,(sss zyx is proportional to the incident -wave amplitude ),,(0 0 kyxUs . This scattered field propagates back to the illuminating/receiving aperture and is collected by the lens , which is again described the convolution of the Green’s function and receiving aperture similar ly to Eq.(5) . Then the amplitude ),,(0 0 kyxBs of the spectral component with wave number k collected at the receiving aperture is     Dyxs ss s s dxdykyyxxHkyyxxL kyxUkyxB ,0 0 0 0 0 0 ), , ( ), , ( ),,( ),,( (6) Here, the pre-integral factor corresponds to the incident -field amplitude, whereas the integral factor describing the collection of the back -scattered wave is identical to integral in Eq.(5) that describes the forward propagation. In view of this, Eq.(6) cab be rewritten in the compact form as the incident field squared : ),,( ),,(0 02 0 0 kyxUkyxBs s  . (7) At this point it should be recalled that Eq.(6) describes the received amplitude ),,(0 0 kyxBs for a single spectral component with the wavenumber nkk , whereas in spectral -domain OCT , the optical -beam contains a large number of illumination spectral components zikn exp with wavenumbers nk . In previous papers16,17,18 (where spectral amplitudes ),,(0 0 kyxBs were found differently than described above ) we discussed in detail that superposition of N n ..1 spectral components with wavenumbers nk forms a pixelated A -scan com posed of N q ..1 pixels. This superposition mathematically correspond s to Fourier transformation over the axial wavenumbers , which is performed in spectral - domain OCT scanners to form in -depth A -scans. Using the notations similar to those used in works16,17,18 we can write the following expression for the complex -valued amplitude ),,(0 0 qzyxA of q-th pixel with the axial (depth) coordinate qz and coordinates ),(0 0yx in the lateral plane:      n sq n n s q zHni kSkyxB zyxA2exp)(),,( ),,(0 0 0 0 (8) In Eq.(8) summation is made for the contributions of all spectral harmonics N n ..1 also taking into account their amplitudes )(nkS in the illuminating -beam spectrum. The second summation in Eq.(8) is made over contribution s of all scatterers located within the beam. In a more compact form Eq.(8) can be written as the Fourier transform over nk      sn n k q kyxBkS FFT zyxA s n,, )( ,,0 0 0 0 , (9) where the spectral form )(nkS of the illuminating beam is explicitly retained . It should be emphasized that although spectral Eq.(3) and spatial -domain Eq. (5) equivalently represent the same amplitude ),,(0 0 kyxUs , the possibility of its efficient evaluation in the spectral form (3) is the central point f or computational efficiency of the described model. Indeed, spectrum ),,( kvug readily allows one to represent an arbitrary phase -amplitude form ),,( kyxL of the illuminating beam via Eq.(1). Therefore, the spectral representation confers the 6 above -described model unprecedented computational speed for treating arbitrary initial beam profiles ),,( kyxL without radical reformulation. In particular, not only Gaussian beams, but other beam forms can be used (e.g., Bessel or Laguerre -Gauss beams that attract high interest in recent years). Also aberrations of the illumination beam at the boundary of the visualized region can easily be introduced, which is of interest for OCT applications in ophthalmology, etc. For any beam profile ),,( kyxL , the transversal spectrum ),,( kvug should be found only once and used for all A - scans and scatteres . Next, it should be pointed out that the spectral representation o f the propagator ),,( kvuhs can be used in the exact form (2) and does not require the use paraxial approximation , which i s often necessary for obtaining analytical result s in other models . Finally, one can readily imitate motions of scatterers by varying coordinates of scatterers ),,(s s szyx for subsequently obtained scans. Notice that for motionless scatte rers, esti mation of Fourier transform given by Eq. (3) can be made using rapid FFT-algorithms to obtain (via Eq s. (7) and (8)) spectral amplitudes ),,(0 0 kyxBs for the entire set of A -scan positions ),(0 0yx within the simulated 3 D volume of OCT -data. For displaced scatterers, the data should be recalculated and the total computation time will depend on the total number of ),(0 0yx -coordinates of the illuminating beam within the 3D data set and the required number of recalculations for moving scatterer s. Depending on the ir velocities , for certain cases, displacements may be sufficient to consider either between entire 3D volumes , in other cases variations in ),,(s s szyx may be significant even on inter -A-scan intervals requiring larger amount of recalculations. In the next section we present examples of utilization of the proposed model to illustrate its main new features /capabilities . 3. Simulation e xamples to demonstrate main model capabilities In the b elow -presented examples we assume that the central wavelength of the OCT system is 1.3 m and spectral width is 90 nm, which are “typical” parameters for OCT . Similarly to18 we simulate the OCT data for rectangul ar volumes with sizes 128x128x128 pixels corresponding to biologically -relevant physical sizes of 512 m in depth and 192 m in the lateral directions . Other parameters of the illuminating beam and scatterers are specified for each of the examples. First, for validation of the proposed model, it makes sense to compare the result of the new spectral formulation given by Eqs. (1) , (3), (7) and (9) with the reference model18, because the latter is based on exact analytical solutions for the spectral amplitudes ),,(0 0 kyxBs obtained for a non-diaphragmed Gaussian beam ( i.e., assuming that the lens aperture is sufficiently big and can be considered infinite). For comparison of the present model with18, we accepted identical parameters of a foc used Gaussian beam and generated scans for a set of 7 scatterers with i dentical lateral coordinates and different depths, including the focal depth 128z µm. Figure s 2a and 2b demonstrate that the described full spectral - domain formulation based on discrete Fourier transforms give s the spatial form s of point -spread functions very close to the results obtained using the exact analytical model18. The depth profile s shown in Fig. 2c confirm that the two methods demonstrate excellent coincidence of the signal amplitude , so that minor differences can be observed only in low - amplitude zones (at a level about -60 dB below the maxima and lower) sufficiently far from the scatterers . This coincidence with the exact analytical results for a pronouncedly focused Gaussian beam confirms the reliability of the proposed model , for which, however, the main new feature is the ability to describe OCT signal s formed by arbitrary beam s, not necessarily axially symmetric ones . 7 Fig.2 . Comparison of OCT images for the same set of scatterers and the same pronouncedly focused Gaussian bea m simulated using: (a) reference model18 (based on exact analytical solutions obtained in spatial domain for Gaussian s beam s) and (b) the present universal model formulated in spectral domain. (c) shows the depth profiles along the beam axis for the two models (a sterisks are for model18 and circles correspond to the present model). The following example s demonstrate that the ability of the model to describe arbitrary beams opens the possibility to efficiency simulat e various aberrations . Such aberration s may characterize , for example, the eye cornea , which is o f interest for ophthalmic applications of OCT . In the following examples, the aberrations are modeled by introducing phase distortions of an initially Gaussian b eam. Again, to clearer demonstra te how various aberration alter the PSF form, we intentionally took a small number (ten) scatterers located at a depth of the focal plane , which may be considered as imitation of focusing at the eye retina. Figure 3 demonstr ates the simulated scans in (x,y) -plane (C-scans) without (panel (a)) and with aberrations of several characteristic types (Figs 3b, 3c an d 3d). S uch aberrations may arise in OCT - imaging of the eye retina surface through cornea and eye lens with imperfect refraction [ 24], phase distortions of which are often approximated by a set of Zernike polynomials [ 25]). It is clear from Fig. 3 that both worsening of the PSF -locali zation in the focal plane and aberration -induced translational shifts of the entire C -scan can readily be simulated. Such OCT image s with a priori known aberrations can be used for development and comparison of aberration -correction procedures , which is im portant, e.g., for improving the quality of OCT -imaging for ophthalmic application s [24,26,27,28] (a special discussion of utilization of simulated OCT images for such applications will be published later). Figure 3 . Example s of simulated OCT C-scans in ),(yx plane corresponding to the focus depth 128z m for a focused Gaussian beam similar to that, for which ),(zx images in Fig. 2 were shown. Panel (a) with well -localized scatterer images is for the absence of aberrations; (b) is for the aberration corresponding to the Zernike p olynomial 1 3Z (coma); (c) is for the Zernike polynomial 2 2Z (2nd astigmatism) ; (d) is for Zernike polynomial 0 6Z (2nd spherical aberration) . Further, t he ability of the model to easily simula te non -Gaussian beams is demonstrated for a Bessel beam (such beam attract significant interest in OCT due to nearly non -diverging character of propagation29). We imitate d formation of such a beam due to propagation of an initially planar beam with a Gaussian amplitude profile )/ exp(2 2ar through a 8 cylindrical axicone with the opening angle 2 (ref. 30). To this end we introduced in the input plane dzz a radially - dependent factor )]( exp[ tgrikn describing the phase delay acquired by the spectral components of the beam after propagating through the axicone . In the B -scans shown in Fig 4a -1 and 4b -1 we again intentionally used a small number (six) scatterers located at three depths . Figure 4 clearly demonstrate s the difference between the propagation of the simulated Bessel beam ( row (b ) in Fig. 4) and a Gaussian beam with a similar radius of the focal are a (row (a ) in Fig 4 ). The Gaussian beam demonstrates pronounced focusing and strong broadening of the PSF outside the focal depth. In contrast , the Bessel beam demonstrates good conservation of the beam -kernel radius , which is accompanied by slow evolution of the surrounding rings corresponding to the Bessel -function sidelobes (compare the C -scans in rows (a) and (b)). Figure 4. Example s of OCT B -scans in ),(zx plane (first column) and C-scans in ),(yx plane (2nd throu gh 4th colum ns) at three depths simulated for a Gaussian beam (panels (a -1) through (a -4)) and a Bessel beam (panels (b -1) through (b -4)). It is clear that , for the Gaussian beam, the minimal size of the PSF corresponds to the focal depth 248z m, whereas outside the focal depth the PSF sizes are much greater, such that for 452z m, the two scatterers are hardly resolved. In contrast, the Bessel beam at all three depth has nearly invariable size of the central kernel, w hereas the outer ring -shape sidelobes described by the Bessel function gradua lly evolve remaining at least 19 dB lower than the amplitude at the axis. The n ext example s demonstrate the intrinsic ease of the proposed model to simulate 3D data formed by illuminating beams with arbitrary degree of focusing, which can be used for development/verification of methods of numerical refocusing . In addition we demonstrate the ease of simulation of various noises (to imitate the noise at the receiving array). Figures 5a and 5b show examples of the initial OCT image for a highly -focused Gaussian beam before and after numerical refocusing (with the refocusing algorithm used in31,17. The density of scatterers is intentionally chosen fairly low to clearer demonstrate different forms of the point -spread function (PSF) near the focal plane and outside it . It is seen that outside the focal depth the weakly localized scatte rer images are essentially masked by the measurement noise. However, after performing numerical refocusing one obtains nearly identical highly localized P SF-forms independently of the depth, i.e., maximal possible resolution the same as at the focal depth is enables everywhere . Consequently, the scatterers become well visible against the noise even outside the physical -focus depth ( 480z m). Besides simulations of 3D volumes of OCT data for motionless scatterers, t he model allows on e to readily re-simulate OCT images for arbitrarily displaced positions of scatterers . Such possibilities are important for testing algor ithms of OCT -based elastography and angiography in highly controllable conditions with flexibly modifiable parameters . In physical experiments, even using phantoms, it is very challenging to realize such imaging in highly controllable conditions . In this context, t he next examples demonstrate simulations of digital refocusing of 3D OCT -data in the 9 presence of a scatterer flow to imitat e the blood flow in a micro -vessel surrounded by the “solid” tissue. The point is that even a fairly slow motion of scatterers may break conditions of constructive interference of scatterer images during refocusing for various positions of broad scanning beam s with multiple overlapping. Consequently, after refocusing these moving scatterers produce regions of very low intensity highly contrast ed with motionless scatterers. T his effect may be used for singling out moving scatterers against the b ackground si gnal from motionless scatterers . As illustrated by Fig. 5c, the moving and motionless scatterers are indistinguish able before refocusing and Fig. 5d demonstrates that after refocusing the simulated vessel becomes well visible due t o strongly re duced signal level . Notice that near the depths of the physical focus, where the PSF has the minimal size and scatterer images nearly do not overlap for neighboring A - scans, refocusing only slightly changes the OCT image. Consequently, the images of moving scatterers do not change their intensity near the focal depth . By this reason , the flow after refocusing becomes visible as a dark channel only sufficiently far from the focal plane, where the initial size of PSF is much larger than in the focal area. How ever, after refocusing the PSF size become comparable with the focal -spot size everywhere, so that even quite narrow channels with flow can be visualized with a high resolution as illustrated in Fig. 5d. Figure 5. Examples of refocusing of simulated 3D scans for a focused Gaussian beam in the presence of noises and moving scatterers (vertically oriented flow of scatterers ). (a) is a 3D image of a rarified set of scatterers to clearer show the difference of highly localized images near the focal plane located closer to the bottom (the focus depth 480z m and the focus radius is 3 m); notice that above the focal depths , the scatterer images are strongly broadened and are nearly masked by the imitated Gaussian “measurement noise ” corresp onding to SNR~12 dB near the focal plane . (b) is the same image after refocusing, due to which the initially defocused scatterer images above the focus depth became well localized and clearly seen against the noise . (c) is a 3D image for a density of scatte rers closer to density of biological cells in tissues with an axial flow 10 m in diameter imitating a blood capillary which is initially indistinguishable against the background; (d ) the same refocused image, in which the capillary outside the focal depth becomes well visible with the maximal resolution (corresponding to the focus radius) and looks as a dark channel , which is due to flow-related breaking of the constructive interference in the refocused image ; particle displacement between the subsequent A -scans is 410z m, which corresponds to the vertical velocity 3104z m/s for 40 kHz A -scan rate . Finally, the next example relates to elastographic imaging of local s trains in the simulated region, in which a prescribed strain -induced variation in the positions of scatterers is introduced . Strain mapping is the basic step in realization of compression optical coherence elastography (C -OCE)23 which in recent years open ed a broad range of biomedical applications: ela stographic differentiation of tu morous versus normal tissue32,33,34, demonstrations of detailed morphological segmentation of heterogeneous tumor tissues using differen ces in elasticity among various morphological components35,36; studying complex 3D str ain fields o f thermo -mechanical origin in cornea l tissue subjected to laser treatment in the context of refraction correction37; detection of insufficient stability of laser -fabricated cartilaginous implants38; studying deformations of osmotic origin during tissue impregnation by non -isotonic liquids39,23. In view of these very different applications , numerical simulations of exactly controllable arbitrarily complex deformations with possibility to flexibly var y OCT -system parameters, si gnal-to-noise ratios (SNR) , etc. , open s unprecedented p rospects for detailed multi -parameter studies of elastographic processing in OCT. For demonstration, in Fig. 6a we show a simulated 3D st ructural image of a multi -layer structure with various scatteri ng strengths of the layers separated by titled boundar ies. Figure 6b shows the prescribed distribution of interframe strains that define axial displacement of scatterer s between the initial and deformed state s, in which the sample was imaged . The 10 OCT imagi ng is simulated for the same highly -focused illuminating beam as used for the examples shown in Fig. 5 with the focus depth 480z m. The corresponding rather complex distribution of inter -frame phase variations is show n in Fig. 6c. These interframe variations are found in the absence measurement noises other than the strain -induced “decorrelation noise” . However, a very noisy structure of the phase -variation s is clea rly seen above the focal depth in regions of pronouncedly laterally -inhomogeneous strains . In these regions, large groups o f differently displaced scatterers contribute to formation of individual A -scans, so that superposition of waves back -propagated from these scatterers produce very noisy resultant phase variation . The next Fig. 6d shows the phase -variation map obtained after preliminary refocusing , which suppresses this degrading interference of multiple scatterers and yields a much clearer phase -variation map. In ph ase-sensitive OCT, t he axial -strain distribution can be reconstructed b y estimating axial gradient s of interframe phase -variation s using . These gradients can be found either using conventional least-square fitting40, or more advanced “vector” method22,41 that is applied here. T he so-reconstructed strain distribution is shown in Fig. 6e that corresponds to directly found noisy interframe phase -variation Fig. 6c. The strain field reconstructed from the phase - variation shown in Fig. 6c demonstrates strong distortion s/errors in the regions outside the focus, where the phase variations in the regions of strong lateral inhomogeneity looks very noisy in Fig. 6c. In contrast, t he strain map of much better quality shown in Fig. 6f is found using the phase -variation map shown in Fig. 6d after preliminary refocusing. The reconstructed strain map in Fig. 6f is much closer to the exact (assumed) strain distribution (F ig. 6b) with much better quality and uniform lateral resolution over the entire 3D scan. Figure 6. Examples of elastography -related simulations of OCT scans. (a) simulated structural image 192x192x512 microns i n size i n which the brighter lines separate regions subjected to different straining. (b) the assumed distribution of strains according to which displacements of scatterer s were defined (the upper -right inset shows the cross section of the central part of the volume through the plane labeled by the black dashed line. (c) the phase variation between t he initial state and final state corresponding to the strain distribution s hown i n (b). (d) the similar phase -variation map, but obtained with the preliminary performed refocusing of the OCT images in the i nitial and final states. (e) reconstructed strain distribution based on the noisy phase -variation map (c) found without refocusing. (f) reconstructed strain map based on the phase -variation map (d) found with preliminary performed refocusing. The total amount of scatteres is 5105.2 and the Gaussian beam has the focal depth 480z m with the focal -waist radius 3 m like for Fig. 5. (g) schematically shown a real experimental configuration in which thermomechanical strain distribution similar to the simulated example (b) w as experimentally produced (like in works37,42): 1- the OCT probe; 2 - buffer silicone layers between which the studied samples of cartilaginous or corneal tissues (layer 3) were placed; 4 – infrared laser operating at 1.5 m that illuminated the samples and produced thermo -mechanical deformations that were partially irreversible and were used for reshaping the tissue. (h) a real example of interframe strains resembling the simulated strain distribution and visualized by proces sing phase -resolved OCT scans using the vector method41. 11 To demonstrate that strain distributions similar to the simulated ones do occur in real situations, Fig. 6g schematically shows such an experimental situation, in whic h thermo -mechanical strains were produced by irradiating a sample of collagenous tissue (a piece of cartilage or eye cornea) by a beam of infra -red laser. Such studies relate to the development of new methods for non -surgical reshaping of cartilaginous or corneal tissues (see a more detailed discussion in37,38). The buffer silicone layers shown in Fig. 6d serve for performing quantitative compression elastography23 to assess the irradiation -produced variations in the tissue elasticity and make conclusions about the microstructural alterations in the tissue42. Figure 6h demonstrates a real experimental example of depth -distribution of interframe laser -induced thermo -mechanical strains with alternate sign s in a sandwich sample consisting of rabbit cornea and buffer sili cone layers . The strain was reconstructed using the vector method41. Such high -resolution OCT -based elastographic imaging based on approach41 attracts much attention is recent years43,44, so that the possibility to realistically simulate OCT images of deformed tissues should be very helpful for further perfection of the OCT -based strain/ elasticity imaging . 4. Discussion In the previous section we presented a new form of the OCT -scan formation model and gave several instructive examples demonstrating the main possibilities of the proposed model. Probably, its key new feature is that , despite the very compact formulation, the new model allows for rigorous simulations of arbitrary beam profiles (including sharply diaphragmed beams) and does not require th e paraxial approximation. This unique flexibility in utilization of arbitrary phase -amplitude beam profiles suggests very convenient possibilities for modeling arbitrary aberrations (see Fig. 4), which is important in the context of the development/compari son of various methods for their numerical correction. Also various degrees of focusing and qualitatively different types of illumination beams (like Gaussian or Bessel ones in Fig. 3) can efficiently be studied in a unif orm manner using Eqs.(1),(3), (7) a nd (9 ), in which only the beam profile ),,( kyxL should be changed . On no less importance are the possibilities of the described new rigorous model to validate the applicability of previous ly proposed versions of OCT -scan formation models based on various approximations and/or special forms of the beam profile ( such as the neglect of the beam -front curvature , assumption of non-diaphragmed Gaussian beams , paraxial approximation , etc.) The point is that such simplified models can enable incom parable computational efficiency in simulations of large sequences of OCT data . Such large sequences are required, e.g., for studying OCT -scan evolution in the presence of flows or Browni an motion of scatterers, or for the development of numerical methods to compensate masking bulk motions of living tissues in OCT -angiography45,46. For such applications, more advanced models may either require unacceptable in practice computation time (days and weeks) or utilization of unreasonably expensive (and not necessary in many cases) supercomputing facilities. However, the drawback of approximate/ simplified models is that , remaining within the framework of the used approximations , it is impossible to reliably verify the range of parameters, for which the se simpler models are applicable. The above -described model readily allows one to validate the applicability of various simplified model formulation s such as16 or model18 that have already been used for studies of OCT -scan evolution due to variation in the OCT -setup parameters and tissue deformation . For example, they were applied in comparative s tudies of the correlation -based and phase -resolved elastographic processing of OCT -scans21 and perfection of strain and displacement reconstruction using phase -sensitiv e OCT20,41. Now the applicability of such models can be reliably validated in fairly arbitrary conditions. Other possibilities, like intrinsic ease in imitating measurement noises (illustrated in Fig. 5) and intr oducing in the described model the optical -field losses due to absorption and scattering are quite straightforward . By this reason we do not specially focus on these issues. The main limitation of the described , otherwise rigorous , model is the use of single -scattering (ballistic) approximation. It should be pointed out, h owever, that the very principle of OCT -visualization is based on this approximation and the multiple -scattering influence in OCT is only one of various degrading factors . Multiple -scattering effect in many cases are of secondary importance, although even degrading of PSF due to multiple forward scattering by refraction -index inhomogeneities can also be introduced in the described approach by analogy with study15. In fact for the majority of other model formulations, the neglect of multiple scattering is typical except for Monte -Carlo -based models. The latter 12 are intrinsically better suitable for considering “snake -like” photon trajectories incl uding multiple scattering, but computationally the Monte -Carlo approach is known to be very demanding. In this context, rather high computational efficiently of the present model can be specially emphasized. Namely, the main E qs. (1),(3), (7) and (9 ) rely on Fourier transformations that can be performed using efficient algorithms of discrete fast Fourier transform (FFT) . The discrete transversal spectrum ),,( kvug of an arbitrary beam profile ),,( kyxL should be found only once a ccording to Eq. (1) for a chosen mesh on the ),(yx -plane to represent the positions of A -scans. Then for all of the chosen positions ),(0 0yx , the 2D Fourier transformation over wavenumbers ),(vu in Eq. (3) gives the amplitudes ),,(0 0 n s kyxU of the waves incident onto the scatterers. Consequently, all spectral amplitudes ),,(0 0 kyxBs at the receiving array are immediately found ( see Eq.(7) ). Then the computational ly efficient FFT over t he wave -numbers nk and straightforward summation over all scatterers are performed in Eq. (9), which yields the entire 3D OCT -image. The required computation time linearly depends to the number of spectral harmonics in the source spectrum )(nkS (determining the number o f pixels in the axial direction ). Similarly , the linear proportionality is for the total number of A - scan positions in ),(yx -plane (i.e., the number of pixels in lateral directions) and the total number of scatterers . To give an idea about the require d time, simulation of a volume 128x128 x128 pixels (corresponding to 192x192 microns laterally and 512 microns in dept h) with 250000 scatterers (which is close to the density of biologica l cells ~ 5-6 microns in size) requires ~100 min using 10 cores of Intel I -7 6950X CPU. This is about ~twice faster than simulation for the same configuration using the reference model18 based on the exact analytical expressions found for the special case of non - diaphragmed Gaussian beams . In that model spectral amplitudes )(nskB were analytically found and the simulation of A - scans i s reduced to FFT over wave -numbers nk and summation over contributions of all scatterers like in Eqs .(8) and (9). Nevertheless, for simulation of 3D OCT -scans, the present model has demonstrated comparable and even higher computation efficiency than model18, but without the requirement to use non -diaphragmed Gaussian beams . 4. Conclusions In this paper we presented a new form of computationally efficient full -wave model of OCT -scan formation. The model formulation is straightforward and compact with a n unique set of capabilities/features: (i) the illuminating beam may have arbitrary phase -amplitude profile with allowance of a sharp diaphragm and presence of aberrations; (ii) the model does not use paraxial approximation that limits the degree of focusing/divergence of the illuminating beam; (iii) although the broadly used approximation of single (ballistic) scattering by discrete scatterers is assumed , there are no other limitations on the density of scatterers, their distribution in space and sca ttering strengths with possible frequency - dependence; (iv) besides rigorous accounting for focusing/divergence, factors describing the wave decay can readily be introduced by analogy with Monte -Carlo approaches; (v) arbitrary measurement noises also can ea sily be added to simulate required signal -to-noise ratios. The model readily also allows one to re-simulate the data for arbitrary changed scatterer position (e.g., random or flow/deformation -induced) . Thus, within the framework of ballistic scattering by discrete scatterers, the above -listed intrinsic abilities of the model give reasons to characterize it as comprehensive. The unprecedented flexibility and high computational efficacy of the model open unique possibilities for studying OCT -scan properties a nd developing novel diagnostic OCT -based modalities for biomedical diagnostics. Acknowledgements. Basic formulation of the OCT -scan formation model was supported by RFBR grant No 19-02- 00645 ; introducing in the model aberration s and refocusing procedures were developed with support of the R SF grant 17- 72-20249 . Conflict of interest: The authors declare no conflict of interest. Authors’ contributions: ALM, GVG and VYZ derived the background equations of the model ; AAM and GVG proposed modifications to pe rform digital refocusing and introduce aberrations ; ALM and LAM wrote the code and performed simulations ; VYZ, ALM, and AAS prepared the manuscript . 13