Symmetry constraints for vector scattering and transfer matrices containing evanescent components: energy conservation, reciprocity and time reversal

In this work we study the scattering and transfer matrices for electric fields defined with respect to an angular spectrum of plane waves. For these matrices, we derive the constraints that are enforced by conservation of energy, reciprocity and time reversal symmetry. Notably, we examine the general case of vector fields in three dimensions and allow for evanescent field components. Moreover, we consider fields described by both continuous and discrete angular spectra, the latter being more relevant to practical applications, such as optical scattering experiments. We compare our results to better-known constraints, such as the unitarity of the scattering matrix for far-field modes, and show that previous results follow from our framework as special cases. Finally, we demonstrate our results numerically with a simple example of wave propagation at a planar glass-air interface, including the effects of total internal reflection. Our formalism makes minimal assumptions about the nature of the scattering medium and is thus applicable to a wide range of scattering problems.


I. INTRODUCTION
Scattering and transfer matrices are important mathematical tools that have been applied in a variety of fields, including the transport of electrons in wires [1,2], telecommunications [3], acoustics [4,5] and photonic crystals [6,7]. Within both formalisms, a scattering object or medium is viewed as a 'black box' and the scattering and transfer matrices describe the coupling between different modes in the surrounding regions.
In optics, the scattering matrix, which describes the reflection and transmission of electromagnetic modes by a system, is particularly powerful in the study of disordered media for which the underlying scattering potential is practically unknowable, but for which the reflected and transmitted fields are experimentally measurable [8]. In many modern scattering experiments, the degrees of freedom of an electromagnetic field are explored through wavefront shaping using spatial light modulators, which has facilitated experimental measurements of scattering matrices of complex media for both scalar and vectorial light [9][10][11]. Knowledge of these matrices and their statistical properties has revealed the existence of highly transmitting 'open eigenchannels', even for optically thick systems [12][13][14], and has greatly enhanced the prospect of imaging through multiple scattering media, such as biological tissue, through careful wavefront control and the exploitation of scattered field correlations, such as the memory effect [15][16][17][18][19][20][21]. Scattering matrices have also been employed in theoretical studies of random laser modes [22], information transfer through random media [23], coherent backscattering [24] and Anderson localization [25].
Alternatively, the transfer matrix describes coupling between modes on either side of a scattering medium, * Corresponding author: matthew.foreman@imperial.ac.uk provided that the geometry of the system permits a meaningful identification of two opposite sides. Examples of such systems include optical waveguides [26][27][28] and stratified media consisting of a series of contiguous slabs [29,30], such as superlattices [31,32] and multilayer thin-films [33,34]. The primary advantage of transfer matrices over scattering matrices for layered systems is that a system's overall transfer matrix can be easily computed by taking the correctly-ordered matrix product of the transfer matrices of each separate layer [35]. In contrast, the corresponding composition law for scattering matrices is more computationally demanding [36] and, in some instances, requires unphysical assumptions, such as neglecting multiple reflections between surfaces within a stratified medium [37]. This property makes transfer matrices particularly well suited to numerical simulations [38][39][40] and theoretical studies, such as in photonic band structures [41] and mesoscopic scattering [42,43].
Traditionally, scattering and transfer matrices are defined only for modes that propagate to the far field. The exact values of the elements of these matrices depend strongly on the type of system being considered and may vary significantly from one example to another. Nevertheless, under rather general conditions, both matrices can be shown to obey certain mathematical constraints valid for large classes of scattering media. For example, it has long been known that a system that conserves energy, i.e. does not absorb or generate light, regardless of its microscopic configuration, must possess a scattering matrix that is unitary [44]. More recently, techniques such as scanning near-field optical microscopy have enabled studies in which the scattering of evanescent fields play a critical role, such as in single molecule near-field imaging [45], scattering from plasmonic nanoantennas [46], particle tracking [47] and near-field speckle imaging [48]. Such studies have motivated the introduction of an extended version of the scattering matrix that is capable of describing scattering to and from evanescent field components. In a notable paper by Carminati et al., the mathemat-ical constraints obeyed by the extended scattering matrix were explored for scalar waves under the conditions of conservation of energy, reciprocity and time reversal symmetry [49]. The corresponding scattering matrix constraints due to reciprocity for vector evanescent waves has also been considered separately [50].
Matrix constraints such as those mentioned place limits on the set of all physically possible scattering and transfer matrices, and hence can serve as useful guides in determining whether a given experimental or simulated matrix satisfies the corresponding physical law. In addition, these constraints may also be useful in designing matrix-based models and simulations for scattering in complex media. Furthermore, random matrix theory, in which the set of constraints satisfied by a matrix is generally the only assumption made, has proven to be very fruitful at uncovering universal properties of random scattering media [23,44,51]. We therefore believe that an accurate knowledge of the constraints satisfied by both the scattering and transfer matrices is important for future studies of scattering problems.
Compared to the scattering matrix, theoretical analysis of the transfer matrix seems to have received less attention in the optics literature. Moreover, while previous works have explored matrix constraints pertinent to a continuous decomposition of an electric field containing an infinite set of modes, such as in a continuous angular spectrum decomposition [52,53], the corresponding constraints satisfied by the scattering and transfer matrices defined with respect to a finite set of modes are equally important, particularly for experiments and simulations in which only a finite description is physically possible. The purpose of this paper is therefore to present a selfcontained derivation of the set of constraints imposed upon the scattering and transfer matrices by conservation of energy, reciprocity and time reversal symmetry. Notably, our treatment takes full account of the vector nature of light and allows for evanescent components. A treatment of vector fields is essential in optics, as scattering typically gives rise to polarization mixing and depolarization, neither of which can be described within a scalar wave formalism. We describe arbitrary fields using a vector angular spectrum decomposition, first over a continuous range encompassing all possible wavevectors and then for a discrete angular spectrum containing a finite set of modes. The angular spectrum decomposition is particularly useful as it is able to conveniently discriminate between propagating and evanescent modes. This work builds upon previous results and, for the sake of literary continuity, we adopt similar notation and follow a similar format to that of Ref. [49].
In Section II, we define the scattering and transfer matrices for a continuous angular spectrum of an electric field and derive the constraints imposed by the aforementioned conditions. In Section III, we introduce the scattering and transfer matrices for a discrete angular spectrum and derive the associated constraints from the continuous case. We then compare our results to previ-ously reported results and show that the latter follow as special cases. We end Section III by giving a simple numerical example of wave propagation at a glass-air planar interface, including the effects of total internal reflection. Finally, in Section IV, we summarise and conclude our work.

