Mélanie M Bay

Silvia Vignolini

Kevin Vynck
email: kevin.vynck@univ-lyon1.fr

PyLlama: a stable and versatile Python toolkit for the electromagnetic modeling of multilayered anisotropic media

Keywords: multilayers, anisotropic optical materials, optical modelling, photonic crystals, cholesterics, surface phonon polaritons

PyLlama is a handy Python toolkit to compute the electromagnetic reflection and transmission properties of arbitrary multilayered linear media, including the case of anisotropy. Relying on a 4 × 4-matrix formalism, PyLlama implements not only the transfer matrix method, that is the most popular choice in existing codes, but also the scattering matrix method, which is numerically stable in all situations (e.g., thick, highly birefringent cholesteric structures at grazing incident angles). PyLlama is also designed to suit the practical needs by allowing the user to create, edit and assemble layers or multilayered domains with great ease. In this article, we present the electromagnetic theory underlying the transfer matrix and scattering matrix methods and outline the architecture and main features of PyLlama. Finally, we validate the code by comparison with available analytical solutions and demonstrate its versatility and numerical stability by modelling cholesteric media of varying complexity. A detailed documentation and tutorial are provided in a separate user manual. Applications of PyLlama range from the design of optical components to the modelling of polaritonic effects in polar crystals, to the study of structurally coloured materials in the living world.

Introduction

