Permittivity tensor imaging: modular label-free imaging of 3D dry mass and 3D orientation at high resolution

The dry mass and the orientation of biomolecules can be imaged without a label by measuring their permittivity tensor (PT), which describes how biomolecules affect the phase and polarization of light. Three-dimensional (3D) imaging of PT has been challenging. We present a label-free computational microscopy technique, PT imaging (PTI), for the 3D measurement of PT. PTI encodes the invisible PT into images using oblique illumination, polarization-sensitive detection and volumetric sampling. PT is decoded from the data with a vectorial imaging model and a multi-channel inverse algorithm, assuming uniaxial symmetry in each voxel. We demonstrate high-resolution imaging of PT of isotropic beads, anisotropic glass targets, mouse brain tissue, infected cells and histology slides. PTI outperforms previous label-free imaging techniques such as vector tomography, ptychography and light-field imaging in resolving the 3D orientation and symmetry of organelles, cells and tissue. We provide open-source software and modular hardware to enable the adoption of the method.


DERIVATION OF THE IMAGING MODEL
In this note, we derive the imaging model that expresses the measured intensities in terms of the permittivity tensor of the specimen and the configuration of the microscope.Our model accounts for the following: • Scattering of the vector electric field by the permittivity tensor [1,2] of the specimen.
• Illuminating the specimen with incoherent source and finite illumination aperture, i.e., partially coherent illumination [3,4,5,6].• Polarization of electric field along the imaging axis due to high-NA illumination and imaging.
We assume that the imaging path is free of scalar or vector aberrations.We also assume that the specimens we image are weak scatterers and multiple scattering is negligible.
Some notes on the nomenclature: • We track the dimensions of physical quantities through the development of the imaging model to clarify the properties of the specimen that we recover from the measured intensities.• In optics and the following derivation, the term polarization can imply two distinct quantities: polarization of material and polarization of light.The former refers to the induction of electric dipoles in a material as it interacts with the light, and the latter refers to the orientation of the electric field of light.Unless explicitly stated, the word polarization refers to the polarization of light.
1.1.Permittivity tensor.Permittivity is a measure of how much a material is polarized by an incident electric field.The more easily a material is polarized, the more it delays the electromagnetic wave traveling through it, and the higher the permittivity of the material.Polarization of the material induces electric dipoles, which radiate scattered field.Permittivity of an isotropic material such as water is a scalar, i.e., the induced electric dipole moment aligns and scales linearly with the applied electric field.The polarization of an anisotropic material depends on the angle between the electric field and the symmetry axis of the material.As a result, the strength and the orientation of the induced dipole moments may not align with the direction of the incident electric field.The permittivity of an anisotropic material hence takes a form of a second-order tensor [7,8].
At the spatial scale of visible wavelengths, a material can be considered to be isotropic, uniaxial, or biaxial depending on the symmetry of the molecular distribution [7].Here we focus on biological structures that are uniaxial, i.e., materials that have one anisotropic axis.The symmetry axes of the material are called optic axes.
If a uniaxial material's symmetry axis is aligned with the z-axis, its relative permittivity tensor is a diagonal matrix: where n o and n e are refractive indices experienced by the ordinary and extraordinary wave, respectively.Refractive indices are generally complex numbers, whose real part represents phase delay and imaginary part represents absorption.
Permittivity tensors can be parameterized in several ways: with independent entries of the matrix representation, with eigenvalues and eigenvectors of the matrix representation, or with parameters of the angular distribution of biomolecules.To facilitate physical interpretation of the material properties, we parameterize the permittivity tensor in terms of angular distribution of biomolecules using mean (ε r ) and difference (∆ε r ) of the ordinary and extraordinary permittivity.They are related to the refractive indices as follows: 1 Permittivity tensors have the same dimensions as the scalar permittivity of vacuum, ε 0 .Microscopes are sensitive only to the relative permittivity tensor, which is the ratio of the material's permittivity to the vacuum's permittivity.The relative permittivity and the refractive indices in (1) and ( 2) are dimensionless quantities.
In terms of the mean and differential permittivities, the permittivity tensor in (1) is expressed as A uniaxial material can have arbitrary 3D orientation.We represent the 3D orientation of the permittivity tensor with in-plane orientation, ω, and inclination, θ , as shown in fig.1B.The permittivity tensor in the microscope's frame of reference is computed by applying two rotational transformations [9] to ε r and is expressed as ∆ε r sin 2 θ sin 2ω ∆ε r sin 2θ cos ω ∆ε r sin 2 θ sin 2ω ε r − ∆ε r (cos 2 θ + sin 2 θ cos 2ω) ∆ε r sin 2θ sin ω ∆ε r sin 2θ cos ω ∆ε r sin 2θ sin ω ε r + ∆ε r cos 2θ where T (ω, θ ) is the overall coordinate rotation matrix that is composed of two consecutive rotations around z and y axes as Note that all the variables in the above equations can be a function of 3D space,⃗ r = [x, y, z] T .