A. Preliminaries
We consider the scattering problem depicted in Figure 1. A dielectric scattering medium is situated within the region −l ≤ z ≤ l. We denote by R − and R + the regions −L < z < −l and l < z < L surrounding the scattering medium, which are assumed to be dielectric with constant permittivity 1 . All sources are contained in the regions z < −L and z > L, and may produce both propagating and evanescent incident fields. We assume that the scattering medium is linear and all fields are monochromatic with angular frequency ω. The scattering medium can be described by a spatially inhomogeneous, complex-valued permittivity function (r, ω), where r = (x, y, z) T is the position vector and the superscript T is used to denote the transpose of a vector or matrix. We also assume that both the scattering and background media are non-magnetic and have magnetic permeabilities equal to the vacuum permeability µ 0 . We denote by E(r) and H(r) the complex phasor representations of the real electric and magnetic fields with a suppressed time factor of exp(−iωt).
The frequency-domain Maxwell equations for the entire region −L < z < L are given by [54] ∇ · E(r) = 0, From Eqs. (3) and (4), the vector wave equation for the electric field can be shown to be where k 2 = ω 2 µ 0 (r, ω). In the regions R − and R + , we may write the electric field as a sum of plane wave components using the well-known angular spectrum representation [55] E ± (r) = a ± (κ)e i(κ·ρ+kzz) dκ where the plus or minus sign is chosen according to the region under consideration and dκ = dk x dk y . In Eq.  For each plane wave, the z component of the associated wavevector k z = k z (κ) is given by The vector a ± (κ) denotes the amplitude of the righttravelling plane wave with wavevector k = (k x , k y , k z ) T . Similarly, the vector b ± (κ) denotes the amplitude of the left-travelling plane wave with wavevector k = (k x , k y , −k z ) T . In Eq. (6) and all integrals that follow, the domain of integration is assumed to be from −∞ to ∞ for all integration variables unless specified otherwise. The corresponding angular spectrum representation for the magnetic field can be obtained by taking the curl of Eq. (6) and using Eq. (3). It is also necessary that the electric field satisfies the divergence condition in Eq. (1). By taking the divergence of Eq. (6) and using Eq. (1), we find that for all κ.
In R − and R + , k 2 is a constant since k 2 = n 2 k 2 0 = n 2 (2π/λ 0 ) 2 , where λ 0 is the wavelength in vacuum, n = 1 / 0 is the refractive index and 0 is the vacuum permittivity. Clearly when |κ| 2 ≤ k 2 , k z is real and the plane waves are homogeneous, or propagating. When |κ| 2 > k 2 , k z is imaginary and the plane waves are inhomogeneous, or evanescent. For convenience, we define the sets Γ p and Γ e , where Γ p = {κ : |κ| 2 ≤ k 2 } is the set of all transverse wavevectors corresponding to propagating plane waves and Γ e = {κ : |κ| 2 > k 2 } is the set of all transverse wavevectors corresponding to evanescent plane waves.
For linear scattering, we may relate the amplitudes of plane waves travelling towards and away from the scattering medium using the scattering matrix S. Specifically, if we define the column vectors I(κ) = [a − (κ), b + (κ)] T and O(κ) = [b − (κ), a + (κ)] T , then the continuous scattering matrix is defined to be the matrix that satisfies For a given κ and κ , the scattering matrix S(κ, κ ) is a 6 × 6 matrix of complex entries, but it is useful to write it as a 2 × 2 block matrix in the form where r, r , t and t are 3 × 3 matrix generalizations of transmission and reflection coefficients. An alternative description of the scattering problem is possible using the transfer matrix M, which relates the amplitudes of plane waves on the left and right hand side of the medium. Letting L(κ) = [a − (κ), b − (κ)] T and R(κ) = [a + (κ), b + (κ)] T , the continuous transfer matrix is defined to be the matrix that satisfies As with the scattering matrix, it is useful to write the transfer matrix as a 2 × 2 block matrix in the form Unlike the scattering matrix, however, the sub-matrices α, β, γ and δ, do not have such an obvious physical interpretation.