Multilayered media made of anisotropic materials are widespread in nature [1] and in nanotechnologies, especially for optoelectronic [2], optical and photonic applications [3,[START_REF] Mackay | Electromagnetic anisotropy and bianisotropy: a field guide[END_REF]. For instance, chiral nematic (or cholesteric) structures consist in birefringent units arranged into periodic helicoidal architectures, which can selectively reflect circularly-polarised light in a specific wavelength range [START_REF] Oseen | The theory of liquid crystals[END_REF][START_REF] Vries | Rotatory power and other optical properties of certain liquid crystals[END_REF][START_REF] Belyakov | Optics of cholesteric liquid crystals[END_REF][START_REF] Werbowyj | Optical Properties of (Hydroxypropyl)cellulose Liquid Crystals. Cholesteric Pitch and Polymer Concentration[END_REF][START_REF] Sixou | Cholesteric mesophases[END_REF][START_REF] Frka-Petesic | The angular optical response of cellulose nanocrystal films explained by the distortion of the arrested suspension upon drying[END_REF]. In living organisms, these structures are responsible for the vibrant blue of Pollia condensata fruits [START_REF] Vignolini | Pointillist structural color in Pollia fruit[END_REF] and the shiny green of Scarabeid beetles [START_REF] Carter | Variation in the circularly polarized light reflection of lomaptera (scarabaeidae) beetles[END_REF]. They can also be fabricated through self-assembly mechanisms from bio-compatible materials, such as cellulose nanocrystals [START_REF] Parker | Hierarchical Self-Assembly of Cellulose Nanocrystals in a Confined Geometry[END_REF], cellulose derivatives [START_REF] Gray | Chiral nematic ordering of polysaccharides[END_REF] or their composites [START_REF] Saraiva | Flexible and structural coloured composite films from cellulose nanocrystals/hydroxypropyl cellulose lyotropic suspensions[END_REF][START_REF] Boott | Cellulose Nanocrystal Elastomers with Reversible Visible Color[END_REF], with applications in photonic pigments, color printing and optical sensors [START_REF] Parker | Hierarchical Self-Assembly of Cellulose Nanocrystals in a Confined Geometry[END_REF][START_REF] Bardet | Engineered pigments based on iridescent cellulose nanocrystal films[END_REF][START_REF] Yi | Ultra-Adaptable and Wearable Photonic Skin Based on a Shape-Memory, Responsive Cellulose Derivative[END_REF][START_REF] Zhao | Printing of Responsive Photonic Cellulose Nanocrystal Microfilm Arrays[END_REF][START_REF] Boott | Cellulose Nanocrystal Elastomers with Reversible Visible Color[END_REF]. In optical engineering, multilayered anisotropic media have been used to realize linear polarizers, waveplates and birefrigent filters for laser and telecommunication technologies, imaging systems, displays or gas sensing [3,[START_REF] Yeh | Optics of Liquid Crystal Displays[END_REF]. More recently, they opened new perspectives for subdiffraction wave focusing [START_REF] Dai | Subdiffractional focusing and guiding of polaritonic rays in a natural hyperbolic material[END_REF][START_REF] Li | Hyperbolic phonon-polaritons in boron nitride for near-field optical imaging and focusing[END_REF] thanks to surface phonon polaritons [START_REF] Dai | Tunable phonon polaritons in atomically thin van der waals crystals of boron nitride[END_REF] and low-loss surface wave guiding [START_REF] Pulsifer | Observation of the Dyakonov-Tamm wave[END_REF][START_REF] Mackay | Dyakonov-voigt surface waves[END_REF].

Analytical solutions for the optical properties of anisotropic multilayers are only available in specific situations (e.g., for periodic cholesterics, at normal incidence [START_REF] Vries | Rotatory power and other optical properties of certain liquid crystals[END_REF][START_REF] Chandrasekhar | Optical rotatory power of liquid crystals[END_REF] or at oblique incidence with certain restrictions [START_REF] Belyakov | Optics of cholesteric liquid crystals[END_REF][START_REF] Fergason | Cholesteric Structure-1 Optical Properties[END_REF]). A rigorous and general theoretical framework was established with the seminal contributions of Billard [START_REF] Billard | Contribution à l'étude de la propagation des ondes électromagnétiques planes dans certains milieux matériels[END_REF], Teitler and Henvis [START_REF] Teitler | Refraction in Stratified, Anisotropic Media[END_REF], and Berreman [START_REF] Berreman | Optics in Stratified and Anisotropic Media: 4x4-Matrix Formulation[END_REF], who expressed the electromagnetic problem as a 4 × 4 matrix ordinary differential equation. 4 × 4 matrix formalisms experienced many developments over the years [START_REF] Yeh | Electromagnetic propagation in birefringent layered media[END_REF][START_REF] Lin-Chung | 4 x 4 Matrix formalisms for optics in stratified anisotropic media[END_REF][START_REF] Ko | Scattering matrix method for propagation of radiation in stratified media: attenuated total reflection studies of liquid crystals[END_REF][START_REF] Eidner | Light propagation in stratified anisotropic media: orthogonality and symmetry properties of the 4x4 matrix formalisms[END_REF][START_REF] Oldano | Electromagnetic-wave propagation in anisotropic stratified media[END_REF][START_REF] Schubert | Polarization-dependent optical parameters of arbitrarily anisotropic homogeneous layered systems[END_REF][START_REF] Schubert | Explicit solutions for the optical properties of arbitrary magneto-optic materials in generalized ellipsometry[END_REF][START_REF] Stallinga | Berreman 4x4 matrix method for reflective liquid crystal displays[END_REF], including faster algorithms [START_REF] Palto | An algorithm for solving the optical problem for stratified anisotropic media[END_REF] and a correct treatment of degeneracies causing singularities [START_REF] Xu | Optical degeneracies in anisotropic layered media: Treatment of singularities in a 4x4 matrix formalism[END_REF][START_REF] Passler | Generalized 4x4 matrix formalism for light propagation in anisotropic stratified media: study of surface phonon polaritons in polar dielectric heterostructures[END_REF]. Besides commercial solutions, such as the powerful program WVASE [42] dedicated to ellipsometry analysis [START_REF] Schubert | Infrared ellipsometry on semiconductor layer structures: phonons, plasmons, and polaritons[END_REF], freely-available and open-access codes [START_REF] Castany | Berreman4x4: Python implementation of berreman's 4x4 matrix method[END_REF][START_REF] Passler | Generalized 4 x 4 Matrix algorithm for light propagation in anisotropic stratified media[END_REF] generally exploit the elegant transfer matrix method [START_REF] Mackay | The Transfer-Matrix Method in Electromagnetics and Optics[END_REF] to propagate the solution from layer to layer, as originally proposed by Berreman [START_REF] Berreman | Optics in Stratified and Anisotropic Media: 4x4-Matrix Formulation[END_REF]. Unfortunately, such an approach can become numerically unstable for large systems in presence of evanescent modes due to the coexistence of exponentially decaying and growing waves (in the forward propagating direction) in the multilayer. This shortcoming was addressed by Ko and Sambles [START_REF] Ko | Scattering matrix method for propagation of radiation in stratified media: attenuated total reflection studies of liquid crystals[END_REF], who suggested to use a scattering matrix to treat all evanescent modes as forward or backward decaying waves. Scattering matrix formalisms have gained popularity in the framework of the rigorous coupled wave analysis (RCWA) for periodically-corrugated multilayered media [START_REF] Cotter | Scattering-matrix approach to multilayer diffraction[END_REF][START_REF] Li | Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings[END_REF][START_REF] Whittaker | Scattering-matrix treatment of patterned multilayer photonic structures[END_REF][START_REF] Tikhodeev | Quasiguided modes and optical properties of photonic crystal slabs[END_REF][START_REF] Liscidini | Scattering-matrix analysis of periodically patterned multilayers with asymmetric unit cells and birefringent media[END_REF] to avoid numerical instabilities. Some freely-available RCWA codes [52,[START_REF] Liu | S4: A free electromagnetic solver for layered periodic structures[END_REF][START_REF] Germer | pySCATMECH: A Python interface to the SCATMECH library[END_REF] can be used for uncorrugated multilayered anisotropic media, yet with an unnecessary complexity. The scattering matrix method has been implemented recently for the study of cholesterics but at normal incidence only [START_REF] Kragt | 3D Helix Engineering in Chiral Photonic Materials[END_REF].

Here, we present a freely-available and open-source Python toolkit named "PyLlama"1 that implements both the transfer matrix and the scattering matrix methods to compute the reflection and transmission properties of arbitrary (linear and uncorrugated) multilayered media in all possible situations. While the transfer matrix method proves being faster, the scattering matrix method ensures numerical stability in all cases. PyLlama can thus deal with systems composed of an arbitrary number of layers of any thicknesses, mixing isotropic and anisotropic, absorbing and non-absorbing materials, studied at all possible incident angles and wavelengths for both linear and circular polarizations. Surface waves may be excited by modelling the Otto-Kretschmann configurations [START_REF] Sambles | Optical excitation of surface plasmons: an introduction[END_REF][START_REF] Passler | Generalized 4x4 matrix formalism for light propagation in anisotropic stratified media: study of surface phonon polaritons in polar dielectric heterostructures[END_REF].

In addition to robustness, we designed the code to be convenient to use and included tools that are particularly suitable for the modelling of cholesteric liquid crystals. First, layers are treated as building blocks. The optical calculations occur internally at the layer level and the user can write their own routines, in form of new classes, to pile up layers into stacks according the parameters of their choice without having to handle optical parameters. Second, the code is organised in modules. Periodic structures can be easily combined into master structures with sub-periodicities, and different pre-defined structures can be stacked onto each other, in order to enable the modelling of complex situations, such as beetle cuticles made of a combination of cholesterics and absorbing layers [START_REF] Carter | Variation in the circularly polarized light reflection of lomaptera (scarabaeidae) beetles[END_REF][START_REF] Mendoza-Galván | Pitch profile across the cuticle of the scarab beetle Cotinis mutabilis determined by analysis of Mueller matrix measurements[END_REF] or stacks of cellulose nanocrystal cholesteric structures separated by a nematic layer [START_REF] Fernandes | Mind the Microgap in Iridescent Cellulose Nanocrystal Films[END_REF]. Third, once the user has defined which structures they need to model, a level of automation has been added to allow for the calculation of spectra and the export in MATLAB and/or Pickles with only a few lines of code. Lastly, more specific to cholesterics, we incorporated a physical model [START_REF] Frka-Petesic | The angular optical response of cellulose nanocrystal films explained by the distortion of the arrested suspension upon drying[END_REF] that describes the cholesteric helicoid in various situations (choice of handedness, tilt, vertical compression of tilted helicoids). In particular, this enables the optical modelling (at normal and oblique incidences) of distorted cholesterics for which the distortion impacts the angular response, polarisation selectivity, and overall reflection spectrum of the structure [START_REF] Dreher | Reflection properties of distorted cholesteric liquid crystals[END_REF][START_REF] Frka-Petesic | The angular optical response of cellulose nanocrystal films explained by the distortion of the arrested suspension upon drying[END_REF][START_REF] Kragt | 3D Helix Engineering in Chiral Photonic Materials[END_REF][START_REF] Chan | Visual Appearance of Chiral Nematic Cellulose-Based Photonic Films: Angular and Polarization Independent Color Response with a Twist[END_REF][START_REF] Boott | Cellulose Nanocrystal Elastomers with Reversible Visible Color[END_REF].

PyLlama addresses a need from the scientific community for a simple, robust and flexible program, with many possibilities for future developments and improvements thanks to its object-oriented organisation and detailed documentation.

x y k s,→ E s,→ H s,→ k s,← E s,← H s,← x y k p,→ E p,→ H p,→ k p,← E p,← H p,← s-polarisation p-polarisation θ in θ in a b c x y Entry Multilayer stack Exit z k → k ← ... z 0 z 1 z N ψ 0 ψ 1 ψ N z N-1 Layer 0 Layer 1 Layer N-2 Layer N-1 E x -H x H y E y z 2
... The following sections are organized as follows. Section 2 describes the electromagnetic theory underlying the 4 × 4 matrix formalism, as well as the transfer matrix and scattering matrix methods to propagate the solution throughout the multilayer stack. Section 3 outlines the architecture and main features of PyLlama. Several practical examples of optical computations and comparisons with analytical solutions are given to illustrate the method and show its versatility. The performance and stability of the transfer matrix and scattering matrix methods are compared in Section 4 on the example of a cholesteric structure. Finally, we provide concluding remarks in Section 5. The Appendices contain details about the numerical analysis of the so-called partial waves in optical layers, a description of the discrete model for cholesterics, and all the material and wave parameters used to construct the figures.

ψ i E p, → E s, ← E s, → E p, ← E p, → 0 E s, → 0 i n c i d e n t

Theory

Description of the problem

We consider a multilayer stack composed of N layers perpendicular to the z axis and translationally-invariant in the (xy) plane, see Fig. 1. The layers are indexed from 0 to N -1 to match Python's indexing convention. Each layer i is bounded by two interfaces at z i and z i+1 , with an arbitrary thickness h i , and described by an arbitrary permittivity tensor ǫ i , which may contain complex values to describe absorbing materials. The multilayer stack is enclosed between two semi-infinite isotropic media called "entry" and "exit" with refractive indices n entry and n exit , respectively.

The multilayer stack is illuminated from the entry half-space by a planewave at frequency ω propagating in the (xz)-plane at an angle θ in from the surface normal z. The wavevector k of the incident planewave is thus given by

k =   k x k y k z   = k 0   K x K y K z   , (1)
with k 0 = 2π λ . Here, K x = n entry sin(θ in ) and K y = 0 are constant throughout the multilayer stack.

In the isotropic half-spaces and layers, planewaves can be decomposed into s and p polarisations, see Fig. 1(a,b). The s-polarised wave has its electric field perpendicular (s from the German senkrecht, perpendicular) to the plane of incidence (xz), along y, and the p-polarised wave has its electric field in the plane of incidence (xz). Of particular relevance for helicoidal structures, the waves can also be described in a circular polarisation basis denoted as R and L for the right and left-handed polarisations, respectively. Our objective here is to compute the complex reflection and transmission coefficients of the multilayer stack, defined as r kj and t kj , where j indicates the polarisation of the incident wave, either in the linear or the circular polarisation basis, and k indicates the polarisation of the reflected or transmitted wave, respectively. Note that for multilayers containing anisotropic layers, the cross-terms in the linear polarisation basis (r ps , r sp , t ps and t sp ) can be non-zero. The reflectance and the transmittance can then be calculated respectively as R kj = |r kj | 2 and T kj = |t kj | 2 cos(θ in )/ cos(θ out ) [START_REF] Yeh | Optical waves in layered media[END_REF], where θ out = asin [n entry sin(θ in )/n exit ] is the angle of the planewave transmitted in the exit half-space.

4 × 4 matrix formalism

Maxwell's equations set relationships between the 3 components of the electric field and the 3 components of the magnetic field of a lightwave in a given medium, which, in a matrix form, read

         0 0 0 0 -∂ ∂z ∂ ∂y 0 0 0 ∂ ∂z 0 -∂ ∂x 0 0 0 -∂ ∂y ∂ ∂x 0 0 ∂ ∂z -∂ ∂y 0 0 0 -∂ ∂z 0 ∂ ∂x 0 0 0 ∂ ∂y -∂ ∂x 0 0 0 0                  E x E y E z H x H y H z         = 1 c ∂ ∂t         D x D y D z B x B y B z         . ( 2
)
Introducing the permittivity ǫ, permeability µ and optical activity ρ and ρ ′ tensors of the material to relate the D and B fields to the E and H fields, applying translational-invariance in the x and y directions and considering harmonic fields (the exp (-iωt) convention is used hereafter), Eq. ( 2) can be rewritten as

        0 0 0 0 -∂ ∂z 0 0 0 0 ∂ ∂z 0 -iK x 0 0 0 0 iK x 0 0 ∂ ∂z 0 0 0 0 -∂ ∂z 0 iK x 0 0 0 0 -iK x 0 0 0 0                 E x E y E z H x H y H z         = - iω c         ǫ xx ǫ xy ǫ xz ρ xx ρ xy ρ xz ǫ yx ǫ yy ǫ yz ρ yx ρ yy ρ yz ǫ zx ǫ zy ǫ zz ρ zx ρ zy ρ zz ρ ′ xx ρ ′ xy ρ ′ xz µ xx µ xy µ xz ρ ′ yx ρ ′ yy ρ ′ yz µ yx µ yy µ yz ρ ′ zx ρ ′ zy ρ ′ zz µ zx µ zy µ zz                 E x E y E z H x H y H z        
.

(3) Expressing the normal (z) components as a function of the tangential (x, y) components, taking out the third and six rows and rearranging the matrices finally leads to the famous 4 × 4 matrix ordinary differential equation [START_REF] Billard | Contribution à l'étude de la propagation des ondes électromagnétiques planes dans certains milieux matériels[END_REF][START_REF] Teitler | Refraction in Stratified, Anisotropic Media[END_REF][START_REF] Berreman | Optics in Stratified and Anisotropic Media: 4x4-Matrix Formulation[END_REF] ∂Ψ(z) ∂z

= iω c ∆(z)Ψ(z), (4)
where

Ψ = [E x , H y , E y , -H x ]
T is a vector describing the field at position z. 2The elements of the matrix ∆ depend on the wavevector constant component K x and on the material permittivity, permeability and optical activity. When the material is non-magnetic and non-optically active, the matrix ∆(z)

explicitly reads ∆(z) =      -K x ǫzx ǫzz 1 -K 2 x ǫzz -K x ǫzy ǫzz Kx ǫzz ǫ xx -ǫxzǫzx ǫzz -K x ǫxz ǫzz ǫ xy -ǫxzǫzy ǫzz ǫxz ǫzz 0 0 0 1 ǫ yx -ǫyzǫzx ǫzz -K x ǫyz ǫzz ǫ yy -K 2 x -ǫyzǫzy ǫzz ǫyz ǫzz      , (5)
where the z-dependence is in the permittivity tensor ǫ.

Propagation in a homogeneous layer

Let us now consider a homogeneous layer i located between two interfaces at positions z i and z i+1 , see Fig. 1(c). The corresponding 4 × 4 matrix ∆ i is then constant throughout the layer thickness and Eq. ( 4) directly leads to

Ψ(z i+1 ) = exp (ik 0 h i ∆ i ) Ψ(z i ) = R i Ψ(z i ) (6)
where h i = z i+1 -z i is the layer thickness and R i describes the propagation from z i to z i+1 . Discarding the layer label i for the sake of legibility, there are different ways to calculate the propagator R:

1. Expanding R as an infinite sum, as exp (ik 0 h∆) = 1 + ik 0 h∆ -1 2 k 2 0 h 2 ∆ 2 +... as suggested by Berreman [START_REF] Berreman | Optics in Stratified and Anisotropic Media: 4x4-Matrix Formulation[END_REF], or via Pad's decomposition as in Ref. [START_REF] Castany | Berreman4x4: Python implementation of berreman's 4x4 matrix method[END_REF]. This is directly implemented in the matrix exponential function expm in Pythons package Scipy (version 1.4.1) and in MAT-LAB (version R2019b). 2. Factorizing ∆ = PQ ∆ P -1 where Q ∆ is a diagonal matrix containing the eigenvalues and P the eigenvectors of ∆. This leads to an eigendecomposition of R as R = PQP -1 , with identical eigenvectors and eigenvalues straightforwardly related those of ∆, see below. This was also suggested by Berreman [START_REF] Berreman | Optics in Stratified and Anisotropic Media: 4x4-Matrix Formulation[END_REF] and is implemented in Passler and Paarmann's code [START_REF] Passler | Generalized 4x4 matrix formalism for light propagation in anisotropic stratified media: study of surface phonon polaritons in polar dielectric heterostructures[END_REF][START_REF] Passler | Generalized 4 x 4 Matrix algorithm for light propagation in anisotropic stratified media[END_REF]. 3. Using an algorithm based on Sylvester's formulas applied to biaxial crystals, as proposed by Palto et al. [START_REF] Palto | An algorithm for solving the optical problem for stratified anisotropic media[END_REF], which is expected to converge faster and requires less calculation steps than Pad's decomposition.

In the second approach, which we will follow here, the fields in the layer are decomposed into four partial waves [START_REF] Yeh | Optical waves in layered media[END_REF]. The eigendecomposition has a clear physical meaning, as sketched in Fig. 2a. The matrix P contains the eigenvectors p j of ∆ as

P = p 0 p 1 p 2 p 3 , (7)
where p j = [E j,x , H j,y , E j,y , -H j,x ] T in P thus represent the partial wave field components in the homogeneous layer. The diagonal matrix Q = exp(ik 0 hQ ∆ ) depends on the corresponding eigenvalues q j of ∆ as

Q =     e ik 0 hq 0 0 0 0 0 e ik 0 hq 1 0 0 0 0 e ik 0 hq 2 0 0 0 0 e ik 0 hq 3     . (8)
It thus describes the coherent propagation of the partial waves in the layer, where the eigenvalues q j represent the z-components of the corresponding wavevectors K j,z . The four partial waves can generally be represented as a pair of waves travelling forward (towards +z, subscript →) and a pair of partial waves that travel backward (towards -z, subscript ←). The two partial waves within each pair can then be identified according to their polarisation. In PyLlama, the partial waves are sorted as follows • j = 0: forward direction (→), mostly polarised along the x axis,

• j = 1: forward direction (→), mostly polarised along the y axis,

• j = 2: backward direction (←), mostly polarised along the x axis,

• j = 3: backward direction (←), mostly polarised along the y axis but this is only one possibility 3 . Details about the numerical analysis of the partial waves is provided in Appendix A.

Let us emphasize that the sorting is quite irrelevant for obtaining the transfer and scattering matrices of the multilayer stack (when one does not wish to extract transfer and scattering matrices for subsets of the system), the only constraint being that the intermediate matrices used in the calculation be invertible. The sorting allows avoiding such situations.

z i z i+1 Layer i, h i , ε i E i+1, x H i+1, y E i+1, y -H i+1, x ψ i+1 propagation transition P i -1 Q i E i, x H i, y E i, y -H i, x ψ i → ← E i, → (z i ) E i, ← (z i ) → ← E i, → (z i+1 ) E i, ← (z i+1 ) P i z i z i+1 Layer i, h i , ε i P i, out P i, in Q i, → Q i, ← → ← ← → → ← E i, → (z i ) E i, ← (z i ) → ← E i, → (z i+1 ) E i, ← (z i+1 ) P i → ← E i+1, ← (z i+1 ) E i+1, → (z i+1 ) z y x z i , plane of the interface i k i-1 E i-1 H i-1 propagation transition E i,y E i,x E i, → (z i ) → Layer i h i , ε i k i E i H i z i+1 ,

Transfer matrix method

The computation of R i for layer i allows us to relate, via Eq. ( 6), the tangential components of the electric and magnetic fields at the interface z i to those at the interface z i+1 .

Imposing the continuity of the field tangential components at each interface naturally leads to a generalization of the propagator for a multilayer stack of N layers (N + 1 interfaces) as [START_REF] Berreman | Optics in Stratified and Anisotropic Media: 4x4-Matrix Formulation[END_REF][START_REF] Mackay | The Transfer-Matrix Method in Electromagnetics and Optics[END_REF]

Ψ N = RΨ 0 , (9)
with

R = N -1 i=0 R i = N -1 i=0 P i Q i P -1 i , (10)
and

N -1
i=0 is a left-side matrix product. 4It is usually more convenient to express the electromagnetic field in terms of s and p polarisations in the entry and exit isotropic half-spaces. For a planewave in an isotropic medium, the s-and p-polarised components of the electric and magnetic fields are related to the tangential components of the fields as [START_REF] Yeh | Optical waves in layered media[END_REF] 

   E x H y E y -H x     = L     E p,→ E s,→ E p,← E s,←     , (11)
with

L =     cos(θ) 0 cos(θ) 0 n 0 -n 0 0 1 0 1 0 n cos(θ) 0 -n cos(θ)     . (12)
Here, n = {n entry , n exit } and θ = {θ in , θ out } for the entry and exit media, respectively. Combining the corresponding matrices, L entry and L exit , at both ends of Eq. ( 9), we can now define the transfer matrix T of the multilayer stack relating the ingoing and outgoing fields at the first interface z 0 to the ingoing and outgoing fields at the last interface z N ,

    E p,→ (z N ) E s,→ (z N ) E p,← (z N ) E s,← (z N )     = T     E p,→ (z 0 ) E s,→ (z 0 ) E p,← (z 0 ) E s,← (z 0 )     . (13)
with

T = L -1 exit RL entry . (14)
The complex reflection and transmission coefficients are finally obtained as

r pp = T 30 T 23 -T 20 T 33 T 22 T 33 -T 32 T 23 r ps = T 31 T 23 -T 21 T 33 T 22 T 33 -T 32 T 23 r sp = T 20 T 32 -T 30 T 22 T 22 T 33 -T 32 T 23 r ss = T 21 T 32 -T 31 T 22 T 22 T 33 -T 32 T 23 t pp =T 00 + T 02 r pp + T 13 r sp t ps = T 01 + T 02 r ps + T 03 r ss t sp =T 10 + T 12 r pp + T 03 r sp t ss = T 11 + T 12 r ps + T 13 r ss (15)
These coefficients can conveniently be assembled into a Jones matrix and expressed on the circular polarisation basis from [START_REF] Yeh | Optical waves in layered media[END_REF]

r RR r RL r LR r LL = 1 1 i -i -1 r pp r ps r sp r ss 1 1 -i i (16) t RR t RL t LR t LL = 1 1 -i i -1 t pp t ps t sp t ss 1 1 -i i (17)

Scattering matrix method

The scattering matrix method was proposed as a solution to the numerical instabilities encountered by the transfer matrix method in certain situations, such as cholesteric liquid crystals [START_REF] Ko | Scattering matrix method for propagation of radiation in stratified media: attenuated total reflection studies of liquid crystals[END_REF]. The scattering matrix relates the waves coming in the multilayer stack from both sides to the waves coming out from it, in contrast with the transfer matrix which relates the waves at the first interface to the waves at the last interface, see Eq. [START_REF] Parker | Hierarchical Self-Assembly of Cellulose Nanocrystals in a Confined Geometry[END_REF].

Following the transfer matrix formalism above, the transition from the set of partial waves in layer i to those in layer i + 1 can simply be expressed as

P i+1 E i+1,→ (z i+1 ) E i+1,← (z i+1 ) = P i Q i E i,→ (z i ) E i,← (z i ) (18)
Knowing the travel direction and the z coordinate at which the electric fields are evaluated enables us to determine whether they are ingoing or outgoing waves, as schematised on Fig. 2b. Equation ( 18) can be rearranged with linear operations to dispatch the forward-and backward-propagating waves into ingoing and outgoing waves, leading to

P i,out E i,→ (z i+1 ) E i+1,← (z i+1 ) = P i,in E i+1,→ (z i+1 ) E i,← (z i+1 ) , E i,→ (z i+1 ) = Q i,→ E i,→ (z i ), E i,← (z i+1 ) = Q i,← E i,← (z i ), (19)
where the matrices P i,out , P i,in , Q i,→ and Q i,← are defined as

P i,out = p i,0 p i,1 -p i+1,2 -p i+1,3 , P i,in = p i+1,0 p i+1,1 -p i,2 -p i,3 , Q i,→ =     e ik 0 hq i,0 0 0 0 0 e ik 0 hq 1,0 0 0 0 0 1 0 0 0 0 1     , Q i,← =     1 0 0 0 0 1 0 0 0 0 e ik 0 hq i,2 0 0 0 0 e ik 0 hq i,3     . ( 20
) E i,→ (z i+1
) and E i,← (z i+1 ) can then be eliminated from Eq. ( 19),

P i,out Q i,→ E i,→ (z i ) E i+1,← (z i+1 ) = P i,in Q i,← E i+1,→ (z i+1 ) E i,← (z i ) . ( 21
)
This eventually leads to the definition of the scattering matrix S i,i+1 between interfaces i and i + 1, as

E i+1,→ (z i+1 ) E i,← (z i ) = S i,i+1 E i,→ (z i ) E i+1,← (z i+1 ) (22)
with

S i,i+1 = Q -1 i,← P -1 i,in P i,out Q i,→ . (23)
Evidently, the combination of scattering matrices to describe multilayer stacks cannot be as straightforward as matrix multiplications in the case of the transfer matrix method. For example, the construction of S i,i+2 would require removing E i+1,→ and E i+1,← . As above, this can be done by linear operations.

Rewriting the scattering matrices in terms of their 2 × 2 quadrants,

S i,i+1 = S (1) 00 S (1) 01 S (1) 10 S (1) 11
,

S i+1,i+2 = S (2) 00 S (2) 01 S (2) 10 S (2) 11
,

S i,i+2 = S (0) 00 S (0) 01 S (0) 10 S (0) 11 , (24)
eventually leads to the following expressions for the quadrants of S i,i+2 ,

S (0) 00 = S (2) 00 C -1 S (1
with

C = 1 0 0 1 -S (1) 01 S
(2) 10

Carrying this combination for all layers of the stack, including the isotropic entry and exit media (with analytical formulas for their eigenvectors and with their propagation matrix set to identity), and expressing the forwardand backward-propagating electric fields in terms of incident, reflected and transmitted waves eventually leads to the total scattering matrix of the multilayered stack,

    E p,→ (z N ) E s,→ (z N ) E p,← (z 0 ) E s,← (z 0 )     = S     E p,→ (z 0 ) E s,→ (z 0 ) E p,← (z N ) E s,← (z N ),     . (27)
The reflection and transmission coefficients can be then obtained straightforwardly from the scattering matrix elements, r pp = S 20 r ps = S 21 r sp = S 30 r ss = S 31 t pp = S 00 t ps = S 01 t sp = S 10 t ss = S 11 [START_REF] Billard | Contribution à l'étude de la propagation des ondes électromagnétiques planes dans certains milieux matériels[END_REF] The reflection and transmission coefficients in the circular polarisation basis can finally be obtained from Eqs. ( 16) and [START_REF] Bardet | Engineered pigments based on iridescent cellulose nanocrystal films[END_REF], respectively.

Implementation

Code architecture

PyLlama is practically structured to separate the implementation of the theory from the user, and at the same time to allow the user to access details about the optical calculations on demand. As schematised on Fig. 3, the code is organised into:

• Classes that carry the optical calculations: the class Structure represents the multilayer stack and lists the layers it consists in; the class Layer represents one layer and its partial waves, its child class Half-Space being the equivalent for the entry and exit semi-infinite media; and the class Wave represents optical waves.

• Classes that contain routines to construct Structures from appropriate parameters. They are children of a master class Model and solely serve the purpose of building Structures in a way that is easier to manipulate for the user than the class Structure itself. A few models are already implemented in PyLlama. • A class Spectrum that allows for a further level of automation for common operations such as calculating a full spectrum (in linear or circular polarisation) and exporting the calculated data. This class enables the use of the code as a black-box -though any detail on the optics remains accessible.

Optical calculations

Each partial wave in PyLlama is represented by the class Wave that contains fields for the electric field (Wave.elec), the magnetic field (Wave.magnet), the Poynting vector (Wave.poynting) as well as the constant component of the wavevector K x = q j (Wave.Kx).

The class Layer contains fields for the layer permittivity ǫ (Layer.eps), its thickness h (Layer.thickness), and the wavevector constant projection on the x axis k x = K x k 0 (Layer.Kx and Layer.k0). With these parameters, the 4 × 4 matrix ∆ (Layer.D) is calculated with the function Layer.build D() from the other fields of Layer. The class Layer also implements the calculation of the partial waves described in Section 2.3 as a list of four Wave instances stored in a field Layer.partial waves. The entry and exit isotropic media are represented by the class HalfSpace, which is a children of Layer with the only differences residing in the calculation of the eigenvectors (calculated analytically) and eigenvalues (set to 1 for no propagation phenomena).

The eigenvalues and eigenvectors are calculated from the matrix ∆ (Layer.D) with Pythons numerical package Numpy or Pythons symbolic package Sympy. The eigenvalues are unique but the eigenvectors are defined up to a normalisation factor that depends on the algorithm (for example, Numpy normalises each eigenvector so that its module is 1). The eigenvectors calculated with one or another numerical method might therefore differ by a factor -1 or i, which will inhibit their sorting in the way that is described in Appendix A. This does not prevent the use of the transfer or scattering matrix methods to go from one layer to the next. However, the eigenvectors of the entry and exit HalfSpaces need to be set in a specific order to allow for the decomposition of the incident, reflected and transmitted electric fields into their appropriate components. In general, if the eigenvectors are calculated analytically (formulas may be found in the literature for specific anisotropic media [START_REF] Yeh | Optical waves in layered media[END_REF][START_REF] Palto | An algorithm for solving the optical problem for stratified anisotropic media[END_REF][START_REF] Eidner | Light propagation in stratified anisotropic media: orthogonality and symmetry properties of the 4x4 matrix formalisms[END_REF][START_REF] Kragt | 3D Helix Engineering in Chiral Photonic Materials[END_REF][START_REF] Oldano | Electromagnetic-wave propagation in anisotropic stratified media[END_REF]), the analytical formulas enable to sort them as forward and backward pairs with distinct polarisations (as p and s in the isotropic entry and exit semi-infinite spaces) and this would allow the user to access intermediate transfer or scattering matrices for specific layers in the stack.

One key point emerging from the 4 × 4 matrix formalism is the possibility to treat all layers of the stack independently: once the matrix ∆ has been calculated for each layer (and the eigenvalues, eigenvectors and partial waves that arise from it) as instances of the Layer class, the Layers can be treated as building blocks to build a stack. In PyLlama, we represent the multilayer stack as a Structure, which contains a list of Layers (field Structure.layers) as well as the entry and exit half-spaces (fields Structure.entry and Structure.exit). Once all Layers of the Structure have been characterised, they can be combined together with the half-spaces, either with the transfer matrix method or with the scattering matrix method, to obtain the reflection and transmission coefficients.

The class Structure contains several functions that implement the optical calculations:

• Structure.build transfer matrix() calculates the transfer matrix of the stack (entry half-space, consecutive layers, and exit half-space)

• Structure. build scattering matrix to next(this layer, next layer) calculates the scattering matrix for the layer this layer and the interface between this layer and the following layer next layer

• Structure. combine scattering matrices(S ab, S bc) combines the scattering matrices S a,b S ab and S b,c S bc into the scattering matrix S a,c

• Structure.build scattering matrix() calculates the scattering matrix for the complete multilayer stack

• Structure.get fresnel(method=<"SM"|"TM">) calculate the reflection and transmission coefficients, which take the form of a Jones matrix. The user may choose the underlying matrix method with the argument method. The reflection coefficients can be converted to circular polarisation with the function Structure.fresnel to fresnel circ.

• Structure.get refl trans(method=<"SM"|"TM"|"EM">, circ=<True|False>) calculates the reflectance and transmittance directly. The argument circ enables to choose between linear and circular polarisation bases.

The transfer matrix can also be calculated with the direct exponential of ∆ (without the diagonalisation) through the choice "EM".

The user manual contains tutorials and a detailed documentation containing all implemented classes and functions as well as their parameters and returns.

Periodic stacks

Many interesting multilayer stacks are periodic (such as Bragg stacks and cholesterics). Their optical properties can of course be modelled by building a periodic structure "by hand" (e.g., a periodic list of permittivities) but they can also be characterised by their repeating unit to speed up the calculations: transfer matrices and scattering matrices can therefore be calculated for a repeating unit only, and combined together in a second step.

The basic situation is straightforward for the transfer matrix method. As shown in Eq. ( 9), the propagator R of a stack of N layers (excluding the entry and exit isotropic half-spaces) is equal to the product of the propagators of all layers, R = N -1 i=0 R i . When the N layers of the stack correspond to N per repeating units made of n layers each, the propagator of the repeating unit is R per = N -1 i=0 R i such that the transfer matrix of the N per repeating units is R = Nper-1 i=0 R per . To speed up the computation time, the decomposition of the stack is different than multiplying one matrix per repeating unit. The total number of repeating units N per is decomposed in powers of two: N per = 2 k for the appropriate list of k p s (the list of k p s associated with N per is the binary form of N per : we obtain it by converting N per to binary). We calculate the transfer matrix for sub-units that contain 2 m repeated units R m per for m between 0 and the higher power of two in the decomposition, and store them in a database. Then, we select the matrices that we need to build up the full system of N per repeated units.

For the scattering method, Eq. ( 25) shows how to combine S i,i+1 (the scattering matrix of the layer i and its next interface with the layer i + 1) with S i+1,i+2 (the scattering matrix of the layer i + 1 and its next interface with the layer i + 2) into the scattering matrix S i,i+2 , which encompasses the layer i, the layer i + 1, and the interface from the layer i + 1 to the layer i + 2. For a periodic stack, the repeating unit starts at the layer 0 and finishes at the interface between the layer N -1 layer and the layer 0 of the next repeating unit. However, the last period of the stack does not end by an interface between the layer N -1 and the layer 0, but with an interface between the layer N -1 and the isotropic half-space. Therefore, the combination of the repeated motive is only done N per -1 times (with the decomposition of N per -1 into powers of two) for the N per -1 first periods, and the last period is added afterwards.

Useful routines and customisation

In principle, a direct implementation of the theory behind the transfer matrix or the scattering matrix methods is sufficient to obtain the reflectance and transmittance for a multilayer stack once the layers have been properly constructed. However, the construction of the layers from scratch can become tedious. We designed the class Model and its children as user-friendly routines that build Structures made of appropriate Layers and HalfSpaces. Model is a parent class in which we set features that are common to all models and we wrote a few children classes, SlabModel, StackModel, StackOpti-calThicknessModel and CholestericModel. The user can therefore write a custom layer-building block of code, without having to re-write any of the underlying theory, with suitable parameters that are relevant to the physical stack being modelled.

A Model essentially represents an empty multilayer stack, consisting in two semi-infinite media separated by no layer. The idea behind the code architecture is that the children classes of Model contain additional parameters used to build a Structure that is relevant to the stack, for example:

• SlabModel represents a slab of homogeneous material and its parameters are its permittivity (optionally calculated from an inputted permittivity and a rotation angle (in radians) rotangle rad around an arbitrary axis rotaxis) eps, its thickness (in nanometers) thickness.

• StackModel represents an arbitrary multilayer stack and its parameters are a list of permittivities eps list and a list of thicknesses thickness list corresponding to its layers, as well as a number of periods N per that could in principle be set to one if the lists of permittivities and thicknesses contain an item for each layer.

• StackOpticalThicknessModel represents a stack of isotropic layers which all have the same optical thickness and its parameters are a list of refractive indices n list and the total thickness of the stack total thickness.

• CholestericModel represents a cholesteric material made of birefringent units and its parameters are a Cholesteric object representing the (potentially distorted) helicoidal architecture cholesteric, the extraordinary and ordinary refractive indices of the birefringent units n e and n o and the number of repeating units where one unit is cholesteric. Our CholestericModel class interacts with our Cholesteric library: a Cholesteric instance is a parameter to build a Cholester-icModel.

Each children class of Model contains a redefinition of the parent function Model.build structure() that overrides it and creates Layers in a routine that uses the model parameters. The user manual contains tutorials explaining how to create custom children classes appropriately.

Practical examples

Figure 4 shows examples of results obtained with our different models and validated with available analytical formulas. Figure 4a andb display the case of an isotropic Bragg stack. The system consists in a periodic unit of two layers that are repeated 10 times, including the entry and exit halfspaces that are the same as the first and last layers (with no thickness). The reflection spectra of the isotropic Bragg stack calculated with our code and the model StackModel agrees with Yeh's matrix method [START_REF] Yeh | Optical waves in layered media[END_REF] for s and p-polarisations. Here we display spectra for an angle of incidence θ in = 60 • .

Figure 4c and d display the case of an isotropic non-homogeneous slab analysed by Yeh [START_REF] Yeh | Optical waves in layered media[END_REF]. The slab has a thickness L in the order of hundreds of nanometers, a refractive index n 0 at z = 0 and a refractive index z s at z = L. The refractive index profile across the slab is continuous with an exponential profile defined by:

n(z) = n 0 n s n 0 z/L (29)
Yeh gives analytical formulas for the reflection spectrum and also shows analytically that a series of N discrete layers with appropriate thicknesses and refractive indices can be used to model the non-homogeneous slabs, with a precision that increases when the resolution N increases. We used our model discretely with 20 layers of the same thickness (orange) and discretely with 4 layers of the same thickness (green). d) Reflectance calculated with PyLlama in the discrete scenarios and with Yeh's analytical formulas in the continuous scenario [START_REF] Yeh | Optical waves in layered media[END_REF]. e, f) Uniaxial crystals extracted from a cholesteric structure with a rotation angle around the z axis φ i = 15 • and 30 • . The direction of the electric field of the four partial waves calculated numerically with PyLlama matches these calculated analytically [START_REF] Oldano | Electromagnetic-wave propagation in anisotropic stratified media[END_REF]. The angle of incidence is e) 0 • and f) 60

• .
StackModel and input a list of N permittivities (calculated from the N refractive indices) and a list of N thicknesses, all identical. We could also have used StackModel and calculated the individual thicknesses outside, it was more convenient to create a dedicated model. If we were to work often with non-homogeneous slabs, we could also have imagined a model NonHomoge-neousSlabModel taking a function n(z) (lambda z : n0 * (ns / n0) ** (z / L) in our case) as a parameter and automatically handling the splitting in layers of the same thickness. This illustrate the main idea behind the creation of models: the user can write up the classes (children of Model) that they wish in order to model structures by handling the parameters that the user think are convenient. Figure 4e shows the direction of the eigenvectors computed for a single layer of nematic crystal (taken here as an individual Layer in a Structure generated by a CholestericModel) and compares it with analytical predictions [START_REF] Oldano | Electromagnetic-wave propagation in anisotropic stratified media[END_REF]. The agreement is quantitative. It is interesting to note that the electric field follows the cholesteric rotating director at normal incidence, as predicted by De Vries [START_REF] Vries | Rotatory power and other optical properties of certain liquid crystals[END_REF] and confirmed here.

Cholesterics are of special concern in the present work. For modelling purposes, these continuous non-homogeneous materials are assimilated to discrete stacks of slices where each slice i (at depth s i on the helical axis) possesses a director n i and rotation angle φ i , leading to a permittivity tensor ǫ(φ i ). A more general discussion on cholesterics is provided in Appendix B.

The Cholesteric library allows the user to manipulate (distort, compress) the structure in a physical way, and then export the slices directors to build up the corresponding multilayer stack from which to calculate the reflectance. Details about the physical models derivation and how it handles vertical compression may be found in Ref. [START_REF] Frka-Petesic | The angular optical response of cellulose nanocrystal films explained by the distortion of the arrested suspension upon drying[END_REF]. This enables us to easily model different configurations of practical interest. Figure 5 showcases a few configurations that our Cholesteric library can generate and that we can optically model: a right-handed cholesteric, a left-handed cholesteric, oblique incidence of light upon a straight or tilted cholesteric, a vertically compressed cholesteric, and a distorted cholesteric. To account for the tilt of a helicoid, we set the z-axis along the helical axis and we calculate the effective angle of incidence upon the stack with the absolute angle of incidence and the tilt of the Cholesteric.

As we can see from the cases represented here, the handedness, the angle of incidence and the distortion strongly impact the peak wavelength, the peak shape and the polarisation selectivity (reflectance from unpolarised light can : Collection of cholesteric-specific situations implemented in PyLlama (righthanded cholesteric, left-handed cholesteric, tilt, oblique incidence, vertical compression, distortion combining a tilt and a vertical compression), and their reflection spectra spectra calculated numerically. The user can easily mix-and-match these situations. Both the 3D representations of the cholesterics and the optical spectra are automatically computed through PyLlama. even be higher than 50% in some cases). Each time, a Cholesteric object describes the physical architecture and the required parameters are extracted to build up an optical model upon it and enables to calculate reflection spectra in the circular-polarisation basis.

We use this example to illustrate results obtained with PyLlama and to demonstrate how flexible it is to handle custom structures. The user may add supplementary features to our Cholesteric library, or pair an external library to PyLlama in a similar way to model custom stacks.

Mixing multiple stacks

Throughout our code, Layers are treated as independent units and can be added to a Structure provided that the wavevector x-component stays the same. Similarly, Structures can be stacked together to model the behaviour of complex samples such as understanding the self-assembly of dry cellulose nanocrystal films that creates distorted domains [START_REF] Frka-Petesic | Retrieving the Coassembly Pathway of Composite Cellulose Nanocrystal Photonic Films from their Angular Optical Response[END_REF], explaining the polarization independence of hydroxypropyl cellulose films distorted by cross-linking [START_REF] Chan | Visual Appearance of Chiral Nematic Cellulose-Based Photonic Films: Angular and Polarization Independent Color Response with a Twist[END_REF] or tuning the distortion to obtain samples with a particular polarization selectivity [START_REF] Kragt | 3D Helix Engineering in Chiral Photonic Materials[END_REF] or gaining insight on mechanical deformations imposed on cellulose nanocrystal elastomers [START_REF] Boott | Cellulose Nanocrystal Elastomers with Reversible Visible Color[END_REF]. Combining models into a master model allows to take advantage of sub-periodicities within each model and to keep the computation time fast. The class MixedModel, child of Model, takes as additional parameter a list of models models list to be combined (in the order of the list). They are screened in order to add only the models that are compatible with the first one from the list (in terms of the conservation of the wavevectors tangential component). The function MixedModel.get refl() overrides its parent function Model.get refl() and calls for get refl TM multiple structures (respectively, SM), which builds up transfer or scattering matrices for each sub-model and combines them with a cautious handling the interfaces between consecutive models for the scattering matrix method.

An example is given on Fig. 6 where two periodic CholestericModels of opposite handedness and same pitch are stacked upon each other in a different order. The resulting spectrum for the stacked cholesterics under oblique incidence show a difference depending on which cholesteric is on top of the other.

A level of automation

The classes Structure, Layer, HalfSpace and Wave contain every function that is necessary to calculate the reflectance of a mutlilayer stack at a specific wavelength, and the class Model and its children enable to build Structures easily from appropriate parameters.

In order to enable the use of our code without requiring extensive programming skills and/or manipulation of the theoretical concepts to calculate reflection spectra for a given multilayer stack, we embedded our models inside the class Spectrum that provides a level of automation to the user to calculate the reflectance and transmittance for a range of wavelengths (to get a full spectrum) and export the data for further analysis in Python or in MATLAB, while the optics calculations are occurring in the background. The user who decides to write custom models can integrate them to the class Spectrum too.

Spectrum contains a list of wavelengths inputted by the user, a dictionary whose keys correspond to the model parameters and whose values correspond to the values assigned to these parameters inputted by the user, and an empty dictionary which will contain the calculated reflection spectra. The function Spectrum.calculate refl trans() enables to calculate reflection spectra in the linear polarisation basis or in the circular polarisation basis with a choice of the matrix method with the parameter, and with a display of the calculation progress with the parameter talk=<False|True>. Results for (potentially both) polarisation bases are then added to Spectrum.data. The function Spectrum.export() enables to export the content of Spectrum.data for further processing in Python or in MATLAB. Its argument is a filename path out which extension should be .pck for an export with Pickles and further use in Python or .mat for an export in MATLAB-compatible format.

When pairing Model and its children with the class Spectrum, spectra from complex structures are straightforward to obtain even without much experience in programming. The custom Models created by the user can also be with the class Spectrum. The user manual contains tutorials explaining how to incorporate custom children classes appropriately.

Comparison between the transfer matrix and scattering matrix methods

The calculations behind the transfer matrix and the scattering matrix methods have been detailed in Section 2.4 and 2.5, and show that the combination of subsequent Layers with the transfer matrix method simply consist in matrix multiplications, Eq. ( 9) while their combination with the scattering matrix method requires to unpack matrix elements, to combine them externally and to reassemble them into a new scattering matrix, Eq. ( 25). The scattering matrix method is consequently slower than the transfer matrix method, but it is also more robust and stable.

Figure 7a shows the evolution of selected (representative) matrix elements from the transfer and scattering matrices when the number of layers in a cholesteric structure increases (Appendix Appendix B provides details about the structure). They are plotted in the complex plane. As one can see, the values obtained with the transfer matrix method are continuously growing, while the values obtained with the scattering matrix method are converging. In fact, at each calculation step, the values of the scattering matrix represent transmission and reflection coefficients, Eq. ( 28), and their modulus is therefore always under 1, while with the transfer matrix method, reflection coefficients are calculated as ratios between given matrix elements, Eq. [START_REF] Saraiva | Flexible and structural coloured composite films from cellulose nanocrystals/hydroxypropyl cellulose lyotropic suspensions[END_REF].

Usually, reflection spectra calculated from both methods lead to the same result (Figure 7c: a cholesteric structure with a birefringence n e -n o = 0.015 and 375 pitches). However, when the number of layers becomes too large or the birefringence becomes too high, the transfer matrix method fails dramatically (Figure 7d: a cholesteric structure with n e -n o = 0.035 and 1125 pitches) while the scattering matrix remains stable. Figure 7b displays which combinations of birefringences and number of pitches are not compatible with the transfer matrix method.

The transfer matrix method remains useful as it is faster than the scattering matrix method: computing a spectrum on the visible range (400 to 800 nm with a resolution of 1 nm) corresponding to the cholesteric structure on Figure 7c takes on average 6.5s with the transfer matrix method against 14.0s with the scattering matrix method (averaged over 10 runs with HP EliteBoox x360 1030 G2, with a standard deviation under 0.1).

Concluding remarks

We have presented a stable and versatile Python toolkit for the electromagnetic modelling of anisotropic multilayer stacks. PyLlama relies on known concepts and numerical techniques, but it answers a need from the scientific community for a simple, freely-available and open-source program for the purpose. It is accessible to use without extended programming experience even for complex multilayer stacks (distorted cholesterics, stacks containing periodic sub-structures...), and it is also easy to customise without having to re-write any optical calculation.

The strength of PyLlama is undoubtedly its capability to deal with arbitrary multilayers even in extreme situations (large number of layers, high birefringence, grazing angles of incidence) thanks to the implementation of the scattering matrix method, and to switch seamlessly to the transfer matrix method for faster computation speed, when suitable. PyLlama thus allows performing a side-by-side comparison between the transfer matrix and the scattering matrix methods.

Many additional functionalities could be added in the near future without considerable efforts, like the possibility to incorporate magneto-optic effects [START_REF] Schubert | Explicit solutions for the optical properties of arbitrary magneto-optic materials in generalized ellipsometry[END_REF], and to compute spatial field distribution in the stack [START_REF] Passler | Generalized 4x4 matrix formalism for light propagation in anisotropic stratified media: study of surface phonon polaritons in polar dielectric heterostructures[END_REF] and Bloch modes for periodic media [START_REF] Cao | Stable and efficient bloch-mode computational method for one-dimensional grating waveguides[END_REF]. Classes dedicated to the acquisition of datasets other than reflection and transmission spectra could be created quite easily, such as for ellipsometry measurements [START_REF] Schubert | Infrared ellipsometry on semiconductor layer structures: phonons, plasmons, and polaritons[END_REF] or critical angle measurements [START_REF] Rivière | Determination of liquid crystal refractive indices from critical angle measurements[END_REF], which are useful tools to characterise samples. The necessary calculations are already carried out in the core of PyLlama, which paves the way to the implementation of a dedicated framework similar to our class Spectrum.

The electromagnetic modelling of space-time-modulated media -a rapidly growing topic in wave physics [START_REF] Galiffi | Photonics of time-varying media[END_REF] -would be an exciting future development of PyLlama. The transfer matrix method has recently been generalized to consider time-varying systems [START_REF] Li | Transfer matrix method for the analysis of space-time-modulated media and systems[END_REF] by writing the fields as a sum of harmonics of the modulation frequency and setting up time-varying boundary conditions. This results in matrices of larger sizes depending on the number of harmonics considered. The implementation of this idea for anisotropic materials may be of significant help to explore the multitude of exotic phenomena in such systems.

It is our hope that the scientific community will adopt PyLlama and contribute to its evolution. pair consists in a partial wave rather polarised along the x-axis and a partial wave rather polarised along the y-axis [START_REF] Yeh | Electromagnetic propagation in birefringent layered media[END_REF].

The main scope of the partial waves analysis is to identify the exponentially decaying and growing waves and sort them accordingly, which is important for the scattering matrix method. This can be done easily from the imaginary component of the eigenvalues (or wavevector) q j ≡ K j,z .

Here, we chose to rely instead on the Poynting vector along the propagation direction z, S j,z (we discard here the layer label i), defined as

S j =   S j,x S j,y S j,z   =   E j,y H j,z -E j,z H j,y E j,z H j,x -E j,x H j,z E j,x H j,y -E j,y H j,x   (A.1)
with E j,z and H j,z expressed in function of E j,x , E j,y , H j,x and H j,y ,

E j,z = - ǫ zx ǫ zz E j,x - ǫ zy ǫ zz E j,y - K j,x ǫ zz H j,y H j,z = K j,x E j,y (A.2)
The direction of the partial wave is then given by the sign of the real part of S j,z , when real, and by the sign of the imaginary part of S j,z , when complex.

For the sorting of the exponentially decaying and growing waves, we have observed that this procedure was fully equivalent to that with K j,z5 . An inspection of the two waves within each pair can reveal the partial wave polarisation. If the material is anisotropic (without the crystal axes along the (x, y) axes), each partial wave p j is analysed by calculating the ratio C(p j ):

C(p j ) = |S j,x | 2 |S j,x | 2 + |S j,y | 2 (A.3)
and comparing it to its pair. If C(p 1 ) > C(p 2 ), p 1 describes a wave rather polarised along the x axis and p 2 describes a wave rather polarised along the y axis. If the material is isotropic or anisotropic with the crystal axes aligned with the laboratory axes, this corresponds to p-and s-polarised waves, respectively. According to Ref. [START_REF] Passler | Generalized 4x4 matrix formalism for light propagation in anisotropic stratified media: study of surface phonon polaritons in polar dielectric heterostructures[END_REF], in these cases, the two partial waves within each pair can be analysed with their electric field E j instead of their Poynting vector S j :

C iso (p j ) = |E j,x | 2 |E j,x | 2 + |E j,y | 2 (A.4) Figure A.
8 illustrates the analysis and sorting of the partial waves in three particular cases (all encountered when dealing with a cholesteric structure): an isotropic layer, a uniaxial crystal (anisotropic) whose principal axes are oriented along the laboratory axes x and y, and a uniaxial crystal whose axes are not aligned with the x and y axes. First, the analysis of the direction of the Poynting vector S j,z (left column) enables to sort between forward (plotted in blue) and backward (plotted in orange) partial waves. Note that this is equivalent to sorting according to the direction of the wavevector K j,z (fifth column). Second, the direction of the electric field (third column) or of the Poynting vector (second column) along the x or y axes enables to identify the partial waves individually (plotted with different dashes). In each plot, the positive direction of the z axis is displayed in the bottom-left corner for clarity. Each segment has been normalised to the same length. The elements that are effectively used to sort the partial waves have a thick black border and the others are displayed to provide a more complete picture of the situation. The wavevector K j,z and the Poynting vector S j,z offer different information on the partial waves, indicating respectively the directions of the wave and of the energy flux. This is illustrated in ). Rivière provides analytical formulas for these critical angles [START_REF] Rivière | Determination of liquid crystal refractive indices from critical angle measurements[END_REF], reported here as vertical dashed lines. In particular, when the optical axis of the crystal lays at 45 • in the (xz)-plane, the critical angle of total reflection for p-polarisation is 49.9 • . Interestingly, at slightly smaller angles (gray shaded area on because the imaginary part of the corresponding eigenvalue is zero.

Appendix B. Discrete model for cholesterics

A cholesteric structure is an assembly of elongated birefringent units organised in a stack whose planes twist continuously around a helical axis at a periodicity called the pitch. At a given depth s along the helical axis m, the elongated units point on average towards a preferred direction n(s), which corresponds to a rotation angle φ(s) around the helical axis [START_REF] Vries | Rotatory power and other optical properties of certain liquid crystals[END_REF][START_REF] Frka-Petesic | The angular optical response of cellulose nanocrystal films explained by the distortion of the arrested suspension upon drying[END_REF], even if the direction of individual units may differ. Schematics of the birefringent units is displayed on Optically, a cholesteric structure corresponds to a permittivity tensor ǫ(φ(s)) rotating at an angle φ(s) around the helical axis [START_REF] Oseen | The theory of liquid crystals[END_REF][START_REF] Belyakov | Optics of cholesteric liquid crystals[END_REF], where ǫ(0) corresponds to a birefringent unit aligned with the plane of incidence (xz):

ǫ 0 = diag(n 2 e , n 2 o , n 2 o )
where n e and n o are the refractive indices of a birefringent unit. The optical response of the structure therefore directly relates to the periodic modulation of directors around the helical axis.

Under normal incidence of light, cholesterics reflect perfectly circularly polarised light of the handedness that matches the rotation of the helicoid (contrary to a mirror that swaps the polarisation) and in a narrow wavelength region. The central wavelength λ B of the reflection peak obeys Braggs law and is proportional to the cholesterics pitch p and its average refractive index n av such that λ B = n av p [START_REF] Vries | Rotatory power and other optical properties of certain liquid crystals[END_REF]. Under oblique incidence of light θ in , the reflection peak shifts such that λ Bragg = n av p cos(θ in ) and the structure is less polarisation-selective than at normal incidence [START_REF] Vries | Rotatory power and other optical properties of certain liquid crystals[END_REF][START_REF] Belyakov | Optics of cholesteric liquid crystals[END_REF].

Cholesteric structures that are distorted also reflect both polarisations [START_REF] Kragt | 3D Helix Engineering in Chiral Photonic Materials[END_REF][START_REF] Frka-Petesic | The angular optical response of cellulose nanocrystal films explained by the distortion of the arrested suspension upon drying[END_REF]. Recently, Frka-Petesic et al. developed a physical model to physically represent cholesteric structures, as well as the perturbation it encounters upon external constraints that affect the periodic directors and modify the optical response [START_REF] Frka-Petesic | The angular optical response of cellulose nanocrystal films explained by the distortion of the arrested suspension upon drying[END_REF].

The optics of cholesterics is not trivial and no analytical formula can be found in the literature to predict the reflectance of a distorted cholesteric, for instance.

For modelling purposes, these continuous non-homogeneous materials are assimilated to discrete stacks of slices where each slice i (at depth s i on the helical axis) possesses a director n i and rotation angle φ i , leading to a permittivity tensor ǫ(φ i ), see The partial waves are extracted from the 15th anisotropic slice and b) the 30th anisotropic slice that have a rotation angle of 15 • and 30 • around the z axis, respectively.

Figure 1 :

 1 Figure 1: a, b) Schematic of an incident and a reflected s-polarised waves (a) or p-polarised waves (b) upon the stack. The plane of incidence is the (x, z) plane. The wavevectors k i and k r form an angle θ with the z-axis. c) Schematic of the multilayer stack bounded by the entry and exit half-spaces. This schematic clarifies the notations used in the main text.