1.2.
Vector wave equation and scattering potential tensor.We will now discuss how the permittivity tensor is encoded by the scattered electric field.The scalar diffraction theory [7] has enabled development of the optical diffraction tomography methods that measure complex refractive index (RI).The key insight that enables optical diffraction tomography is that the diffracted field is a 3D Fourier transform of the scattering potential of the material [10].
We extend the scalar scattering potential formalism to scattering potential tensors.
We start by considering the vector wave equation derived from Maxwell's equations as per [1,2] where ⃗ E(⃗ r) = [E x (⃗ r), E y (⃗ r), E z (⃗ r)] T is the electric field in 3D space⃗ r = [x, y, z] T .The dimensions of x, y, and z are [L], where L means length.k 0 = 2π/λ 0 is the free-space wavenumber with the dimension of [L −1 ], and λ 0 is the free-space wavelength of the light with the dimension of [L].The dimensions of the ∇ operator are [L −1 ].
We now consider that the material is surrounded by a medium with scalar permittivity, ε.The vector wave equation in the presence of the surrounding medium is given by where ε rm is the relative permittivity of the surrounding media and is a dimensionless scalar.f (⃗ r) = k 2 0 ε r (⃗ r) − ε rm is the scattering potential tensor with the dimension of [L −2 ] whose components are: When the field propagates in homogeneous media, i.e., ε r = ε rm , (7) reduces to the homogeneous wave equation ∇ × ∇ × ⃗ E inc (⃗ r) = 0, whose solution is the incident propagating wave ⃗ E inc .
The total field in the presence of the specimen is the sum of the incident field and the scattered field ( ⃗ E scat ), which in turn depends on the total field.The total field is expressed as an inhomogeneous vector wave equation: f (⃗ r) ⃗ E(⃗ r) represents the dipoles induced in the specimen by the total electric field.The G(⃗ r) is the Green's tensor that describes the radiation pattern of an anisotropic point scatterer as derived in [11,2].Notice that the scattered field is the 3D convolution of the Green's tensor and the induced dipole distribution.
The Green's tensor is related to the scalar Green's function G(⃗ r) [12,13] that describes radiation due to an isotropic point scatterer.( 10) The scalar Green's function is a spherical wave.Typical imaging systems, however, collect only the forwardpropagating component [12].
The component G pq of the Green's tensor describes the electric field along the p = x, y, z axes due to a dipole oriented along q = x, y, z axes [11].These components are expressed [2,Ch. 10] as follows: The dimensions of G(⃗ r) and components of G(⃗ r) are [L −1 ].In the integral of (9), the dimension of G(⃗ r) f (⃗ r) is [L −3 ], which results in the dimensionless quantity when integrated over the volume.
1.3.Born approximation.The inhomogeneous vector wave equation (9) shows that the interaction between the scattering potential tensor ( f (⃗ r)) and the total electric field ( ⃗ E(⃗ r)) results in the scattered field.In turn, the total field is the sum of the incident and the scattered fields.Thus, there is a recursive functional dependence between the total field and the scattered field.Recursive algorithms [14,15] have been developed to compute the total and scattered fields, but they are computationally expensive.The computation is especially prohibitive when the forward model is used to solve the 3D inverse problem that we have set out to solve.
Here we consider biological specimens that are relatively weakly scattering.The dominant variations in the total field arise from the incident field and the field that has scattered only once.Hence, we can replace ⃗ E(⃗ r) with ⃗ E inc (⃗ r) in ( 9) by assuming the first Born approximation [7] and write the total field as This equation describes how the material interacts with the incident electric field under single scattering conditions.1.4.3D vectorial partially coherent imaging model.In this section, we derive the relationship between the vector field in the image plane and the permittivity tensor when the specimen is illuminated with a partially coherent source.We express the field in the image plane in terms of directly measurable Stokes parameters [3,4].The relationship between the measurable Stokes parameters and the permittivity tensor is nonlinear in general.To develop a computationally efficient inverse algorithm for reconstructing the 3D permittivity tensor, we sacrifice some accuracy by linearizing the model with the approach similar to [16,17,18,19,20].We will show that this linearization leads to a bank of transfer functions that map the components of the scattering potential tensor to the Stokes parameters.We derive 3D transfer functions that describe imaging of 3D specimens that are thicker than the depth of field.At the end, we reduce them to 2D transfer functions that describe imaging of the specimens that are thinner than the depth of field.
1.4.1.Electric field of an incident plane wave.Our polarization microscope measures Stokes parameters of scattered light from the specimen illuminated by a quasi-monochromatic light.Köhler illumination [7] with asymmetric source aperture results in a spatially partially coherent illumination.Each point of the source plane (aperture plane of a illumination condenser lens) generates a plane wave with 3D spatial frequency of T is the transverse spatial frequency, which is proportional to the location of the source point off the imaging axis.ν z is the axial spatial frequency that satisfies the constraint λ is the wavelength of light inside a media with permittivity of ε rm .The dimensions of ν x , ν y , and ν z are The incident plane wave due to a point of the source is where T are the components of the vector electric field.Notice that we have accounted for the axially (z) polarized component that arises due to high-NA illumination.
The electric field at the aperture of the condenser is typically generated with low-NA optics and does not have a component along the z-axis.The electric field at the source plane is: We use high aperture illumination to achieve high 3D resolution, which introduces a z-component in the electric field.The components of the incident electric field are [21]: 1.4.2.Scattered vector field when the specimen is illuminated by a point.To simplify the math, we assume the imaging system has magnification of 1.Hence, the spatial coordinate ⃗ r represents the coordinate in both the specimen space and the image space.
In a Köhler geometry, each point in the source plane (⃗ ν ⊥ ) generates a coherent plane wave that interacts with the specimen and generates scattered field through (12).The camera in the image plane samples a slice of the 3D field that is blurred by the imaging system, i.e., the 3D field is captured by scanning the sample along z axis.The blur due to the imaging system is expressed as a convolution with a 3D coherent point spread function (PSF) h c (⃗ r) of the imaging path, which accounts for the finite aperture of the imaging system and the defocus.
We can now combine the effects of the illumination, scattering, and imaging as follows.
• Compute the incident field due to the source point of choice in the specimen volume using (14).
• Compute the total field in the specimen volume by setting ⃗ E inc = ⃗ E 0 in (12) and computing the scattered field as a convolution of the Green's tensor and the induced dipole distribution.• Filter each component of the specimen field with the 3D coherent PSF to account for the blurring by the microscope.The operations above lead to the following expression for the vector electric field in the image space due to a point at ⃗ ν ⊥ in the source plane: where p = x, y, z, ⊗ denotes 3D convolution.Similar to (9) when the product of the Green's tensor (dimensions [L −1 ]) and the scattering potential tensor (dimensions [L −2 ]) are integrated over the volume, the overall integral of the scattered electric field is dimensionless.
The 3D coherent PSF is: where h 0 (⃗ r ⊥ ) is the 2D in-focus PSF that accounts for the finite aperture of the imaging system, h z (⃗ r ⊥ , z) is the 2D "propagation spread function" that models the blur caused by the propagation of the electric field from defocus z to the image plane, and ⊗ 2D denotes 2D convolution.
To understand the image formation model (15), we take its 3D Fourier transform along⃗ r dimension: where, a quantity x(⃗ u) is the 3D Fourier transform of x(⃗ r).Note that the Fourier transform of the 3D coherent PSF, i.e., hc (⃗ u) is the 3D coherent transfer function (CTF).We visualize this coherent image formation model in Fourier space in Figure 1.It is worth noting that the output electric field has zero bandwidth along the u z direction.This implies that there is no axial sectioning when the image is formed with a single plane wave from a point source (a coherent mode).The axial resolution of PTI is built up by probing the specimen simultaneously with multiple plane waves generated by the partially coherent illumination source, which we discuss next.1.4.3.Scattered vector field when the specimen is illuminated by a source pattern.We just derived the scattered electric field at the detector plane from a coherent point source at ⃗ ν ⊥ .Different partially coherent source patterns (α) are specified by different sets of ⃗ ν ⊥ .We represent the vector electric field in terms of the Stokes parameters as we did in our earlier work [22,23].By integrating the contribution of all the coherent scattered modes from one illumination pattern, we obtain the Stokes parameters [3,4] as Up to this point, we have arrived at a rigorous vectorial single-scattering imaging model termed, vector Born model, which we use to generate all of our simulated data in this paper.Since the Stokes parameters undergo the integration along ⃗ ν ⊥ , the dimension of  (17).The incident field (blue) is a 3D delta function located at the spatial frequency of illumination, ⃗ ν.The scattered field spectrum (red) encodes the information of the scattering potential tensor spectrum ( 1 ⃝) localized on the shell of the Ewald sphere defined by the Green's tensor components( 2⃝) shown at the bottom.The spectrum of the scattering potential tensor is filtered by the spectrum of Green's tensor.Since we only consider forward scattering, this shell is only half of the Ewald sphere [12].Summing the incident field and the scattered field, the total field is also located on the half sphere.The imaging path filters the spectrum of the total field with 3D CTF, hc (⃗ u), which effectively truncates transverse spatial frequencies beyond NA/λ 0 .
1.4.4.Linearized vector Born model.For the subsequent analysis, we take a closer look at the Fourier transform of S pq,α (⃗ r) where ã(⃗ u) denotes the Fourier transform of a function a(⃗ r) at the 3D spatial frequency From this equation, we observe that the spatial frequencies of Stokes parameters that we can measure, Spq,α (⃗ u), are non-linearly related to spatial frequencies of the scattering potential tensor components, fpq (⃗ u).Rigorous estimation of the scattering potential tensor from measured Stokes volumes requires an iterative optimization algorithm as in the classic phase retrieval problem [24], which can be time consuming and subject to convergence instability [25].Instead, we adopt the approach of linearizing the equation for simpler, faster, and more robust solution.
To linearize the relationship between the measured Stokes components and the scattering potential tensor components, we expand (19) into three terms that describe the major contributions of scattered light Spq,α,DC is the constant component of a Stokes parameters, Spq,α (⃗ u), due to the incident electric field.Spq,α,L is the term that linearly relates the scattering potential tensor with the Stokes parameters.The linear term is due to the interference of the scattered electric field and the incident electric field.Spq,α,NL is due to interference between scattered fields, which is the weakest among three terms and non-linearly related to the scattering potential tensor.
When imaging transparent biological specimens, most of the contrast in the Stokes parameters is contributed by the constant and the linear terms.We assume the contribution from the interference of the scattered field with itself is negligible.Neglecting the nonlinear term and subtracting the constant term from Spq,α (⃗ u), we work only with the linear term as After this linearization, we find that the Stokes parameters are a sum of filtered scattering potential tensor components.These filters can be thought of as a filter bank of weak object transfer functions widely used in the scalar wave scattering theory [16,17,18,19,20,26].
We write the scattering potential tensor components in terms of the mean permittivity, differential permittivity, and 3D orientation of the material: The interpretation of the above parameterization of the scattering potential tensor is as follows: ) is a complex function representing density (complex-valued mean permittivity) of the material.The real part (Re {ε r (⃗ r) − ε rm }) delays and the imaginary part (Im {ε r (⃗ r) − ε rm }) absorbs the light.• f 1c (⃗ r), f 1s (⃗ r), f 2c (⃗ r), f 2s (⃗ r), and f 3 (⃗ r) are real-valued functions that encode differential permittivity (∆ε r (⃗ r)), in-plane orientation (ω(⃗ r)), and inclination (θ (⃗ r)).We ignore the effect of diattenuation.Each scattering potential tensor component is linked to the Stokes components via a set of transfer functions.However, the actual Stokes parameters in (18) are a sum or a difference of the Stokes components.
We summarize the linear relationship between the Stokes parameters and the scattering potential tensor in one equation as where S′ m,α (⃗ u) is the DC-subtracted Stokes parameter and Hm,ℓ,α (⃗ u) is the transfer function mapping from each tensor component ℓ to the m-th Stokes parameter under illumination pattern α.The transfer functions are dimensionless.Thus, the dimensions throughout the derivation are consistent because the scattering potential tensor and the Stokes parameters have the same dimensions ([L −2 ]).
The equation above describes the linearized vector Born model that we use for reconstruction of the scattering potential tensor.1.4.5.Reduction to 2D imaging model.When the thickness of a specimen is smaller than the depth of field of a certain imaging condition, acquiring data at one plane by changing the partially coherent illumination patterns enables extraction of 2D phase and anisotropic information at that plane.This 2D information extraction requires a slight modification of the transfer functions.In this 2D imaging case, we consider the scattering potential tensor to be concentrated in a single axial layer and is expressed as f (⃗ r ⊥ ) with the Fourier transform fℓ (⃗ u ⊥ ).
We then substitute this expression for fℓ (⃗ u) in (23) and integrate over the spatial frequency u z to get a corresponding relationship in the 2D imaging case between 2D Stokes parameters and the 2D scattering potential tensor as The transfer functions in 2D imaging case are simply a u z -integration of the previously established 3D transfer functions.