B. Conservation of Energy
We now consider the constraints placed upon the continuous scattering and transfer matrices by energy conservation. In order to enforce conservation of energy, we consider the time-averaged Poynting vector associated with the fields, which is given by The average net rate W at which electromagnetic energy flows out of any closed surface A is given by wheren is the outward unit normal vector to the surface. If energy is not absorbed within the volume enclosed by the surface, then conservation of energy demands that W = 0 and the integral in Eq. (14) must vanish. We shall henceforth assume that this is the case. Consider now the surface ∂V shown in Figure 1 bounding the volume V . The surface consists of several parts: two planar sections at z = z − and z = z + , and a spherical section of radius R in between the two planar sections. In the limit R → ∞, the two planar sections expand to infinite planes and the energy flow through the spherical section becomes negligible. The total energy flow through the planes can be written as (15) where dρ = dxdy. The vector S can be expressed in terms of the angular spectrums of the electric and magnetic fields defined previously. Doing so and evaluating the integrals in Eq. (15) yields where the dependence of each term on κ has been temporarily omitted for brevity.
To derive relations pertaining to the scattering matrix, we notice that Eq. (16) can be recast in the form Introducing the scattering matrix into Eq. (17) using Eq. (10) and simplifying the resultant expressions, we ultimately arrive at where δ is the Dirac delta function, the superscript † denotes the conjugate transpose and we use I n to denote the n × n identity matrix. These so-called extended unitarity relations are hence an expression of conservation of energy and generalize the better known unitarity condition on the scattering matrix so as to include evanescent wave components [49].
The transfer matrix can also be shown to obey a similar set of equations to Eq. (18). Returning to Eq. (16), we note that the integrands can be written in the alternate form where we have introduced the generalized Pauli matrices where O n denotes the n × n zero matrix. Similarly to before, we introduce the transfer matrix into Eq. (19) using Eq. (12), which yields otherwise.
Eq. (21) therefore represents the constraint imposed upon the transfer matrix by energy conservation.

C. Reciprocity
In this section, we consider the constraints imposed upon the scattering and transfer matrices due to the reciprocity principle. Roughly speaking, reciprocity describes a relation between scattering matrices that are related by an exchange of input and output modes. For a more detailed review, see Ref. [56].
Consider again the volume V shown in Figure 1 and its bounding surface ∂V . Let E 1 and E 2 be two arbitrary fields that satisfy Maxwell's equations. Applying the vector analogue of Green's second identity to these two fields gives [57] Since both fields satisfy Eq. (5) in V , it can be shown that the integral on the left hand side of Eq. (22) vanishes. Furthermore, by writing the electric fields in the far field as sums of incoming and outgoing spherical waves (see e.g. Ref. [54]), it can be shown that the integral on the right hand side of Eq. (22) over the spherical section of ∂V vanishes in the limit R → ∞. Therefore, we conclude that the integral on the right hand side of Eq. (22) over the infinite planes z = z − and z = z + is equal to zero. Expressing E 1 and E 2 using angular spectrum representations, we arrive at the equation I z − = I z + , where Written in terms of the vectors I and O, the equation By now introducing the scattering matrix into Eq. (24) using Eq. (10), we are able to derive the reciprocity relation for the scattering matrix which is valid for all pairs of transverse wavevectors κ and κ . Reciprocity relations for the constituent transmission and reflection matrices can be obtained by considering the block form of the scattering matrix as in Eq. (11). Comparing blocks on either side of Eq. (25), we obtain which are consistent with those previously reported [50]. Alternatively, writing Eq. (23) in terms of the vectors L and R, we find which, after introducing the transfer matrix using Eq. (12), gives Eq. (30) is thus the reciprocity relation for the transfer matrix.