Figure 2 :

 2 Figure2: a, b) Schematic representations of the variables used to describe the light interaction with the layer i a) with the transfer matrix method and b) with the scattering matrix method. Both methods are based on the phenomena of transition and propagation, which are schematised with wide red or blue arrows labelled with the corresponding matrices used in the equations. c) 3D representation of a layer and its interfaces. Transition occurs at the interface between two layers and propagation occurs within the thickness of the layer. d) Schematic representation of the transfer matrix method and the scattering matrix method.

Figure 3 :

 3 Figure 3: Schematic of the code organisation. The optical calculations are carried out in the classes Structure, Layer with its subclass HalfSpace and Wave. The children classes of Model generate Structures from chosen parameters according to the user's need and are more convenient to handle. Some models are predefined (a slab with SlabModel, a periodic stack with StackModel, a cholesteric with CholestericModel and a combination of several models with MixedModels) and the user may write their own additional routines. Lastly, the class Spectrum give an additional level of automation and allows the user to use the code as a black box, when and if needed.

Figure 4 :

 4 Figure 4: Comparison of PyLlama's results with other methods in the case of isotropic layers and uniaxial crystals. a, b) Isotropic Bragg stack. a) Schematic of the situation:the system consists in a periodic pattern of two layers that are repeated 10 times (the entry and exit half-spaces are the same as the first and last layers). b) Reflection spectrum of the isotropic Bragg stack calculated with PyLlama and with Yeh's matrix method[START_REF] Yeh | Optical waves in layered media[END_REF] for an angle of incidence θ in = 60 • . c, d) Isotropic non-homogeneous slab with exponentially varying refractive index. c) Refractive index profile: index varying continuously (blue), discretely with 20 layers of the same thickness (orange) and discretely with 4 layers of the same thickness (green). d) Reflectance calculated with PyLlama in the discrete scenarios and with Yeh's analytical formulas in the continuous scenario[START_REF] Yeh | Optical waves in layered media[END_REF]. e, f) Uniaxial crystals extracted from a cholesteric structure with a rotation angle around the z axis φ i = 15 • and 30 • . The direction of the electric field of the four partial waves calculated numerically with PyLlama matches these calculated analytically[START_REF] Oldano | Electromagnetic-wave propagation in anisotropic stratified media[END_REF]. The angle of incidence is e) 0 • and f) 60 • .

  Figure5: Collection of cholesteric-specific situations implemented in PyLlama (righthanded cholesteric, left-handed cholesteric, tilt, oblique incidence, vertical compression, distortion combining a tilt and a vertical compression), and their reflection spectra spectra calculated numerically. The user can easily mix-and-match these situations. Both the 3D representations of the cholesterics and the optical spectra are automatically computed through PyLlama.