COMPUTING STOKES PARAMETERS FROM INTENSITIES
In fig. 1, we have demonstrated that the raw images from the PTI setup contain intensity modulation induced by the specimen's permittivity tensor.In Supplementary Note 1, we established a imaging model to describe the relationship between components of specimen's permittivity tensor and the Stokes parameters of the output scattered light.Here, we relate the Stokes parameters of the scattered light to the measured intensities.
2.1.Calibration.Ideally, our polarization sensitive camera perfectly detects four channels of linearly polarized light oriented at 0°, 45°, 90°, and 135°.According to the definition of the Stokes vector [27], the measured intensity is linearly related to the Stokes parameters in the following form where ⃗ I is the intensity vector of each polarization measurement, A is the instrument matrix [28] of the optical system, and ⃗ S is Stokes vector of the light.Applying the inverse of the known instrument matrix, A, on the measured intensity vector, ⃗ I, allows us to get the first three Stokes parameters for further deconvolution.However, the actual optical setup has multiple components (e.g.beam splitting prisms) with non-ideal polarization effects.These components are modeled by the non-ideal values in the instrument matrix.Thus, calibration of the instrument matrix is essential.To calibrate the unknown instrument matrix, we adopt the approach in [28,29].We place a linear polarizer mounted on a rotational stage at the specimen plane and rotate it to collect corresponding intensity variations in four polarization channels.We model the intensities collected at angle ω LP of the linear polarizer and retrieve the instrument matrix using Mueller matrix calculus.In addition, we calibrate the polarization state of the incident light by extending the calibration process described by earlier methods.
Assuming the incident light has the Stokes parameters of ⃗ S s = [S s0 , S s1 , S s2 , S s3 ] T , the generated polarization state of light after the linear polarizer aligned with angle ω LP in Stokes parameters is calculated as where S LP,0 (ω LP ) = S s0 + S s1 • cos 2ω LP + S s2 • sin 2ω LP (26) M LP (α) is the Mueller matrix of a linear polarizer aligned at angle ω LP .This formula shows that the total energy of the generated state, S LP,0 (ω LP ), will be a constant along ω LP if the input light is unpolarized or circularly polarized.This total energy can vary if our input light has some linear polarization components.We estimate the total energy of the generated light by summing over intensities of four polarization channel.Through a Fourier series fitting, we obtain the Stokes parameters of the incident light and use it for later deconvolution.
The above polarization state is analyzed by the detection path.The measured intensity is a matrix multiplication of the instrument matrix and the generated Stokes parameters where ⃗ a i is the i + 1-th column of the instrument matrix A. If we normalize the measured intensity with the total energy of the generated light, S LP,0 (ω LP ), we observe that the first three columns of instrument matrix are the Fourier series coefficients of the normalized intensities along ω LP .We use this model to compute the instrument matrix by curve fitting.
With the first three columns of the instrument matrix, we obtain the first three Stokes parameters of the light as where A + denotes a pseudo-inverse of the matrix A.