D. Time Reversal Symmetry
Let us now consider time reversal symmetry, which describes an invariance under the transformation t → −t. Let E 1 be an arbitrary, real electric field and let E 2 be its time reversed counterpart, i.e. E 2 (r, t) = E 1 (r, −t). It can be shown that the Fourier transforms of the fields satisfy E 2 (r, ω) = E * 1 (r, ω) [58]. Making use of this property and comparing the angular spectrums of the two fields, we find for κ ∈ Γ p and for κ ∈ Γ e . In terms of the scattering matrix, a system is time reversal invariant if both E 1 and E 2 satisfy Eq. (10) for the same scattering matrix S. Combining this property with Eqs. (31)-(34), we ultimately find which is therefore an expression of time reversal symmetry for the scattering matrix.
If we instead require that both fields satisfy Eq. (12) for the same transfer matrix, we find where Eq. (36) is therefore the time reversal symmetry constraint for the transfer matrix.
To conclude this section, we note briefly that for the scattering matrix the time reversal symmetry equations also follow from the conservation of energy and reciprocity equations, i.e. Eq. (18) and Eq. (25) imply Eq. (35). Similarly, for the transfer matrix, the reciprocity relation can be derived from conservation of energy and time reversal symmetry, i.e. Eq. (21) and Eq. (36) imply Eq. (30). We therefore conclude that, in accordance with Ref. [49], for a system that conserves energy, reciprocity and time reversal symmetry are equivalent.

III. THE SCATTERING AND TRANSFER MATRICES FOR A DISCRETE ANGULAR SPECTRUM
The scattering and transfer matrices we have considered so far are, in principle, defined for all pairs of conceivable transverse wavevectors κ and κ , which form a continuous spectrum and are infinite in number. In both experiments and numerical simulations, however, fields cannot be resolved with infinite precision and must be described using some finite set of modes [59]. Beyond merely being a practical limitation, the number of independent modes a system can support in reality must be finite due to the wave nature of light and is often constrained by the geometry of the scattering system. This is particularly relevant for waveguides, such as optical fibers, where the number of modes is finite and is determined by the fiber's radius and refractive indices [60]. Even for waves in free space, however, diffraction places a lower limit on the resolution to which a field can be discretely sampled [61]. Sampling beyond this limit would result in a set of modes that would not be independent and would therefore not yield additional information.
In this section, we shall consider a pixel-wise discretization of the continuous spectrum of transverse wavevectors. This will allow us to define scattering and transfer matrices that describe coupling between plane waves in a discrete angular spectrum. We shall then derive the constraints that must be satisfied by these matrices.