Figure 6 :

 6 Figure 6: a, b) Schematic representation and reflection spectra under 45 • incidence of a) a right-handed and b) a left-handed cholesteric with the same pitch (500 nm) and with the same number of pitches (30). c, d) Reflection spectra for the stacks of c) the right-handed cholesteric on top of the left-handed cholesteric and d) the left-handed cholesteric on top of the right-handed cholesteric.

Figure

  Figure Comparison between the scattering matrix (SM) and transfer matrix (TM) methods a) Evolution of two matrix of indices (0, 0) and (2, 0) in the plane. size the system increases from one layer a full pitch made of 20 layers to 125 pitches made of 125×20 layers. b) The colormap displays the combinations that cannot be calculated with the TM method. The configurations a), b) and c) are pinpointed colormap. c, d) The TM and SM methods give identical results at low and for thin systems (c) while the TM method becomes numerically unstable at high birefringence and for thick systems (d). The pitch and the Bragg wavelength λ Bragg are given in Appendix B: the pitch p is the periodicity and the Bragg wavelength equals np cos(θ in ), where n is the average refractive index of the cholesteric.

Figure A. 8 :

 8 Figure A.8: Schematic of the analysis of the partial waves (eigenvalues K j,z and eigenvectors [E j,x , H j,y , E j,y , -H j,x ]). a. An isotropic entry half-space. b. A uniaxial crystal (anisotropic) whose principal axes are oriented along the laboratory axes x and y. c. A uniaxial crystal whose principal axes are oriented at an arbitrary angle from the laboratory axes; the rotation angle of the crystal is 30 • here.