Background correction.
The Stokes parameters estimated with the instrument matrix calibration may still be affected by the residual polarization effects in the light path that introduce slowly varying background.In addition, the deconvolution starts with the DC-subtracted Stokes parameters according to (6).
We image an empty field of view to characterize the slowly varying background.We correct the background and obtain the DC-subtracted Stokes parameters as follows.
where S bg,i is the i-th background Stokes parameters obtained when imaging at an empty field of view.When the transmission modulation of the specimen is relatively weak, S 0 /S bg,0 ≈ 1 and this is simply removing the DC components of the Stokes parameters.These operations are done to all sets of polarization acquisitions with different illumination patterns and then the deconvolution takes these DC-subtracted Stokes parameters for the retrieval of uniaxial permittivity tensor components.
2.3.Derivation of the Stokes background correction.Since background variation is usually slowly varying, we adopt a non-diffraction-based Mueller calculus to develop our background removal scheme.An anisotropic sample with transmittance t s = |e µ+iφ | 2 , in-plane orientation of ω, out-of-plane orientation of θ , apparent retardance ρ ap (ω, θ ) can be modeled with an Mueller matrix as Illuminating with the circularly polarized light with the Stokes parameter of ⃗ S s = S s0 • [1, 0, 0, 1] T , where S s0 is the energy of the incident light, the output Stokes parameters are obtained as These are the Stokes parameters we get assuming we do not have additional background polarization effects from the optics in the microscope.When the optical path is not perfectly free of polarization effects, slowly varying polarization modulation can mask the weak anisotropy signal.To remove this effect, we model the background polarization effect as an additional transparent retarder with the Mueller matrix of M bg (t bg , ω bg , θ bg , ρ bg ).With this background Mueller matrix, the overall output Stokes parameters are These overall Stokes parameters contain contributions from the specimen and the background.To separate these effects, we acquire the background Stokes parameters ( ⃗ S bg ) by imaging an area without the specimen, Comparing (32) with (33), we can compute corrected Stokes parameters as With this correction, we also remove the DC components of the Stokes parameters in S out,1 and S out,2 .The only remaining Stokes parameters with the DC-component is the S 0 component of ⃗ S c .With the additional subtraction, we arrive at the final DC-subtracted Stokes parameters that we used for the deconvolution We apply these corrections on all sets of intensities from each illumination patterns before deconvolving the data.

INVERSE ALGORITHM
3.1.Solving for scattering potential tensor.According to the imaging model ( 6), the DC-subtracted (background subtracted) Stokes parameters of the scattered light, S′ m,α (⃗ u), under different illuminations are linearly related to the entries of the scattering potential tensors.Thus, we can write down an inverse problem to retrieve the entries of the scattering potential tensor as where τ ℓ is the Tikhonov regularization parameter for each independent term in the scattering potential tensor.In our implementation, this optimization is done in the Fourier space, where each voxel (or pixel) in Fourier space is independent from other voxels (or pixels).