A. Definitions
Let us consider a finite set of plane waves indexed by their transverse wavevectors. We first define two sets, K p and K e , of transverse wavevectors for propagating and evanescent waves respectively. We form K p by choosing N p transverse wavevectors κ p i corresponding to propagating waves together with their additive inverses −κ p i and possibly the zero vector. As we shall demonstrate, it is necessary to include modes in inverse pairs to fully explore the effects of reciprocity and time reversal symmetry. Similarly, we construct K e by taking N e transverse wavevectors κ e i corresponding to evanescent waves and their additive inverses −κ e i . Thus, we have which, with the inclusion of the zero vector in K p , contain 2N p + 1 and 2N e elements respectively. The set of all wavevectors K is then the union of these sets, i.e. K = K p ∪ K e . A sample of several modes is shown in Figure 2. Note that while in Figure 2 we have, for simplicity, distributed the modes at points on a rectangular lattice in k-space, the choice of other geometries, such as a hexagonal lattice, may have practical advantages [62]. Note also that although k-space is unbounded, there is a practical upper limit to the size of |κ e i |. This is because when |κ e i | is sufficiently large, its corresponding wave amplitude, even at positions very close to the scattering medium, will have decayed to the point of being practically unmeasurable.
As a result of Eqs. (8) and (9), the amplitude associated with each plane wave has only two degrees of freedom. Therefore, in order to remove extraneous information and simplify the ensuing mathematics, it is useful to introduce the vectors Note that when κ ∈ Γ p , k z is real and the vectors e k , e φ and e θ are the standard unit basis vectors in spherical polar coordinates. When κ ∈ Γ e , e k and e θ become complex vectors and the latter may be interpreted using a complex polar angle. The vectors e φ and e θ are also classically referred to as s and p modes in polarisation theory [55], but we have chosen to reserve the letter p in this work for 'propagating'. By definition, it follows that which means we may express a ± (κ) and b ± (κ) in Eq. (6) in terms of their θ and φ components. Explicitly, we have It is now possible to reduce the matrices r, r , t, t , α, β, γ and δ in Eqs. (11) and (13) to 2 × 2 matrices that couple the θ and φ components of a ± (κ) and b ± (κ). As an example, consider r(κ i , κ j ), where κ i and κ j are any two vectors taken from K. This matrix describes the reflection at the left side of the medium of an incident wave with wavevector (κ j , k zj ) T and amplitude a − (κ j ) to a final wave with wavevector (κ i , −k zi ) T and amplitude b − (κ i ), where k zj = k z (κ j ). We define the 2 × 2 reduced reflection matrix for the pair of modes κ i and κ j as the matrix r (i,j) , where and where m and n may be chosen to be either θ or φ. Note that for all κ i and κ j , but, since a − (κ j ) has no k component, r(κ i , κ j )e k (κ j , k zj ) is undefined. We therefore assign which justifies the construction of the reduced reflection matrix as only the four components in Eq. (47) are unconstrained. The other sub-matrices of S and M can be treated similarly and a brief summary is given in Appendix A.
We can now construct the discrete scattering matrix for our finite set of modes by considering a discretized version of Eq. (10). To achieve this, we partition k-space into a series of regions, each of which is centred on a mode in the set K. We denote the area of these regions by ∆κ, which, for convenience, we assume is the same for each region. If, for example, a rectangular partitioning is used, then ∆κ = ∆k x ∆k y , where ∆k x and ∆k y are the distances in k-space between adjacent modes. We assume that the choice of modes and k-space partitioning are such that the scattering and transfer matrices are approximately constant over each region.
Let κ i ∈ K be any transverse wavevector. We may replace the integral in Eq. (10) with a sum and write where κ j in the sum ranges over all modes in K. By now letting κ i vary over the set K, we obtain a system of equations, each of which have the same form as Eq. (51), which can be combined into a single matrix equation. To facilitate this, we first introduce the notation where u denotes either a or b and q is either p or e. Regardless of the choice of q, we order the transverse wavevectors within u from left to right in the same way as they are presented in Eqs. (38) and (39). Using the notation in Eq. (52), we define the four vectors with which we define the discrete scattering matrix S as the matrix satisfying Given our ordering of the modes, we may partition the matrix S in the block form where S qq = r qq t qq t qq r qq (59) and q and q are either p or e. For each sub-matrix, the right subscript denotes the type of incident mode (i.e. propagating or evanescent) and the left subscript denotes the type of outgoing mode. For example, r pe describes the reflection at the left hand side of the system of incoming evanescent modes to outgoing propagating modes. It is formed by concatenating 2 × 2 reduced reflection matrices of the form in Eq. (47). All other sub-matrices of S can be understood in an analogous manner. Similarly, we define the discrete transfer matrix M to be the matrix that satisfies This matrix can also be partitioned in an analogous way to the discrete scattering matrix. Explicitly, In order to simplify the equations in the following section, it is useful to normalize the scattering and transfer matrices. We use a bar to indicate the normalized version of a matrix and define the normalized continuous scattering and transfer matrices bȳ The corresponding normalized discrete scattering and transfer matrices are given bȳ where In the case that ∆κ is different for different k-space regions, we can instead incorporate its different values into the matrices η ± 1 2 . To end this section, we note that it is possible to convert between the discrete scattering and transfer matrices. It can be shown from Eqs. (57) and (60) that where, for example, and the other matrices are defined analogously. These equations also hold for the normalized scattering and transfer matrices.