Figure A. 9 :

 9 Figure A.9: Illustration of the necessity to sort forward-and backward-propagating partial waves according to the direction of the Poynting vector. a) Schematic of the situation. b)Reflectance plots for p and s polarisations and associated critical angles of total reflection according to Ref.[START_REF] Rivière | Determination of liquid crystal refractive indices from critical angle measurements[END_REF] (vertical dashed lines). c, d) Below the critical angle for p-polarised waves (gray shaded area), the Poynting vector (describing the direction of the energy flux) and the wavevector (describing the direction of the wave front) have opposite signs.

  Fig. A.9, where we consider a slab of uniaxial crystal whose director lays in the (xz) plane (see Figure A.9a). In such situations, p-and s-polarised incident waves can have different angles of total reflection (see Figure A.9b

Figure B. 10 :

 10 Figure B.10: Schematic representation of a cholesteric structure in the (xz) plane. The cholesteric structure repeats itself every half-pitch because the directors oriented towards 0 • and 180 • are equivalent; the red and green side of the elongated units and the tips of the arrows are a guide for the eye. a) Refractive indices of a birefringent unit. b) The materials periodic twisted structure. c) A physical model consists in planes whose directors rotate around the helical axis m.

  Figure B.10a (single unit) and b (assembly of units) while Figure B.10c shows the directions of the rotating planes.