H †
α denotes the Hermitian of the matrix Hα .With this fully determined inverse problem, we apply the Cramer's rule on it to obtain the analytical solution as where H ℓ denotes that the ℓ-th column of the H matrix is replaced with ∑ α H † α ⃗ S′ α .We set the regularization τ ℓ as a multiple of the mean absolute value of the ℓ-th diagonal entry of H over all spatial frequency ⃗ u to reduce its dependency on different imaging parameters.Usually the multiple factor of τ ℓ is ranging from 1 to 100 (smaller for less noisy data to get sharper reconstruction and larger for more noisy data to smooth out the noise).It should also be noted that, for simplicity of the math, we only show the single point solution of this inverse problem.In the actual implementation, we vectorize the whole 3D array of the data and the transfer functions for fast parallel processing with GPU.

3.2.
Mean permittivity, differential permittivity and 3D orientation from scattering potential tensor.After solving (36), we obtain components of scattering potential tensor as shown in (3).We report the 2D or 3D mean and differential permittivities in terms of the angular average of permittivity (ε r − ε rm ) and the angular difference of the permittivity (∆ε r ), respectively.Because of this difference in 2D and 3D imaging, we present solutions for both situations.
We start this step of the reconstruction by estimating the differential permittivity and 3D orientation from the 1c, 1s, 2c, 2s terms in the scattering potential tensor.This step requires information of the optic sign (n e > n o or n e < n o ).Here, we first compute two analytical solutions for differential permittivity and 3D orientation assuming two optic signs.These two sets of solutions will be used to estimate the probability of the optic signs in the next part of the algorithm.Using trigonometry relations, we express these two solutions of differential permittivity and 3D orientation in 2D and 3D imaging modes as 2D : , when when n e ≷ n o 3D : , when when where ω ± , θ ± , and ρ ± are in-plane orientation, inclination, and differential permittivity of the positive or negative uniaxial material.Note that ρ 2D,± (⃗ r ⊥ ) is a 2D image with a unit of a distance at each pixel (difference of optical path length), while ρ 3D,± (⃗ r) is a 3D volume without a unit for each voxel.σ is a small number to prevent noise amplification in the estimation of differential permittivity.Typically, σ = 10 −2 ∼ 10 −3 is a good choice to balance between accuracy and noise amplification.
The 0 term in the scattering potential tensor is related to the difference of average permittivity between the specimen and the environment, which corresponds to accumulated optical phase and absorption information.However, this term is coupled with the differential permittivity, ∆ε r .Removing the contribution of the differential permittivity, we express two solutions of phase and absorption in 2D and 3D imaging modes as where φ ± and µ ± are phase and absorption of the positive or negative uniaxial material.Note that φ 2D,± (⃗ r ⊥ ) is a 2D image with a unit of a distance at each pixel (optical path length), while φ 3D,± (⃗ r) is a 3D volume without a unit for each voxel.
At the end of this computation, the range of the in-plane orientation is ω ∈ [0, π) and the range of the inclination is θ ∈ [0, π).These ranges correspond to the front half of the hemisphere of unit sphere (y ≥ 0).Depending on how we visualize the 3D orientation, we can also transform (ω, θ ) coordinates to span the range of the top hemisphere (Extended Data Fig. 2B, z ≥ 0).

3.3.
Optic sign from scattering potential tensor.The optic sign reports symmetry of underlying structure.If we know of the type of the material imaged, we pick one set of solutions from ω ± , θ ± , and ρ ± .More often, the optic sign of a biological structure is not known and can be spatially variable.
For cases where the optic sign is unknown, we employ an algorithm that models the scattering potential as a mixture of positive uniaxial and negative uniaxial material.The algorithm starts with constructing the scattering potential tensor components with positive (ω + , θ + , ρ + ) and negative (ω − , θ − , ρ − ) set of solutions according to (3).We then want to know which set of solution is more favorable in our data by solving the following optimization problem min where w + (⃗ r) and w − (⃗ r) are the weights for positive and negative uniaxial solutions (we only consider positive values of the weight).When the positive material is more favorable (the structure within a voxel is denser along axis of symmetry), w + (⃗ r) is larger than w − (⃗ r).On the other hand, w − (⃗ r) is larger than w + (⃗ r) when negative material is more favorable (the structure is denser perpendicular to the symmetry axis).When the material is isotropic, w + (⃗ r) ≈ w − (⃗ r).We implement a gradient descent iterative algorithm to solve this optimization.To identify material type with these two weights, we define the probability of being positive material to be where w c is a cut-off weight to remove noisy weight estimate below this value for smooth probability reconstruction.w c is usually set around 0.05 ∼ 0.2.The higher the value, the stronger we suppress the noisy estimate.

IMAGE PROCESSING: MULTI-SCALE ANALYSIS, DENOISING, STRUCTURE TENSOR ANALYSIS
4.1.Multi-scale analysis of optical properties by averaging the permittivity tensor.In fig.4, we describe measurements of optical properties of the mouse brain tissue at multiple scales.Computing physically meaningful measurements of mean permittivity, differential permittivity, 3D orientation, and optic sign at lower spatial resolution from the high-resolution acquisition requires a specific filtering approach.The reconstruction steps that lead to components of the scattering potential tensor (39) are linear.However, the computation of optical properties from the scattering potential tensor components involves non-linear operations (( 40) and ( 43)).Therefore, the low resolution measurements of optical properties are computed by linear filtering of the high resolution volumes of the scattering potential tensor, which are transformed into mean permittivity, differential permittivity, 3D orientation, and optic sign according to (40) and (43).This approach ensures that the computed optical properties at lower resolution are representative of an acquisition with a light path of lower spatial resolution.