B. Conservation of Energy
We now aim to derive the constraints imposed uponS andM by energy conservation. We begin by discretizing the conservation of energy equation for the continuous scattering matrix, namely Eq. (18). As before, we assume that the continuous scattering matrix is constant over each k-space region. The Dirac delta function δ(κ − κ ), whose integral is by definition unity, can be replaced by the normalized Kronecker delta δ ij /∆κ. Furthermore, by using the normalized scattering matrix, it is possible to remove all k z terms from the equation. Note that a factor of i is introduced whenever a k z factor corresponds to an evanescent wave. Replacing the integral with a sum and rewriting κ, κ and κ as κ l , κ i and κ j respectively, we obtain (74) We can consider each of the four cases in Eq. (74) separately. In each case, we may form four matrix equations by considering different sub-matrix blocks individually. Since there are a lot of equations and the algebra involved is rather lengthy and repetitive, we shall only present a single example here. Suppose κ i , κ j ∈ K p and we consider the first case of Eq. (74). Equating the top-left blocks of the matrices on either side, we obtain We may further extract four equations from Eq. (75) by pre-multiplying and post-multiplying both sides by different combinations of e θ and e φ . We also make use of the fact that, since e k (κ, ±k z ), e θ (κ, ±k z ) and e φ (κ, ±k z ) form an orthonormal basis of C 3 , we have [63] e k e T k + e θ e T θ + e φ e T φ = I 3 .
This means that multiplying any matrix by the combination of vectors on the left hand side of Eq. (76) will leave the matrix unchanged. Inserting this combination of vectors in between the matrix products in Eq.
(75), pre-multiply the equation by e T θ (κ i , k zi ) and postmultiplying the equation by e θ (κ j , k zj ) yields Note that all terms involving either e k or e T k vanish, as previously discussed. The left hand side of Eq. (77) consists of eight bilinear forms, each of which can be simplified separately. For example, the first part of the first term in the sum in Eq. (77) can be shown to be where we use the fact that e θ = e * θ for propagating waves. Simplifying all terms in Eq. (77) in a similar way gives κ l ∈Kp r * (l,i)θθr(l,j)θθ +r * (l,i)φθr(l,j)φθ +t * (l,i)θθt(l,j)θθ +t * (l,i)φθt(l,j)φθ = δ ij Repeating this process for all combinations of e θ and e φ and all sub-blocks in Eq. (74) yields a large system of equations that can be shown to be equivalent to the single matrix equationS † which is the classic result that a scattering matrix for a system that conserves energy and only considers propagating modes is unitary. The other three cases for κ i and κ j can be treated similarly. After a lot of algebra, we obtain the equations By introducing the matrices I p = diag(I 4Np+2 , O 4Ne ) and I e = diag(O 4Np+2 , I 4Ne ), we can combine Eqs. (80)-(82) into the single equation which is the most general form of the conservation of energy constraint for the scattering matrix, incorporating both vector properties of the electric field as well as evanescent components. We note that our result here is consistent with a previously reported equation [64].
To derive the conservation of energy constraint for the discrete transfer matrix, we begin by discretizing Eq. (21) to obtain otherwise.
Performing similar steps to those demonstrated for the scattering matrix, we find, again after a lot of algebra, the equations where the zero matrix in Eq. (86) is of size (4N p + 2) × 4N e . Eqs. (85)-(87) can be combined into the single equationM † ΩM = Ω, where Ω = diag(Σ z 4Np+2 , Σ y 4Ne ). Eq. (88) therefore represents conservation of energy for the discrete transfer matrix.