Fig

  of incidence (rad) theta in rad 0 (e) and π/3 (f) Wavelength (nm)wl nm λ Bragg = e) 716.5 f) 358.25

Python toolkit for multilayered anisotropic media

The field can indifferently be represented by the vector Ψ = [E x , H y , E y , -H x ] T as in Refs.[30,

[START_REF] Oldano | Electromagnetic-wave propagation in anisotropic stratified media[END_REF][START_REF] Palto | An algorithm for solving the optical problem for stratified anisotropic media[END_REF][START_REF] Passler | Generalized 4x4 matrix formalism for light propagation in anisotropic stratified media: study of surface phonon polaritons in polar dielectric heterostructures[END_REF][START_REF] Kragt | 3D Helix Engineering in Chiral Photonic Materials[END_REF], or by Ψ = [E x , E y , H x , H y ] T as in Ref.[START_REF] Castany | Berreman4x4: Python implementation of berreman's 4x4 matrix method[END_REF].

In Ref.[START_REF] Yeh | Electromagnetic propagation in birefringent layered media[END_REF], the partial waves are sorted according to their polarisation first and their direction second (which would correspond to the order 0, 2, 1, 3, here).

We define here the matrix Q i similarly to Berreman[START_REF] Berreman | Optics in Stratified and Anisotropic Media: 4x4-Matrix Formulation[END_REF], but its inverse (Q ′ i = Q -1i ) is sometimes used in the literature[START_REF] Passler | Generalized 4x4 matrix formalism for light propagation in anisotropic stratified media: study of surface phonon polaritons in polar dielectric heterostructures[END_REF][START_REF] Castany | Berreman4x4: Python implementation of berreman's 4x4 matrix method[END_REF]; in this case, the left-side product in Eq. (9) should be replaced by a right-side productN -1 i=0 P i Q ′ i P -1i , leading to a propagator that is the inverse of the one considered here.

The sorting with K j,z is implemented and commented in the PyLlama code for the interested user.

Acknowledgements

This work was supported by ERC grant ERC-2014-STG H2020 639088 and Philip Leverhulme Prize (PLP-2019-271) for S.V. and M. M. B. The authors are indebted to Jean-Paul Hugonin (Laboratoire Charles Fabry, Palaiseau, France) for his critical reading of the manuscript, for verifying our numerical simulations with the RCWA software RETICOLO [52], and more generaly for very insightful discussions. The authors thank Bruno Frka-Petesic (Department of Chemistry, University of Cambridge, UK) for useful discussions about the optical properties of cholesteric structures. M.M.B. thanks Sing-Teng Chua (Department of Chemistry, University of Cambridge, UK) for proof-testing the user manual, and Mathieu Brèthes for his help regarding the documentation. K.V. thanks Akhlesh Lakhtakia (Penn State University, USA) for pointing out the seminal contribution of Jean Billard to the 4 × 4 matrix formalism [28].

repository link: https://github.com/VignoliniLab/PyLlama Licensing provisions: GPLv3 Programming language: Python Supplementary material: User guide and tutorials at https://pyllama.readthedocs.io/ Nature

Appendix A. Analysis of the partial waves

The four partial waves in a layer can generally be separated into a forward (+z) propagating pair and a backward (-z) propagating pair, where each wl nm 400 to 800 For Fig. 4d, the refractive indices and thicknesses used to construct the Bragg stack (with StackModel) are calculated from Eq. ( 29) with the number of periods indicated on Fig. 4c andn The parameters used to construct the cholesteric model (with Cholester-icModel) are the same as Figure 5a.

For Fig. 5c (tilted and normal incidence), the parameters used to construct the cholesteric object (with Cholesteric) are the same as Fig. 5a except:

The parameters used to construct the cholesteric model (with Cholester-icModel) are the same as Figure 5a. For Fig. 5d (tilted and oblique incidence), the parameters used to construct the cholesteric object (with Cholesteric) are the same as Fig. 5a except

The parameters used to construct the cholesteric model (with Cholester-icModel) are the same as Fig. 5a except:

Parameter Variable Value Angle of incidence (rad) theta in rad 70π/180 For Fig. 5e (compressed), the parameters used to construct the cholesteric object (with Cholesteric) are the same as Fig. 5a except:

Parameter

Variable Value Compression ratio alpha 0.7

The parameters used to construct the cholesteric model (with Cholester-icModel) are the same as Fig. 5a. For Fig. 5f (distorted), the parameters used to construct the cholesteric object (with Cholesteric) are the same as Fig. 5a except:

Parameter

Variable Value Arbitrary distortion d 3

The parameters used to construct the cholesteric model (with Cholester-icModel) are the same as Figure 5a. The partial waves are extracted from a) the entry isotropic half-space, b) the first anisotropic slice and c) the 30th anisotropic slice that has a rotation angle of 30