4.2.
Improving the sensitivity of PTI through averaging and denoising.PTI converts the intensity modulation from the density and anisotropy of the specimens into physical properties as shown in fig.1C and 1D.These intensity modulations are built upon a constant intensity background that defines the transmission of the surrounding media.The ratio between the strength of the intensity modulation and the shot noise created by the constant background intensity defines the signal to noise ratio (SNR) of our measurements.Biological specimens that we work with generally have strong phase (average RI) and relatively weak anisotropy.Hence, SNR is typically lower for anisotropy measurements.In the following, we take two approaches to improve the SNR of our measurements.
The first approach is to reduce the shot noise in the raw intensity through averaging multiple frames.Since we cannot change the strength of anisotropy from a biological specimen, we average our raw intensity by N avg frames to improve our SNR by N avg times.N avg varies by experiments.For strongly anisotropic structures such as the anisotropic glass target in fig. 3 and myelin sheaths in fig. 4 (differential permittivity ≳ 0.01 with the high-NA objective), N avg = 1 already gives satisfying results.For weaker anisotropic structures such as iPSC-derived cardiomyocytes in fig. 5 and A549 cells in Figure 2 and 3, we average 32 frames and 100 frames per intensity acquisition, respectively, to increase the SNR.
The second approach is to apply an additional wavelet denoising on the differential permittivity after the inverse algorithm (36).Tikhonov regularizer in the existing inverse algorithm prevents over-amplification of noise in the reconstructed physical properties.Total-variation or wavelet-based denoising algorithms that leverage the continuity of the images could provide additional improvements in the SNR of the image.However, a proper implementation of these algorithms with the deconvolution requires an iterative approach to solve the optimization problem, further increasing the currently heavy computation load.We then adopt a compromised solution, directly performing wavelet denoising algorithm on the differential permittivity images, to balance the need of SNR improvement and computation load.The denoising algorithm is a single-step soft-thresholding operation in the wavelet space on the images.This operation produces a solution of the following the optimization formalism [30] min where x is the input image, τ w is the wavelet regularization parameter, and W is the operator of the wavelet transform.The solution of this optimization problem is analytical and is expressed as To denoise the differential permittivity images, τ w is generally chosen to be 10% to 20% of the average signal level.
In addition to low SNR of the differential permittivity measurements, strong variations in mean permittivity can mask the true anisotropy by introducing edge birefringence.A distinguishing feature of the edge birefringence is two maxima [31] around the edge with orthogonal orientations.The form birefringence of a true anisotropic structure (not the edge birefringence) shows a single maxima of a consistent orientation.To emphasize true anisotropy in the differential permittivity channel, we suppress edge retardance by computing an orientation consistency map.Using this map, we reject features with fast-varying orientation.We apply this additional filtering on the cardiomyocyte dataset we show in fig. 5.
The orientation continuity map is computed through the following steps.First, we compute the scattering potential tensor components and filter these components with uniform filters, U(⃗ r), of kernel size Second, we use these averaged scattering potential tensor components to compute the differential permittivity according to (40).The averaged differential permittivity is higher when the orientation is more continuous along spatial dimensions and lower when the orientation is variable.The orientation continuity map is computed from the average differential permittivity normalized by its maximum.Multiplying this map with the original differential permittivity suppresses edge retardance.This method is effective when the orientation of the structures is consistent (fig.5).Potential artifacts may occur when the orientation of the structure varies fast, e.g., lipids in fig.4D).

Computation of 3D orientations of axons and their cross-sections from structure tensor of phase image.
To obtain the structure tensor for analysis of axons' morphological orientation, we modify the vessel segmentation algorithm from the X-ray angiography community [32] for our data shown in fig.4D.The algorithm first computes the structure tensor, the second-order gradient, of the mean permittivity φ (⃗ r), as where s is the scale parameter and scale-dependent gradient is defined as a convolution between the scaled intensity and the derivative of a Gaussian function written as Depending on the geometry of the interested feature (a shell-like, tubular-like, or blob-like structure), the eigenvalues of the structure tensor may be used in different ways to construct the segmentation map.In our case, we selected a shell-like structure and constructed the segmentation map as where are eigenvalues of the structure tensor at position ⃗ r, a and b are parameters to adjust the level of sharpness of the segmentation map.The eigenvectors corresponding to the largest eigenvalue of the structure tensor on every voxel within this map contain the normal orientation of the shell-like structures.We then use these eigenvectors as a different measurement of lipids' 3D orientation to verify our results from 3D anisotropy.

TRANSFER FUNCTIONS AND CHOICE OF ILLUMINATION PATTERN
Our choice of the illumination patterns for our experiment is different from past literature of 3D differential phase contrast (3D DPC) microscopy [20] and 3D polarized light imaging (3D PLI) [33,34,35], informed by the following considerations.
To get 3D phase information, 3D DPC exploits illumination patterns of 4 rotating half-circles and the illumination NA matched with imaging NA, which is an extension from 2D DPC [17].In 3D PLI, 3D orientation and the differential permittivity of the optically anisotropic material (axons in human brain) are extracted with circular pattern, polarization-sensitive detection, and a tilting stage.3D PLI relies on one untilted acquisition followed by four acquisitions with tilts spaced at 90°orientations.The 4 tilted conditions with circular pattern is analogous to keeping the stage level, but illuminating the sample with 4 off-axis illuminations at 90°orientations as has been done in 3D DPC.The one untilted polarization measurement plays an important role for retrieving the 3D orientation retrieval shown in [33].From these observations, we chose to use combination of circular illumination pattern plus several off-axis illumination patterns with a polarization-sensitive imaging path to encode 3D mean permittivity, differential permittivity, and 3D orientation information.
There are a variety of illumination patterns that could be used to encode permittivity tensor.Different illumination patterns result in different shapes of the transfer functions in (6), which affects robustness to noise at different 3D spatial frequencies.By visualizing these transfer functions, we visualize the strength with which the specimen information is transferred for chosen illumination pattern.
Extended Data Fig. 3 compares the absolute sums of the transfer functions over illumination patterns for two possible combinations of illumination patterns that are deduced from past literature.These summed transfer functions describe how individual scattering potential tensor components couple into the measured intensities as a function of their spatial frequencies.In summary, S 0 carries phase (ℓ = 0r) and absorption information (ℓ = 0i), while S 1 and S 2 carry anisotropy information (ℓ = 1c, 1s, 2c, 2s, 3).For both sets of illumination patterns, phase, absorption, and projected anisotropy information (ℓ = 1c, 1s) transfer relatively well compared to the last three components (ℓ = 2c, 2s, 3).The last component of the scattering potential tensor, which contains information of optic sign, is the most poorly-transferred in both cases.This is why the estimation of the optic sign is the noisiest .Comparing the energy of the transfer function with these two sets of illumination, we find that the second set of illumination patterns generates more uniform transfer function for all the scattering potential components.Off-axis illumination with more orientations tends to transfer the specimen information better and, therefore, be more robust to noise across different spatial frequencies.However, larger number of off-axis illumination patterns require longer acquisition.To balance the trade-off between the noise performance and the acquisition time, we choose the second set of illumination patterns for our experiment, which consists of 8 off-axis sector patterns at the highest NA, followed by 1 on-axis circular pattern.Note that these are not the only choices for proper permittivity tensor retrieval.Other combinations of patterns are possible, but it would require further investigation to determine the optimal one as that is done in phase imaging community [36,37].