C. Reciprocity/Time Reversal Symmetry
As noted previously, when conservation of energy holds, reciprocity and time reversal symmetry are equivalent. In this section we shall therefore henceforth refer to both conditions as 'reciprocity'. We begin by deriving the reciprocity constraint for the discrete scattering matrix. Note that for the continuous scattering matrix the reciprocity relation Eq. (25) is considerably simpler than the time reversal symmetry relation Eq. (35), particularly as it does not involve an integral. We therefore begin our derivation with Eq. (25).
In terms of the normalized scattering matrix, the reciprocity constraint of Eq. (25) has two different cases: Let us consider first the case where κ i , κ j ∈ K p . As in the previous section, we again consider each sub-matrix ofS separately. Comparing top-left blocks, we havē If we pre-multiply Eq. (89) by e T θ (κ i , −k zi ) and postmultiply by e θ (κ j , k zj ), the left hand side becomes r (i,j)θθ and the right hand side, after a bit of manipulation, becomes r (−j,−i)θθ where here the subscript −j refers to the mode with transverse wavevector −κ j . Repeating this for all four combinations of e θ and e φ , we obtain the relations r (i,j)θθ =r (−j,−i)θθ , ,r (i,j)θφ = −r (−j,−i)φθ , (90) Eqs. (90) and (91) are equivalent to the single matrix relationr where we have introduced the reciprocal operator R, which we define such that if [A] mn is the (m, n) element of the matrix A, then This particular symmetry of the reflection matrix is a well-known result in scattering theory for polarized light and is sometimes referred to as the backscattering theorem [65]. The operator R has also been discussed previously in the context of reciprocal Jones matrices [66].
By now carefully considering the structure of the matrixr pp , it follows from Eq. (92) that where is the matrix containing 2N p +1 copies of the 2×2 identity matrix on its anti-diagonal and zeroes elsewhere. The effect of multiplying a matrix on either side by σ p is to reflect the positions of all 2 × 2 sub-matrices horizontally and vertically about the central rows and columns of the matrix, but to leave the sub-matrices themselves unchanged. This is necessary so that Eq. (94) correctly equates sub-matrices of r pp that are related by an inversion of transverse wavevectors. Similarly, by considering the other sub-matrices ofS, we find which can be combined into the single equation where σ pp = diag(σ p , σ p ). If we similarly introduce σ e as the matrix containing 2N e copies of I 2 on its antidiagonal and zeroes elsewhere, and σ ee = diag(σ e , σ e ), we can further derivē Finally, Eqs. (98)-(100) can be combined into the single equationS where ω = diag(σ pp , iσ ee ). Eq. (101) is the reciprocity constraint for the scattering matrix, which generalizes the well known equation S = S T . To derive the reciprocity constraint for the transfer matrix, we begin with the time reversal invariance constraint Eq. (36), which is notably simpler than the reciprocity constraint in Eq. (30), and perform analogous steps to those used in the derivation for the scattering matrix. We eventually arrive at the equations where Eqs. (102)-(103) can be combined into the single equa-tionM where ω = diag(σ pp , σ ee ). Eq. (105) is hence the reciprocity constraint for the discrete transfer matrix. Eqs. (83), (88), (101) and (105) constitute the main results of our work and must be satisfied for any system that conserves energy and is reciprocal/time reversal invariant. We note that all four of these equations are algebraic and are therefore much simpler to work with than the integral constraints for the continuous scattering and transfer matrices.

D. Comparison With Previous Results
Here we briefly compare our results with those that already exist in the literature. In particular, we show that the more commonly presented equations follow from our results as special cases. From Refs [8], [42], [44], [67] and others we see that, for a system that conserves energy and and is reciprocal/time reversal invariant, the matricesS andM obey the equations where I, Σ x and Σ z are of the appropriate size. There are three important differences between the matrices in Eqs. (106)-(109) and those of our results. Firstly, the scattering and transfer matrices in Eqs. (106)-(109) are defined only for propagating modes in the far field. Therefore, in comparing with our results, it is sufficient to only considerS pp andM pp . As previously noted, Eq. (80) is identical to the unitarity condition of Eq. (106). For the transfer matrix, since in the far field we haveM ep = O, we see that Eq. (85) reduces to the same form as of Eq. (108).
The second important difference is that Eqs. (106)-(109) are defined only for scalar waves. A scalar wave formalism is appropriate when there is no change in polarisation state induced by the scattering medium. For this to be the case within our formalism, it is necessary for each 2 × 2 sub-matrix withinS andM to have zero off-diagonal elements. Mathematically, this means that [S] mn = [M] mn = 0 when m + n is an odd number. By now inspecting Eq. (93) we see that, if this is the case, then the reciprocal operator is identical to a regular matrix transpose and in all of our previous equations we may make the notational transformation R → T.
The final difference is that Eqs. (106)-(109) are defined in a quasi-one-dimensional geometry for which each mode has a wavevector with a unique z component. Consequently, there are no two modes whose wavevectors have the same z components, but different transverse components. More specifically, this means that there is no distinction between the two modes with wavevectors (κ, k z ) T and (−κ, k z ) T . If we were to enforce this constraint within our formalism, all transverse wavevectors in the sets defined in Eqs. (38) and (39) containing a minus sign would become extraneous and could be removed. Moreover, the use of the matrix σ p to correctly associate sub-matrices ofS andM that describe scattering between modes with inverse transverse wavevectors would no longer be necessary and σ p would be replaced by the identity matrix of the appropriate size. Making this change, combined with the previous change regarding the reciprocal operator, transforms Eq. (98) into Eq. (107). Finally, we see that under these changes we also have σ pp → Σ x ,M †R →M * and hence Eq. (102) becomes identical to Eq. (109), completing the comparison.

