Multilayer films were analysed long before numerical solvers. Lord Rayleigh published one of the first analyses of their optical properties in 1887, and in 1917 used the traditional approach of summing multiple reflections and refractions at each interface (Joannopoulos et al., 2008, pp. 44–45). Matrix methods make that sum systematic. Berreman’s 1972 paper describes the 2 × 2 matrix technique as developed by Jones and by Abelès (Abelès, 1950), credits Teitler and Henvis with a 4 × 4 technique for anisotropic stratified media, and gives a 4 × 4 differential-matrix formulation that also covers media of low optical symmetry and liquid crystals with continuously varying order (Berreman, 1972). Yeh built a 4 × 4 transfer matrix from dynamical and propagation matrices for arbitrarily birefringent layers, including periodic birefringent stacks (Yeh, 1979).
Two later problems shaped the solvers used on this site. With evanescent waves, transfer-matrix products contain growing exponentials. Ko and Inkson presented a scattering-matrix formalism, for electron transmission in semiconductor heterostructures, that stays stable without truncating the growing states (Ko & Inkson, 1988), and Li stated the same requirement for layered gratings (Li, 1996). Computational electromagnetics borrowed transfer, scattering, hybrid and impedance matrices from transmission-line and network theory, and scattering matrices became the dominant choice (Rumpf, 2011, §1). The same layer bookkeeping was extended to dipoles inside the stack through additive source terms (Benisty et al., 1998), and to thick incoherent layers through intensity-based or coherency-matrix sums (Katsidis & Siapkas, 2002)(Postava et al., 2002). The transfer and scattering matrices page works through the matrices themselves.