THE FORMATION OF LASER-INDUCED MODIFICATIONS IN FUSED SILICA
In this note, we discuss the origin of the various density and anisotropy patterns seen in the photolithographically written anisotropic glass targets.These targets are written with the technology used to store digital information in glass, which provides the advantages of high density and long-term stability.
When an intense femtosecond laser pulse is focused into transparent material, e.g., silica glass, high-order nonlinear absorption allows the energy to be deposited predominantly within the focal volume, producing a local permanent refractive-index modification.Although the process of energy absorption is now well understood, the optical properties of the induced structures are an exciting area of research.
Depending on the level of laser intensity, one can induce any of four qualitatively different types of material modification: positive refractive-index change (type 1) at relatively low intensity; birefringent regions with the evidence of nanograting formation (type 2) [38,39] at intermediate intensity and voids (type 3) at high intensity.Recently new type of birefringent modification has been observed as random arrangement of prolate-shaped nanopores at transition from the type 1 to type 2 and coined as type X modification [40].The transition threshold from one kind of structure to another depends on the laser writing parameters such as pulse duration, wavelength and irradiation time.The nanograting (type 2) modification is observed under SEM as periodic subwavelength arrangement of nanoplanes with features as small as 20 nm.The orientation of both nanogratings and elongated nanopores can be controlled by the polarization of the writing laser and results in form birefringence with slow axis oriented perpendicular to the polarization of writing laser and along the nanoplanes or long axis of prolate-shape nanopores.Thus, the optic sign of the induced birefringence is negative with the birefringence for nanograting modification close to the birefringence of the quartz crystal, and the birefringence of nanopore modification is about 4 times lower.
Birefringent structures produced with nanograting modification reveal ultra-low scattering loss and were observed when writing with infrared 1030 nm laser wavelength.Writing with shorter wavelength, such as 515 nm is advantageous for increasing the spatial resolution.With a commercial quantitative polarization microscope, Abrio, only nanograting structures were revealed, since Abrio can only measure the average change of birefringence along axial direction.The structure of imprinted birefringent modification, e. g. the possible coexistence of different types of modification along the axial direction shown in fig.3C, in principle can be characterized by destructive SEM.However non-destructive method of 3D birefringence imaging, such as PTI, is more convenient and resulted in the observation of two distinct modifications across the depth.This data shows that PTI can be used for reading 3D information storged in these glass storage devices.

CYTOPATHIC EFFECTS OF RSV INFECTION ON A549 CELLS
In this note, we use PTI multiplexed with fluorescence deconvolution microscopy to analyze cytopathic effects of another model of infection: A549 cells infected by respiratory syncytial virus (RSV).RSV is a leading cause of serious lower respiratory tract infection in infants and elderly [41].Although it is known that RSV replicates in cytoplasm [42], many aspects of its impact on cellular architecture remain to be studied.As a lung cell line, A549 cells are commonly used to study RSV infection [43].Many studies have immunostained A549 cells to investigate the organelle rearrangement [44,45] during RSV infection.Quantitative measurements of physical properties of RSV-infected cells can clarify the cytopathic effects and enable image-based staging of the infection cycle.
Figure 2A and the corresponding through-focus video Supplementary Video 11 show images of fixed A549 cells with the label-free channels (mean permittivity, differential permittivity, in-plane orientation, out-of-plane tilt, and optic sign probability) and the fluorescence channels.PT measurements reported in this result are acquired with an 1.2 NA (NA c ) oil immersion condenser and 1.2NA (NA o ) water immersion objective.The mean permittivity shows nuclei, nucleoli, and a mesh of membranous structures surrounding nuclei.The differential permittivity identifies anisotropic fibers around the nucleus and the membranous structures near the the cell-cell junctions.One may interpret birefringence at cell-cell junctions to be simply edge birefringence.However, these structures have high differential permittivity relative to the mean permittivity and are not present at every edge of the cells.Therefore, we think that these structures are true anisotropic structures.The A549 cells are stained with DAPI (blue) to label DNA and GFP to label RSV.We see co-localization of density, anisotropy, and fluorescence of nuclei.Since these are the uninfected A549 cells, we observe no GFP signal.
Figure 2B and 2C zoom in on two anisotropic structures indicated by the blue box and the orange box in Figure 2A. Figure 2B shows that the anisotropic structure is positive uniaxial and the 3D orientation of the symmetry axis of PT is aligned with the symmetry axis of the structure throughout the volume.The alignment of the symmetry axes of the microscale PT and the mesoscale structure suggests that this is a polymeric fiber assembly with the bound electrons more easily polarized in the fiber axis, analogous to myofibrils in sarcomeres.From the projected 2D orientation in the side view, we find that the fiber is slightly inclined out of the focal plane.Figure 2C shows the orthogonal sections of the anisotropic structure at the cell-cell junction along its tangential and normal directions.The orthogonal sections suggest this positive uniaxial structure has negative phase and almost 3× stronger differential permittivity than the fibers around the nucleus.From the projected 2D orientation in both views, we see the anisotropy is perpendicular to the symmetry axis.This orthogonality between the symmetry axes of the microscale PT and mesoscale structure suggests this structure is a membrane structure with the bound electrons more polarized perpendicular to the sheet, analogous to lipid bi-layers in the myelin sheaths.A549 cells are epithelial cells that form tight junctions between cells when the confluence of cells is high, which can align the membranes of neighboring cells to give rise to high anisotropy.Now we compare above phenotypes with the phenotypes of the RSV infected A549 cells in Figure 3A and corresponding through-focus videos (Supplementary Video 12).In the fluorescence image, GFP distribution is visible in few cells, indicating they are infected.We observed three consistent phenotypes across multiple cells.First, the highly anisotropic membranes at the cell-cell junctions that are visible in the uninfected cells disappear in the infected cells, suggesting the reorganization of the cellular functions due to RSV infection suppresses formation of these membranes.Second, the perinuclear fibers seen in the uninfected cells are still visible in zoom-in fields of view shown in Figure 3B.Third, the phase or electron density of nuclei and nucleoli in the infected A549 cells are generally larger than the phase of the uninfected A549 cells.Using the DAPI images, we conduct a segmentation of the nuclei for both uninfected and infected A549 cells.Figure 3C shows the histogram of the average and standard deviation of phase value in the nuclei.In the infection condition, we observe noticeably higher average and standard deviation of phase in the nuclei.This suggests that the infection causes the cells to condense themselves.The optic sign probability estimation of the infected condition is noisier than the estimation of the uninfected condition because the cells in the infection condition are thicker and may introduce more unwanted scattering artifacts.Further, we find that phenotypic changes observed   SUPPLEMENTARY FIG. 3. 3D orientation and optic sign probability of an adult mouse brain section in 2D and in 3D: (A) 2D differential permittivity, 3D orientation and (B) optic sign probability images of the whole brain section show key anatomical landmarks.3D orientation is encoded using a colorsphere (Extended Data Fig. 2).At the imaging and illumination NA of 0.55 (spatial resolution of ∼ 0.5 × 0.5 × 3.2 µm), axons behave like negative uniaxial material as seen in the optic sign probability.Therefore, we assume the negative uniaxial material when displaying 3D orientation over the whole slide.Multiple anatomical landmarks are labeled according to the coronal section at level 51 of the Allen brain reference atlas, aco: anterior commissure (olfactory limb), cc: corpus callosum, CP: caudoputamen, CTX: cortex, MO: motor cortex, SS: somatosensory area, VL: ventricle.(C, left) We selected aco area marked with orange boxes labeled with 1 ⃝ in (A) and (B) for 3D high-resolution imaging with imaging NA of 1.47 and illumination NA of 1.4, i.e., with the spatial resolution of ∼ 0.23 × 0.23 × 0.8 µm.The orthogonal sections of 3D differential permittivity, 3D orientation, and the optic sign probability show the complex 3D axon network.(C, right) We further zoom into two volumes with distinct tilt angles: green (axons tilted 60°to the left from z-axis) and red (axons tilted 20°to the left from z-axis) boxes.3D mean permittivity and 3D orientation resolve the boundaries of individual axons, which behave as the positive uniaxial material with 3D orientation perpendicular to the surface of the membrane.The orientation is now rendered using 3D colorsphere and using two distinct colormaps that show the in-plane orientation and out-of-plane tilt to better visualize axon networks.