E. Numerical Example: Glass-Air Interface
In this section we demonstrate the validity of our results with a numerical example. We consider a planar glass-air boundary with glass on the left and air on the right. We choose n g = 1.5 and n a = 1.0 to be the refractive indices of the glass and air respectively. In the glass, we consider three right-travelling modes with wavevectors k 1 , k 2 and k 3 whose z components are k z1 /k 0 = 1.41, k z2 /k 0 = 0.90 and k z3 /k 0 = 0.73i, where each numerical value is given to two decimal places. In the air, we consider the set of refracted wavevectors k 1 , k 2 and k 3 , whose z components k z1 /k 0 = 0.87, k z2 /k 0 = 0.66i and k z3 /k 0 = 1.33i follow from those of k 1 , k 2 and k 3 by Snell's law. We identify the 'scattering medium' in this example with the planar boundary and note that our choice of wavevectors incorporates different types of scattering. For example, the mode with wavevector k 1 is partially transmitted and reflected, but the mode with wavevector k 2 is totally internally reflected. We describe left-travelling modes in the glass medium using the reflected wavevectors k 1 , k 2 and k 3 and, similarly, we describe left-travelling modes in the air using the wavevectors k 1 , k 2 and k 3 . For each wavevector mode, we consider both s and p polarisations (we shall use s and p to refer to polarization states in this section alone) with the usual definitions for a planar interface and write, k 1s to refer to an s polarized wave with wavevector k 1 and similarly for k 1p .
We calculate the scattering matrix element-wise using the standard Fresnel equations (see e.g. Ref. [55]). We then normalize each element of the scattering matrix according to Eq. (62), which gives where we have indicated the incident and outgoing modes above and to the right of the matrix respectively. We have included horizontal and vertical lines within the matrix as a visual aid to map out the block structure ofS as in Eq. (58). We verify numerically thatS satisfies the conservation of energy constraint given by Eq. (83) with I p = diag(I 6 , O 6 ) and I e = diag(O 6 , I 6 ). In this particular example, the scattering matrix includes evanescent components, but contains no polarization mixing and does not consider pairs of inverse trans-verse wavevector modes. As per the discussion in the previous section, the reciprocity constraint for the scattering matrix can therefore be simplified. Specifically, in Eq. (101) we may make the transformations R → T and σ pp , σ ee → I 6 , which yields the modified equation where, in this instance, ω = diag(I 6 , iI 6 ). Eq. (111) is also satisfied by our scattering matrixS.
We calculate the transfer matrix from the scattering matrix using Eqs. (69)-(72). This gives where, as before, the horizontal and vertical lines show the block structure ofM as in Eq. (61). Since there is an uneven number of propagating and evanescent modes on either side of the boundary,M has an irregular block structure and, notably, the matricesM pp andM ee are no longer square. It is therefore necessary to modify the sizes of the sub-matrices of Ω and ω in Eqs. (88) and (105) to account for this asymmetry. Furthermore, in addition to the transformations made for the scattering matrix reciprocity equation, we also replace σ pp → Σ x . Ultimately, we find thatM must satisfȳ where , Despite these modifications, we note that Eqs. (113)

IV. CONCLUSION
In this paper we have derived the general set of constraints for the scattering and transfer matrices imposed by conservation of energy, reciprocity and time reversal symmetry. Our formalism considers the general case of vectorial light and allows for fields containing evanescent components. We have extended previously known results, such as in Eqs. (106)-(109), and have introduced a more general set of constraints for the discrete scattering and transfer matrices, which are given in Eqs. (83), (88), (101) and (105). We have discussed the differences between our results and Eqs. (106)-(109) and have demonstrated that the latter can be regarded as special cases of the former. In particular, Eqs. (106)-(109) follow from our results when one is able to neglect evanescent field components, polarization mixing and a distinction between inverse transverse wavevector modes. We have given a demonstration of the validity of our results in the simple example of the scattering of polarized light at a planar glass-air interface, including the effect of total internal reflection. In our formalism, we have made minimal assumptions regarding the nature of the scattering medium and our results should therefore be valid for a wide range of systems.
Given the popularity of matrix methods in scattering studies and the increasing ease with which it is possible to experimentally determine scattering and transfer matrices, we believe that our results will serve as useful guides for future experimental research. We also believe that our results will help define the limits for physically realizable scattering and transfer matrices in matrix-based numerical simulations and theoretical studies of optical scattering. Our results should be particularly important when one is interested in a description of a scattering problem that incorporates any combination of evanescent modes, vectorial light and reciprocity/time reversal symmetry in a system of more than one dimension.