FIGURE 1 .
FIGURE 1. Illustration of the coherent vector imaging model in Fourier space: The above figure illustrates quantities in(17).The incident field (blue) is a 3D delta function located at the spatial frequency of illumination, ⃗ ν.The scattered field spectrum (red) encodes the information of the scattering potential tensor spectrum ( 1 ⃝) localized on the shell of the Ewald sphere defined by the Green's tensor components( 2  ⃝)  shown at the bottom.The spectrum of the scattering potential tensor is filtered by the spectrum of Green's tensor.Since we only consider forward scattering, this shell is only half of the Ewald sphere[12].Summing the incident field and the scattered field, the total field is also located on the half sphere.The imaging path filters the spectrum of the total field with 3D CTF, hc (⃗ u), which effectively truncates transverse spatial frequencies beyond NA/λ 0 .

FIGURE 2 .
FIGURE 2. Imaging physical and molecular features in healthy A549 cells: (A) 3D mean permittivity, 3D differential permittivity, fluorescence, in-plane orientation, out-of-plane tilt, and optic sign probability images of the healthy A549 cells.The fluorescence overlay shows DAPI stain in magenta and RSV-GFP signal in green.(B) xy slice of mean permittivity and differential permittivity from the volume in blue-square region in (A) shows dense, anisotropic perinuclear structures.We show orthogonal slices (x-y,y-z) of differential permittivity, optic sign probability, and orientation projected onto the plane of viewing.(C) The orthogonal slices of 3D mean permittivity, 3D differential permittivity, projected orientation, and optic sign probability over the field of view indicated by the orange box in (A) show highly anisotropic structures at the cell-cell junction with a positive optic sign.The orange dashed lines cut at two axes (labeled with 1⃝ and 2 ⃝) tangent and normal to the structure to show its cross-sections.

FIGURE 3 .
FIGURE 3. Imaging physical and molecular changes in A549 cells due to infection by respiratory syncytial virus (RSV): (A) 3D mean permittivity, 3D differential permittivity, fluorescence, in-plane orientation, out-of-plane tilt, and optic sign probability images of the A549 cells infected with RSV.(B) Two fields of view chosen from (A) show the perinuclear fibers observed in Figure 2 with label-free channels.We do not observe the highly anisotropic membranes at cell-cell junctions in the infected condition.(C) Histogram of mean and standard deviation of mean permittivity of about 60 nuclei for the uninfected and infected A549 cells show significant increase in density and variability in density of nuclei in infected cells.These statistics agree with the images of mean permittivity shown in Figure 2, and this figure.
in label-free images are present across the field of view, even when the infected cells are only a small fraction of the cell population.This suggests that the infected cells modulate the architecture and activity of neighboring cells through paracrine signaling.SUPPLEMENTARY FIG.2.Experimental verification of the reconstructed tilt angle using a test target: We imaged (A) a flat and (B) a tilted anisotropic glass targets with 20x objective.x-y and x-z sections of 3D differential permittivity and the 3D orientation histogram of the orange ROI are shown.The x-z sections of both targets show that the tilt angles are measurable geometrically.As seen from 3D orientation histograms, the 3D orientation measured with PTI reconstruction matches well with the experimental tilt.