Born–Kuhn coupled oscillator model for optical activity in ordered media

Abstract

Considering the phenomenological electromagnetic description of optical activity, we establish a unified formal framework, wherein the gyration and permittivity tensors of any ordered medium (a molecule or a crystal) can be derived using the classic Born–Kuhn coupled oscillator model. This model allows for efficient parameterisation of both tensors in terms of configurational and structural parameters (orientations, positions and coupling strengths) of a system of coupled oscillators having the same symmetries as the medium of interest. Thus, Born–Kuhn molecular models of all optically active crystal classes are established.

---

Introduction

Optical activity (OA), also known as circular birefringence, is most commonly defined as the ability of certain media to rotate the plane of polarisation of linearly polarised light upon transmission.^1,2 The presence of OA in a given substance is commonly considered to be due to enantiomorphism (or chirality), i.e., the existence of two non-superimposable forms of a substance that are mirror images of each other. However, detailed analysis shows that enantiomorphism is sufficient but not necessary for OA.¹ Similarly, the rotation of the plane of polarisation upon transmission, or rotatory power, is the most common manifestation of OA; however, certain media, despite being optically active, do not exhibit it.² Although experimentally discovered and phenomenologically explained in the very beginning of the nineteenth century by Biot³ and Fresnel,⁴ OA was formally described almost a hundred years later. Drude⁵ and Voigt⁶ were the first to realise that the common homogeneous constitutive relations (or material equations) complementing Maxwell's equations in matter need to be extended to account for OA, as well as to formulate such extensions. Drude likewise proposed a qualitative microscopic picture to support extended constitutive relations. Some years later, Born⁷ and Kuhn⁸ analysed Drude's qualitative picture on a quantitative basis by considering a system of coupled electronic oscillators whose spatial extent is smaller than the wavelength of the probing light. Born considered the general case of a system of coupled oscillators, whereas Kuhn focused on the most elementary structure exhibiting OA, namely, the system of two coupled oscillators in a chiral arrangement. Born was further able to derive Drude's extended constitutive relations using his microscopic model for isotropic optically active media. However, none of them applied the oscillator model (referred to as the Born–Kuhn model from now on) to ordered structures such as crystals. A first step in this yet unexplored direction was taken by Chandrasekhar,⁹ who successfully used the Born–Kuhn model to describe the OA of quartz at the molecular level. Chandrasekhar's approach was later refined and extended by V. Vyšín¹⁰ and, particularly, I. Vyšín and coworkers.¹¹

The advent of nonlinear optics and, more recently, nanoplasmonics and metamaterials has significantly revived interest in the classic Born–Kuhn oscillator model. Svirko and Zheludev¹² applied the model in its general (Born's) version to derive the gyration and permittivity tensors of an infinite optically active medium. In recent years, the simplest two-oscillator (Kuhn's) version of the model has been successfully used by several authors¹³ to phenomenologically describe the optical response of chiral nanostructures.

Developments in nanoplasmonics have been particularly relevant for resurged interest in the Born–Kuhn model. Plasmonic nanostructures, which exhibit strong interactions with light at the nanoscale, are conceptually analogous to plasmonic systems, where metallic nanoparticles act as oscillators coupled through near-field interactions.¹⁴ The model provides insights into the tunability of optical activity through geometric parameters, inter-particle coupling, and material composition, facilitating the design of plasmonic metamaterials with tailored chiroptical properties for applications in sensing, polarization control, and photonic devices. However, very few studies¹⁵ have explored scenarios beyond the two-oscillator model, particularly in complex three-dimensional arrangements that essentially reproduce crystal structures. These configurations are increasingly relevant owing to advancements in modern design and nanofabrication tools, which allow for precise engineering of intricate geometries. Investigating such multi-oscillator systems could unveil new regimes of optical activity, enabling unprecedented control over light–matter interactions and expanding the functional capabilities of photonic materials.

The purpose of the present paper is twofold. Firstly, it revisits the phenomenological description of OA, creating a formal unified framework within which it derives the gyration and permittivity tensors from the Born–Kuhn coupled oscillator model for finite and infinite ordered media. Secondly, it establishes, as an illustrative application, parametric forms of gyration tensors for all optically active crystal classes. The framework bridges the gap between microscopic descriptions of optical activity based on coupled oscillator models and macroscopic electromagnetic theory of matter based on the gyration and permittivity tensors, both of which can be experimentally probed through polarised light measurements. Since the approach is general, it allows for efficient geometric representation and analytic parameterisation of any natural or artificial optically active structure as a function of its symmetries.

Constitutive relations for optically active media

Herein, we assume plane monochromatic wave dependence, exp(ikr − iωt), for all electromagnetic vectors (k is the wavevector, r is the position and ω is the cyclic frequency of the wave). The usual linear homogeneous constitutive relations for anisotropic nonmagnetic media:¹⁶

D = ε₀εE (1a)

B = μ₀H (1b)

relate the electric (magnetic) induction (D (B)) to the electric (magnetic) field (E (H)) through the vacuum permittivity (permeability) (ε₀ (μ₀)) and the permittivity tensor (ε) of a generally anisotropic medium. From an optical viewpoint, these relations describe only linear polarisation effects such as linear birefringence and dichroism, i.e., the difference in phase velocities and extinction coefficients for polarised light waves with orthogonal linear (e.g. horizontal and vertical) polarisations. To formally describe circular polarisation effects, namely, circular birefringence (or OA) and dichroism, originating from different phase velocities and extinction coefficients of circularly polarised light waves, respectively, with orthogonal (left-handed and right-handed) polarisations, a more general, inhomogeneous form of linear constitutive relations must be used:

D = ε₀εE + P_OA (2a)

B = μ₀(H + M_OA), (2b)

where P_OA and M_OA are additional polarisation and magnetization terms arising from OA, respectively. As shown phenomenologically by Drude⁵ and Voigt,⁶ as well as formally by Boys¹⁷ and Condon,¹⁸ the OA polarisation and magnetization terms entering the electric and magnetic constitutive relation, respectively, depend on quantities of an opposite physical (magnetic or electric) nature, i.e.P_OA = P_OA(B) and M_OA = M_OA(E).

By performing Fedorov's transformation,¹⁹

D′ = D − ω⁻¹k × M_OA (3a)

H′ = H + M_OA (3b)

leaving invariant Ampère's law,

k × H = −ωD (4a)

and not affecting Faraday's law,

k × E = ωB (4b)

the general magnetic constitutive relation (2b) can be made homogeneous, such as eqn (1b), but at the cost of redefining electric induction from D to D′ and magnetic field from H to H′. Indeed, the elimination of D and H from eqn (2) with the help of eqn (3) transforms the former set of constitutive relations as follows.

D′ = ε₀εE + P_OA − ω⁻¹k × M_OA (5a)

B = μ₀H′ (5b)

Contrarily, by expanding the permittivity tensor in powers of the wavevector within the framework of the spatial dispersion picture, Landau and Lifshitz¹⁶ obtained the following set of constitutive relations for optically active media:

D′ = ε₀εE + iε₀cω⁻¹(gk) × E (6a)

B = μ₀H′ (6b)

where g is the gyration tensor describing OA and c is the speed of light in vacuum. Comparison of the electric constitutive relations from the two sets (5) and (6) gives the following relation.

P_OA − ω⁻¹k × M_OA = iε₀cω⁻¹(gk) × E (7)

Eqn (7), first established by Pauli²⁰ in the special case of isotropic media, allows the determination of the gyration tensor when OA-induced polarisation and magnetization of the medium are known. As given in the following section, the latter can be derived from the microscopic description of the medium. Thus, Pauli's eqn (7) provides the link between microscopic and macroscopic cases on the formal basis of the Born–Kuhn approach.

Before ending this section, it should be emphasized that the constitutive relations (6) featuring a homogeneous magnetic relation are valid for infinite media only, as noted by Landau and Lifshitz.¹⁶ When dealing with layered optically active media instead (as is commonly the case in optical characterisation), the following Fedorov's transformation must be performed,

D′ = D + iε₀cω⁻¹k × f^TE (8a)

H′ = H − iε₀cf^TE (8b)

bringing eqn (6) into the inhomogeneous set,

D = ε₀εE +iε₀cfB (9a)

B = μ₀(H − iε₀cf^TE) (9b)

where f is the OA tensor and is related to the gyration tensor (g) as follows:²¹

f = ½trgI − g^T, (9c)

(the notation ‘tr’ stands for the trace, whereas I denotes the identity matrix). Set (9) was established first by Drude⁵ and Boys¹⁷ in the case of isotropic optically active media. It satisfies the common boundary conditions at interfaces (continuity of the tangential components of the fields and the normal components of the inductions), as demonstrated theoretically²² and experimentally.²³ Note that set (9) is similar to generic set (2), its OA-induced polarisation and magnetization terms being dependent on quantities of an opposite physical nature, i.e.P_OA = P_OA(B) and M_OA = M_OA(E).

Coupled oscillator model for optically active media

The Born–Kuhn model for OA, largely inspired by the classic Lorentz microscopic model of electromagnetism,²⁴ considers a system of electrons constrained to oscillate about spatially-separated equilibrium positions under the action of an externally applied electric field. The spatial extent (a) of the system of electrons is sufficiently smaller with respect to the wavelength (λ) of the electric field; therefore, only linear expansions in the a/λ ratio can be considered.

The approach is as follows: First, the “mechanical” part of the model, the oscillation amplitudes of all electrons, is determined from the equations of motion. Electrons are assumed to oscillate along straight lines under the influence of the time-varying electric field of probing light. The second “electromagnetic” part of the model consists of the evaluation of the polarisation and magnetization of the medium to various multipole degrees of approximation, in either finite medium (molecule) or infinite medium (crystal) cases. For simplicity, the induced electric and magnetic multipole densities constituting polarisation and magnetization are evaluated classically; although a quantum-mechanical approach exists,²⁵ it has been shown to yield similar results.²⁶ Finally, eqn (7) is applied to retrieve the gyration tensor from OA-induced contributions to polarisation and magnetization, while the permittivity tensor is obtained with the help of the electric constitutive relation (5a).

Equations of the motion of a system of coupled oscillators

For simplicity and following Kuhn's original derivation,⁸ we shall start with a system of only two coupled oscillators, whose results are readily generalizable to an arbitrary number. Consider two electrons oscillating along straight lines with respective unit vectors u₁ and u₂. The equilibrium positions r₁ and r₂ of the electrons are measured with respect to the centre (r₀) of the system. The time evolutions of the displacements a₁(t) and a₂(t) of the electrons with respect to r₁ and r₂ are given by the following simultaneous coupled equations of motion:where k₁₂ = k₂₁ = k is the coupling between two electrons with charge e and mass m_e. Physically, the coupling originates from the multipole–multipole (typically, dipole–dipole) electromagnetic interaction between the oscillating electrons.²⁷ In eqn (10), E_n(t) = E(t, r₀ + r_n), n = 1, 2, is the probing light electric field acting upon the nth oscillator at the equilibrium position; its projection upon the oscillation direction (u_n), multiplied by the charge, equals the driving force. Finally, ω₁ and ω₂ are the frequencies of free oscillations of the two electrons considered uncoupled, i.e. if E_n(t) = 0 and k = 0.

For a monochromatic electric field (∝exp(−ωt)), eqn (10) can be represented in the matrix form,

which is readily solvable for the two oscillation amplitudes a₁ and a₂,where

D(ω) = (ω₁² − ω²)(ω₂² − ω²) − k², (12b)

Eqn (11) and (12) express the linear relationships existing between the amplitudes of the two oscillating electrons and the magnitude of the electric field driving them out of their respective equilibrium positions. The generalization of these equations in the case of N > 2 oscillators whose oscillation amplitudes are grouped in the vector a is immediate. Eqn (11) then takes the following form:

Da = e′, (13)

where D is the N × N “dynamic” matrix with elements D_nm = D_mn = m_ek_mn/e for m ≠ n and D_nn = m_e(ω_n² − ω²)/e, whereas e′ is the “electric” vector with components , n = 1, 2, …, N. Note that the matrix D is symmetric by virtue of the reciprocity of the coupling within the pairs of oscillators.

The solution of eqn (13) is obviously:

a = D⁻¹e′ (14)

or, in component form (summation over repeated indices being understood, unless otherwise stated),

(Note that D_mn⁻¹ designates the mnth element of the inverse matrix D⁻¹ rather than the reciprocal of D_mn.)

Born–Kuhn model for a finite structure (molecule)

As mentioned in the previous section, the second step in the Born–Kuhn approach consists of evaluating the polarisation and magnetization of the optically active structure (molecule or crystal) by calculating the microscopic electric and magnetic multipole moments to different orders of approximation using the oscillation amplitudes of the system of electronic oscillators. Raab and De Lange²⁸ showed that the first, electric dipole order of approximation is insufficient to account for OA; for a complete description of the phenomenon, calculations should be performed to the second-order electric quadrupole–magnetic dipole. In our calculations, we shall consider both orders consecutively, the first being an essential intermediate step for the second, and critically compare the results obtained.

Electric dipole approximation. Within this first-order approximation of multipole theory, the polarisation P of the molecule is identified with the induced electric dipole moment density p of the system of electronic oscillators,where v ∼ a³ is the volume occupied by a molecule of size a. Substituting eqn (15) for the oscillation amplitudes a_m yields the polarisation.

We have put A_mn = ε₀⁻¹eD_mn⁻¹/v and have expressed the electric field E_n applied on the nth oscillator through the field E = E(0) at the molecular centre r₀ = 0 by means of linear expansion, i.e., E_n = E(r_n) ≈ (1 + ikr_n)E.

Again, within the first-order approximation, the magnetization M, identified by the induced magnetic dipole moment density m of the system of oscillators, is neglected since it is of a higher order with respect to the electric dipole moment density, i.e., M = m = 0. Next, we evaluate the expression P − ω⁻¹k × M = p with the help of eqn (17).

The above expression contains both OA- and non-OA contributions. Comparison of eqn (18) and (7) allows to isolate the OA contribution as the last term:

P_OA − ω⁻¹k × M_OA = i[½ε₀A_mn(u_n × u_m)r^T_n]k × E (19)

and consequently, to identify the gyration tensor of the molecule:

g′ = ½c⁻¹ωA_mn(u_n × u_m)r^T_n (20)

(We have used a primed notation for the tensor to emphasize the first order of the approximation.) To give meaning to the first non-OA term in eqn (18), the redefined electric constitutive relation (5a) can be written in the general inhomogeneous form:

D′ = ε₀E + P′ (21)

where polarisation P′ is redefined as follows.

P′ = ε₀(ε′ − I)E + P_OA − ω⁻¹k × M_OA (22)

The last two terms in eqn (22) represent the OA-contribution already identified in eqn (18). Consequently, the first term corresponds to the non-OA contribution; when put in correspondence with the first term from eqn (18), it provides the permittivity tensor of the molecule.

ε′ = I + A_mnu_mu^T_n + i½A_mnkr_n(u_mu^T_n + u_nu^T_m) (23)

This result is unphysical since the last term explicitly depends on the wavevector. Furthermore, the presence of the last term breaks the essential symmetry property of the permittivity tensor, ε^T = ε, following from the conservation of energy.^16,20 In fact, the Lorentz classic electron theory of the permittivity in the absence of OA yields only the first two physical terms in eqn (23).¹² This points out the deficiency of the first-order electric dipole approximation. To obtain the correct forms of permittivity and gyration tensors, second-order electric quadrupole–magnetic dipole approximation must be considered.

Electric quadrupole–magnetic dipole approximation. In this order of multipole approximation, polarisation (P) of the molecule comprises the induced electric dipole moment density (p) and the induced electric quadrupole moment density (q) contributions.^28,29

P = p − i½kq (24)

Dipole moment density (p) is given by eqn (16), whereas the induced electric quadrupole moment density of the system of electronic oscillators is:

at the first order in the oscillation amplitude (a_m), substituted for from eqn (15). At this order of approximation, magnetization (M) of the molecule is given by the induced magnetic dipole moment density (m) of the system of oscillators,into which expression (15) for a_m has been substituted, as previously mentioned. Note that magnetization depends on the electric field and not on the magnetic field, i.e., M = M(E). This observation implies that the magnetization of a molecule is entirely due to OA, M = M_OA, in accordance with the discussion following the constitutive relations (5). Generally, the Born–Kuhn model based on electrons oscillating along straight lines can model only non-magnetic media, i.e. it formally yields μ = I for the permeability tensor.

The expression from eqn (18) involving polarisation and magnetization is therefore as follows.

P − ω⁻¹k × M = p − i½kq − ω⁻¹k × m (27)

The first term on the right-hand side of eqn (27) is the electric dipole contribution already calculated in eqn (17). The quadrupole contribution due to the last two terms is evaluated using eqn (25) and (26).

−i½kq − ω⁻¹k × m = −iε₀A_mnkr_mu_mu^T_nE (28)

The addition of eqn (17) and (28) finally yields the right-hand side of eqn (27).

As mentioned in the previous section, identification of the above expression with eqn (19) and (22) yields:

g = c⁻¹ωA_mn(u_n × u_m)r^T_n (30)

for the gyration tensor and

ε = I + A_mnu_mu^T_n (31)

for the permittivity tensor of the molecule, respectively.

Comparison of the two expressions for the gyration tensor obtained at different degrees of multipole approximations, as given by eqn (20) and (30), shows that the electric quadrupole–magnetic dipole contribution is equal in magnitude and of the same structure as the electric dipole contribution alone since g = 2g′. This observation confirms the general result from multipole theory²⁸ that taking the multipole order next to the electric dipole order into account is essential for the complete quantitative description of the OA phenomenon. Furthermore, addition of the electric quadrupole–magnetic dipole contribution effectively removes the unphysical dependence of the permittivity tensor on the directional derivatives (compare ε′ from eqn (23) with ε from eqn (31)).

It should be noted that the electric quadrupole contribution (q) can be neglected with respect to the magnetic dipole contribution (m) if an ensemble of randomly oriented molecules, representing an isotropic optically active substance, is considered. It can be easily observed from eqn (28) that if q = 0, then the correct expression (30) for the gyration tensor is still obtained; however, the expression for the permittivity tensor contains extra unphysical directional derivative terms, such as those in eqn (23). These spurious terms vanish upon the orientational averaging of the molecule, as shown by Caldwell and Eyring,²⁷ as well as in the original works of Born⁷ and Kuhn.⁸ Nevertheless, as already discussed, the electric quadrupole contribution is essential for obtaining the correct results for both tensors in case of an oriented (i.e., an orientationally non-averaged) molecule. In ordered systems, where the molecular orientation is not random, the electric quadrupole–magnetic dipole contribution ensures that the model effectively predicts the optical response of the material in practice, irrespective of the propagation direction of the probing light.

Born–Kuhn model for an infinite structure (crystal)

The case of an infinite structure is computationally simpler since it is sufficient to work to the first-order, electric dipole, of multipole theory. In fact, as shown by Robinson²⁹ and Rosenfeld,³⁰ the consideration of the second-order, electric quadrupole–magnetic dipole, is necessary only when the structure is inhomogeneous or, in particular, discontinuous, such as in the case of finite (molecular) structures.

Following Svirko and Zheludev,¹² we modelled the crystal as a three-dimensional regular lattice of N unit cells (N ≫ 1) of volume v. Thus, the induced electric dipole moment density of the unit cell whose centre is located at r₀ is as follows.

p(r₀) = eu_ma_mδ(r − r₀ − r_m) = eD_mn⁻¹u_mu^T_nE_nδ(r − r₀ − r_m) (32)

In accordance with the electric dipole approximation, polarisation (P) of the crystal equals the induced electric dipole moment density averaged over the lattice:

where E(r + r_n − r_m) has been replaced with E_n−m to simplify the notations. Since M = 0 at this order of approximation, the expression P − ω⁻¹k × M appearing in eqn (18) is transformed to the following one:

In eqn (34), electric field (E_n−m) has been linearly expanded about its value E = E(r) at an arbitrary position (r) in the crystal, i.e.E_n−m = E(r + r_n − r_m) ≈ (1 + ikr_n − ikr_m)E.

Identification of eqn (34) with eqn (7) and (22) yields the gyration tensor,

g = c⁻¹ωA_mn(u_n × u_m)r^T_n (35)

as well as the permittivity tensor,

ε = I + A_mnu_mu^T_n (36)

of the crystal.

The Born–Kuhn model provides identical results for finite and infinite optically active media. Indeed, the expressions for the gyration tensors, as well as those for the permittivity tensors, coincide; compare eqn (30) and (31) with eqn (35) and (36), respectively. In the framework of the Born–Kuhn model, the properties of the crystal are entirely determined by the properties of its unit cell, i.e. the crystal is equivalent to an oriented molecule.

Another prominent feature of the Born–Kuhn model is that it provides not only the gyration tensor of the optically active medium but also the permittivity tensor, both tensors being closely related physically. An important practical consequence of this fact is that the same model should be used to describe both tensors when fitting the simulated optical response of the medium or structure under study to experimental data.

Eventually, it should be noted that both gyration and permittivity tensors depend on the cyclic frequency (ω) of the probing light either explicitly or implicitly through the matrix element A_mn. This phenomenon, called time or frequency dispersion,¹⁶ which results in the spectral dependence of both tensors, has been addressed in detail by many authors over the years^9–11 and is still an active topic;³¹ thus, it will not be dealt with here. Indeed, spatial and time (or frequency) dispersions are treated independently in theoretical models considering their distinct physical origins.¹⁶ Alternatively, spatial dispersion results from the non-local response of the medium, while time dispersion arises from the finite-time response of the medium.

Properties of the gyration tensor

The gyration tensor of the molecule (or the crystal) given by eqn (30) (or by eqn (35)) can be readily represented as

g = c_mn(u_n × u_m)v^T_n (37)

by substituting r_n = d_nv_n and c_mn = d_nc⁻¹ωA_mn (no summation over repeated indices being assumed), where d_n is the distance of the nth oscillator from the centre of the molecule (or of the unit cell) and v_n is the unit vector along r_n. It should be noted that while A_mn satisfies A_mn = A_nm, the matrix formed by the coupling coefficients c_mn is not symmetric, i.e., c_mn ≠ c_nm, because d_n ≠ d_m, in general. Note that eqn (30), (35) and (37) present the gyration tensor in the form of a dyadic, i.e., a linear combination of dyads (or outer products of vectors).

In the case of only two coupled oscillators i and j forming Kuhn's pair, eqn (37) reduces to:

g_ij = (u_i × u_j)(c_jiv^T_i − c_ijv^T_j) (38)

(no summation over repeated indices), where u_i and u_j are the oscillation directions of the two oscillators forming the pair, while v_i and v_j are the unit vectors pointing at their positions.

In isotropic disordered optically active media, the gyration tensor must be averaged over all possible orientations. Tensor (37) then becomes a scalar equal to its trace, i.e., g = trgI, where

trg = c_mn(u_n × u_m)v_n = c_mn(u_nu_mv_n) (39)

in which (u_nu_mv_n) = (u_n × u_m)v_n is the scalar triple product of vectors.

Before proceeding with the identification of the properties of gyration tensor g, it should be recalled that strictly speaking, the second-rank tensor g is not a true (or a polar) tensor but rather a pseudo-tensor (or an axial tensor), as shown by Landau and Lifshitz.¹⁶ In fact, when written in component form in the general (non-plane-wave) case, the redefined electric constitutive relation (6a) becomes:

where

γ_ikn = cω⁻¹e_ijkg_jn (41)

is the true (or polar) third-rank gyration tensor obeying the anti-symmetry relation γ_ikn = −γ_kin by virtue of Onsager's principle (e_ijk is the Levi–Civita symbol, i.e., the fully antisymmetric third-rank pseudo-tensor). Note that eqn (40) is the first-order expansion of the redefined electric induction D′ in the position and ∇ is the gradient operator, in agreement with the spatial dispersion picture of OA.¹⁶Eqn (41) defines g as a second-rank axial tensor (or pseudo-tensor), unlike the third-rank polar (or true) tensor (γ). Consequently, gyration tensor g transforms similar to a pseudo-tensor––it behaves similar to a polar tensor under rotation transformations; however, unlike a polar tensor, it changes its sign upon reflection transformations.¹⁶

We can now draw a list of properties satisfied by the gyration tensor that follow directly from eqn (37).

1. If u_n → −u_n, then g → g. Time reversal, i.e., inversion of the oscillation directions of all oscillators, does not alter the tensor.

2. If v_n → −v_n, then g → −g. Space inversion, i.e., inversion of the system of oscillators to its centrosymmetric image, reverses the sign of the tensor, in agreement with g being a pseudo-tensor. Therefore, if the system itself is centrosymmetric, then g = 0 (centrosymmetric media do not exhibit OA, as is well known).¹⁶

3. If u_n → −u_n and v_n → −v_n, then g → −g. Simultaneous time reversal and space inversion reverse the sign of the tensor.

4. If R is a rotation transformation leaving the system of oscillators invariant, i.e., if Ru_m = u_n and Rv_m = v_n, then RgR^T = g, i.e., the tensor is unaltered. This property is in agreement with g being a pseudo-tensor behaving as a true (polar) tensor with respect to rotations.

5. If S is an inversion (e.g., mirror reflection) operation leaving the system of oscillators invariant, i.e., if Su_m = u_n and Sv_m = v_n, then SgS^T = −g, i.e., the tensor reverses its sign. Similar to the previous property, this is in agreement with g being a pseudo-tensor and changing its sign upon inversion.

6. If the system of oscillators decomposes into two subsystems that are uncoupled from one another with respective tensors g₁ and g₂, then g = g₁ + g₂. This is an obvious consequence from g being a dyadic, i.e., a sum of dyads. In particular, if every oscillator from the system is coupled to another one only, i.e., if the coupling is pairwise, then gyration tensor g of the whole system is made up of Kuhn's pair contributions of the form given by eqn (38).

7. The gyration tensor vanishes when

(a) c_mn = 0 for m ≠ n, i.e., all oscillators are uncoupled. Coupling is essential for OA.

(b) v_n = 0: a system of centred oscillators, i.e., oscillators whose equilibrium positions coincide with the centre of the molecule (or of the unit cell) exhibit no OA. Spatial dispersion is essential for OA.¹⁶

8. The trace of the gyration tensor given by eqn (39) vanishes when

(a) u_m, u_n and c_mnv_n − c_nmv_m are coplanar, i.e., all oscillators lie in a plane.

(b) u_m||v_m, i.e., all oscillators are radial, i.e., oscillations take place along the radius vectors joining the oscillators to the centre of the molecule (the unit cell). Note that this condition is equivalent to pairwise coplanar oscillators, i.e., coplanar Kuhn's pairs. Since the gyration tensor of each coplanar Kuhn's pair is traceless, thus, the tensor of the whole system will also be traceless (being a sum of Kuhn's pair contributions by virtue of property 6.)

Consequently, a system of coplanar or radial (i.e., pairwise coplanar) oscillators exhibits no molecular enantiomorphism.

Coupled oscillator models for optically active crystal classes

The formal coincidence of the expressions for the gyration tensor of a molecule and a crystal, given by eqn (30) and eqn (35), respectively, suggests the possibility of constructing molecular models for all optically active crystal classes. The molecular model of a given crystal class consists of a system of coupled oscillators whose spatial configuration obeys the symmetries of the class. While constructing the model, the total number of oscillators and the coupling between them can be minimised. As a result, all optically active crystal classes can be modelled by systems of coupled oscillators containing a certain number of Kuhn's pairs. Consequently, we shall establish the molecular models of all optically active crystal classes. The configurations are shown in figures in which oscillators are represented by red rods. Coupled oscillators are connected by straight lines. Equal symmetric coupling between oscillators within different pairs (i.e. A_ij = A_mn, c_ij = c_ji and c_mn = c_nm for pairs ij and mn) is represented by connecting lines of the same colour, while distinct colours indicate unequal (A_mn ≠ A_nm) or asymmetric (c_ij ≠ c_ji, c_mn ≠ c_nm) coupling.

Biaxial crystals

Class 1. This triclinic class exhibits no symmetries, and therefore, all nine tensor elements are generally non-zero. To model it, three arbitrarily oriented oscillators are considered:

u_k = [sin θ_kcos φ_k sin θ_ksin φ_k cos θ_k]^T, (42)

which are located along the coordinate axes v₁ = [1 0 0]^T, v₂ = [0 1 0]^T and v₃ = [0 0 1]^T, respectively, as shown in Fig. 1. In eqn (42), φ_k and θ_k are the azimuth and the polar angle of the kth oscillator, respectively (k = 1, 2, and 3). The gyration tensor results from the pairwise coupling of the triplet of oscillators:

g₁₂ + g₂₃ + g₃₁ (43)

and is obtained using eqn (38) for g_ij. The model has twelve parameters, nine of which can be considered independent, while the remaining three can be assigned predefined values to match the number of non-zero tensor elements. The nine independent model parameters uniquely determining the non-zero tensor elements are the three polar angles θ₁, θ₂ and θ₃ together with the six coupling coefficients c₁₂, c₂₁, c₂₃, c₃₂, c₁₃ and c₃₁, whereas the remaining three azimuths φ₁, φ₂ and φ₃ may be ascribed (non-uniquely) to φ₁ = 0, φ₂ = π/4 and φ₃ = π/2.

Fig. 1 Molecular model for crystal class 1.

Class m. This monoclinic class features a single mirror plane coinciding (by convention) with the xy-plane. Its gyration tensor is traceless, does not exhibit any symmetries, and has four non-zero elements.

To model it, an arbitrarily oriented oscillator is considered:

u₁ = [sin θ₁cos φ₁ sin θ₁sin φ₁ cos θ₁]^T (45a)

coupled with a second one that is its mirror image with respect to the xy-plane.

u₂ = [sin θ₁cos φ₁ sin θ₁sin φ₁ −cos θ₁]^T (45b)

Both oscillators are located along the z-axis symmetrically with respect to the origin so that v₁ = −v₂ = [0 0 1]^T and d₁ = d₂. A third oscillator lies entirely in the xy-plane, i.e.,

u₃ = [cos φ₂ sin φ₂ 0]^T (46a)

and

v₃ = [cos φ₃ sin φ₃ 0]^T (46b)

and is coupled with the first two. This configuration is shown in Fig. 2. The tensor is given by eqn (43) and (38), with c₁₂ = c₂₁ = c₁, c₁₃ = c₂₃ = c₂ and c₃₂ = c₃₁ = c₃ (because d₁ = d₂ and A₁₃ = A₂₃). The model has seven parameters, four of which can be considered as independent. The four independent parameters are the azimuth φ₃ and the three coupling coefficients c₁, c₂ and c₃, whereas the remaining three parameters may be fixed to φ₁ = θ₁ = π/4 and φ₂ = −φ₁.

Fig. 2 Molecular model for crystal class m.

Class 2. This monoclinic class has a single 2-fold axis taken (conventionally) along the z-axis. Its gyration tensor does not exhibit any symmetry and has five non-zero elements.

It can be modelled with two Kuhn's pairs. The first pair consists of two arbitrarily oriented oscillators:

u₁ = [sin θ₁cos φ₁ sin θ₁sin φ₁ cos θ₁]^T (48a)

u₂ = [sin θ₂cos φ₂ sin θ₂sin φ₂ cos θ₂]^T (48b)

whose position directions are considered symmetric with respect to the x-axis.

v₁ = [sin θ_pcos φ_p sin θ_psin φ_p cos θ_p]^T (49a)

v₂ = [sin θ_pcos φ_p −sin θ_psin φ_p −cos θ_p]^T (49b)

The pair is rotated through π about the z-axis to produce the opposite, z-axis-symmetric second pair,

u₃ = [−sin θ₁cos φ₁ −sin θ₁sin φ₁ cos θ₁]^T (50a)

u₄ = [−sin θ₂cos φ₂ −sin θ₂sin φ₂ cos θ₂]^T (50b)

and

v₃ = [−sin θ_pcos φ_p −sin θ_psin φ_p cos θ_p]^T (51a)

v₄ = [−sin θ_pcos φ_p sin θ_psin φ_p −cos θ_p]^T (51b)

The oscillators within each pair and the first and second oscillators from the opposite pairs are coupled so that all oscillators are coupled “peripherally”; see Fig. 3. The tensor is given by the sum of the pairwise contributions,

g₁₂ + g₃₄ + g₂₃ + g₄₁ (52)

and eqn (38) with c₁₂ = c₃₄ = c₁, c₂₁ = c₄₃ = c₂, c₂₃ = c₄₁ = c₃ and c₃₂ = c₁₄ = c₄ (because d₁ = d₃ and d₂ = d₄, as well as A₁₂ = A₃₄ and A₂₃ = A₁₄). The model has ten parameters, five of which can be fixed. The independent model parameters are the polar angle θ₁ and the four coupling coefficients c₁, c₂, c₃ and c₄, whereas the fixed parameters are φ₁ = π/4, φ₂ = −φ₁, θ₂ = −θ₁, φ_p = π/4 and θ_p = π/2.

Fig. 3 Molecular model for crystal class 2.

Note that assuming coupling between all oscillators, i.e., adding g₁₃ + g₂₄ with c₁₃ = c₃₁ = c₅ and c₂₄ = c₄₂ = c₆ to eqn (52), would yield a tensor of the same form as eqn (47), resulting in an over-parameterised model.

Class 222. This orthorhombic class has three mutually orthogonal 2-fold axes coinciding with the three coordinate axes. Its tensor is diagonal and has three non-zero elements.

To obtain it, two Kuhn's pairs, eqn (48)–(51) from the class-2 model, can be used by adding extra symmetries. The second oscillator within each pair is oriented symmetrically to the first with respect to the x-axis, i.e. φ₂ = −φ₁ and θ₂ = π − θ₁ are used in eqn (48b) and (50b). Similar to that in the class 2 model, the oscillators within each pair and the first and second oscillators from the opposite pairs are coupled (so that all oscillators are coupled “peripherally”). Fig. 4 exemplifies the configuration. The tensor is represented by eqn (52) and (38) with c₁₂ = c₂₁ = c₃₄ = c₄₃ = c₁ and c₁₄ = c₄₁ = c₂₃ = c₃₂ = c₂ (because d₁ = d₂ = d₃ = d₄, A₁₂ = A₃₄ and A₂₃ = A₁₄). By fixing three of the six model parameters, the number of non-zero tensor elements is matched. The polar angle θ₁ and two coupling coefficients c₁ and c₂ are independent whereas φ₁ = φ_p = θ_p = π/4 are fixed.

Fig. 4 Molecular model for crystal class 222.

Class mm2. This orthorhombic class features a single 2-fold axis along the z-axis and two mirror planes coinciding with xz- and yz-planes, respectively. Its tensor is off-diagonal, exhibits no definite symmetry and has only two non-zero elements.

Similar to that in orthorhombic class 222, additional symmetries can be imposed on the two Kuhn's pairs from class 2, given by eqn (48)–(51), to model the tensor. The second oscillator from the pair is the mirror image of the first with respect to the xz-plane, as shown in Fig. 5, i.e. θ₂ = θ₁ and φ₂ = −φ₁ are used in eqn (48b) and (50b). Furthermore, the centres of all four oscillators lie in the xy-plane, and therefore, θ_p = π/2 is substituted into eqn (49) and (51). The tensor is obtained from eqn (52) and (38) when all coupling coefficients equal c_c (because the oscillators are equidistant from the origin). Two model parameters can be fixed, θ₁ = φ_p = π/4, while the two independent ones are the azimuth φ₁ and the common coupling coefficient c_c.

Fig. 5 Molecular model for crystal class mm2.

Uniaxial crystals

Class 4̄. This tetragonal class has a single inversion 4-fold axis along the z-axis. Its tensor is traceless symmetric:and depends on two parameters. Similar to that in the two orthorhombic classes 222 and mm2, the above tensor can be modelled by imposing additional symmetries on the two Kuhn's pairs from class 2, described by eqn (48)–(51). The second oscillators from each pair are oriented and located so as to satisfy both π/2-rotation and inversion, as shown in Fig. 6, i.e., φ₂ = φ₁ − π/2 and θ₂ = π − θ₁ are substituted into eqn (49b) and (51b). Their position directions are taken to be:

v₂ = [sin θ_psin φ_p −sin θ_pcos φ_p −cos θ_p]^T (56a)

in the first pair and

v₄ = [−sin θ_psin φ_p sin θ_pcos φ_p −cos θ_p]^T (56b)

in the second one. The tensor is given by eqn (52), with all coupling coefficients, equal to c_c as they are equidistant from the origin, substituted in eqn (38). The independent model parameters are the azimuth φ₁ and the common coupling coefficient c_c, while the fixed ones are θ₁ = π/4, φ_p = 0 and θ_p = π/2.

Fig. 6 Molecular model for crystal class 4̄.

Class 4̄2m. Besides the single inversion 4-fold axis along the z-axis, this class has two orthogonal 2-fold axes perpendicular to the 4-fold axis and two mirror planes containing the 4-fold axis and bisecting the two 2-fold axes. Its tensor is traceless diagonal:

and depends only on a single parameter. This class can be treated as a higher-symmetry special case of class 4̄. It is therefore sufficient to use φ₁ = 0, θ₁ = π/4 and φ_p = θ_p = π/2 in the model equations for class 4̄ from the previous subsection. The configuration is shown in Fig. 7. The only model parameter is the common coupling coefficient c_c.

Fig. 7 Molecular model for crystal class 4̄2m.

It should be noted that one may also use φ_p = 0 instead of φ_p = π/2. This would yield the alternative traceless symmetric form of the tensor.¹⁶

It is straightforward to show that the two tensors given by eqn (57a) and (57b) transform into one another via a π/4-angle rotation about the z-axis.

Classes 3, 4 and 6. These classes have a single n-fold axis along the z-axis. Their tensors are of the following form:

and therefore, they depend on three parameters. These classes can be modelled by sets of n uncoupled Kuhn's pairs. We shall consider the case n = 3, with the other two being analogous. Two arbitrarily oriented oscillators forming a Kuhn's pair, given by eqn (48), with positions directed symmetrically with respect to the x-axis, given by eqn (49), are consecutively rotated twice about the z-axis through 2π/3 by the rotation matrix,yielding two more pairs, u₃ = R₃u₁, u₄ = R₃u₂, and v₃ = R₃v₁, v₄ = R₃v₂ as well as u₅ = R₃u₃, u₆ = R₃u₄, v₅ = R₃v₃, and v₆ = R₃v₄. The configuration is represented in Fig. 8. The tensor is given by the contributions of the initial and rotated pairs,

g₁₂ + g₃₄ + g₅₆ (60)

with c₁₂ = c₃₄ = c₅₆ = c₁ and c₂₁ = c₄₃ = c₆₅ = c₂ (because d₁ = d₃ = d₅, d₂ = d₄ = d₆ and A₂₃ = A₃₄ = A₅₆) substituted into eqn (38). The independent model parameters are the two polar angles θ₁ and θ₂ together with the coupling coefficient c₁, whereas the fixed ones are φ₁ = θ_p = π/4, φ_p = π/2, φ₂ = −φ₁ and c₂ = c₁ (because all six oscillators can be assumed to be equidistant from the origin).

Fig. 8 Molecular model for crystal class 3.

Similar to class 2, assuming additional couplings between oscillators, does not change the form of the tensor but only over-parameterises it.

Classes 32, 422 and 622. Besides the single n-fold axis along the z-axis, these classes have n evenly spaced coplanar 2-fold axes normal to the n-fold axis. Their tensors are diagonal:

and depend on two parameters. They can be modelled as higher-symmetry special cases of the n-classes by adding extra symmetries. In the case n = 3, the second oscillator from the initial Kuhn's pair is considered symmetric to the first with respect to the x-axis (see Fig. 9) so that φ₂ = −φ₁ and θ₂ = π − θ₁ are put into eqn (48b) and (50b). The tensor is obtained from eqn (60) and (38) with all coupling coefficients put equal to c_c (because all oscillators are equidistant from the origin). The independent model parameters are the polar angle θ₁ and the common coupling coefficient c_c, while the fixed ones are φ₁ = θ_p = π/4 and φ_p = π/2.

Fig. 9 Molecular model for crystal class 32.

Classes 3m, 4mm and 6mm. Besides the single n-fold axis along the z-axis, these classes have n evenly spaced mirror planes containing the n-fold axis. Their tensors are antisymmetric:

and depend only on a single parameter. In the case n = 3, an oscillator lying on the xz-plane:

u₁ = [sin θ₁ 0 cos θ₁]^T (63)

and located along the x-axis, v₁ = [1 0 0]^T, is consecutively rotated twice about the z-axis through 2π/3 by the rotation matrix R₃ given by eqn (59), yielding two more oscillators, u₂ = R₃u₁, v₂ = R₃v₁ and u₃ = R₃u₂, v₃ = R₃v₂, as shown in Fig. 10. The tensor of the triplet is given by eqn (43) and (38), all three coupling coefficients being equal to a common value c_c (because d₁ = d₂ = d₃). The only model parameter is the coupling coefficient c_c, and θ₁ = π/4 is fixed.

Fig. 10 Molecular model for crystal class 3m.

Isotropic crystals

Class 23. This isotropic class presents tetrahedral symmetry: three 2-fold axes along the three coordinate axes together with four 3-fold axes. Six oscillators are oriented along the edges of a regular tetrahedron, u₁ = a[0 1 1]^T, u₂ = a[0 −1 1]^T, u₃ = a[−1 0 1]^T, u₄ = a[1 0 1]^T, u₅ = a[1 −1 0]^T and u₆ = a[−1 −1 0]^T, where , and are located at mid-points between the summits so that v₁ = [1 0 0]^T, v₂ = [−1 0 0]^T, v₃ = [0 1 0]^T, v₄ = [0 −1 0]^T, v₅ = [0 0 1]^T and v₆ = [0 0 −1]^T, as shown in Fig. 11. The tensor is given by the contributions of the three Kuhn's pairs consisting of opposite coupled oscillators with equal pairwise coupling coefficients, i.e.eqn (60) is applied with all c_ij equalling c_c (because all oscillators are equidistant from the origin). It reduces to a scalar proportional to the coupling coefficient c_c, and only molecular enantiomorphism occurs.¹

Fig. 11 Molecular model for crystal class 23.

Class 432. The octahedral (or cubic) symmetry of this isotropic class presents three 4-fold axes along the three coordinate axes, four 3-fold axes and six 2-fold axes. This class requires three couples of Kuhn's pairs, i.e. a total of twelve oscillators. The couples are located at the centres of the opposite faces of a cube. The respective oscillators from the two pairs in a couple are perpendicular in orientations (“cross” structure); the two triplets of pairs thus formed have 222-symmetry (“2 × 3 tilted crosses”). This configuration is represented in Fig. 12. The tensor, evaluated using eqn (37), reduces to a scalar proportional to the common coupling coefficient c_c. Similar to the other optically active class 23, only molecular enantiomorphism occurs.¹

Fig. 12 Molecular model for crystal class 432.

Table 1 summarizes the derived tensors for the molecular models discussed. Note that it coincides with the complete table of optical activity tensors in crystallography.¹

Table 1 Optically active crystal classes in international (and in Schönflies) notations and their associated OA tensors

Crystal system	Crystal class point symmetry	OA tensor	Crystal system	Crystal class point symmetry	OA tensor
Triclinic (biaxial)	1 (C1)		Trigonal, tetragonal, hexagonal (uniaxial)	3, 4, 6 (C3, C4, C6)
Monoclinic (biaxial)	2 (C2)			32, 422, 622 (D3, D4, D6)
	m (Cs≡C1h)			3m, 4mm, 6mm (C3v, C4v, C6v)
Orthorhombic (biaxial)	222 (D2≡V)			4̄ (S4)
	mm2 (C2v)			4̄2m (D2d≡Vd)
			Cubic (isotropic)	23, 432 (T, O)

---

Conclusions

By combining the phenomenological description of OA with the inhomogeneous forms of electromagnetic constitutive relations, it was possible to construct microscopic models based on systems of coupled oscillators for optically active media and structures. Specifically, the evaluation of the OA-induced polarisation and magnetization of a system of coupled oscillators obeying the symmetries of an ordered medium of interest allowed to derive parameterised forms of the gyration and permittivity tensors of the medium. In particular, the coupled oscillator model of the gyration tensor enabled to establish a convenient list of properties common to all optically active media. The proposed approach is very general: the only condition for its applicability is knowing the symmetries of the ordered medium. To illustrate it, parameterised models of all optically active crystal classes were constructed. We believe that this contribution is of interest to modelling theoreticians as it provides them with microscopic models for arbitrary ordered structures and to experimentalists willing to describe the optical responses of optically active media through their gyration and permittivity tensors in terms of a minimal set of meaningful microscopic parameters.

Author contributions

R. Ossikovski performed the analytical calculations and parameterised the molecular models. O. Arteaga performed the symbolic calculations and designed the figures. The authors jointly established the molecular models for the optically active crystal classes.

Data availability

All the data necessary for the analytical calculations and to parameterize the molecular models are provided in the manuscript. The Mathcad symbolic calculations of the molecular models of the optically active crystal classes and the code for figure plotting are freely available from O. Arteaga upon request.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was funded by projects TED2021-129639B-I00, CNS2022-136051, PID2022-138699OB-I00 (Ministerio de Ciencia e Innovación, Spain) and European Metrology Programme for Innovation and Research Grant No. 20IND04ATMOC (10.13039/100014132).

References

1. R. Ossikovski, O. Arteaga and C. Sturm, Adv. Photonics Res., 2021, 2100160.
2. J. K. O'Loane, Chem. Rev., 1980, 80, 41–61.
3. J.-B. Biot, Mémoires de l'Académie de Sciences de l'Institut, 1817, vol. 2, pp. 114.
4. A. Fresnel, in Oeuvres complètes d'Augustin Fresnel, ed. H. de Sénarmont, E. Verdet and L. Fresnel, Imprimerie Impériale, Paris, 1866.
5. P. Drude, Gött. Nachr., 1904, 1, 1–8; P. Drude, The Theory of Optics, Longmans, New York, 1902.
6. W. Voigt, Gött. Nachr., 1903, 4, 155–185.
7. M. Born, Ann. Phys., 1918, 55, 177–240; M. Born, Optik, Springer, Berlin, 1972.
8. W. Kuhn, Z. Phys. Chem., Abt. B, 1929, 4, 14–36; W. Kuhn, Trans. Faraday Soc., 1930, 26, 293–308.
9. S. Chandrasekhar, Proc. – Indian Acad. Sci., Sect. A, 1953, 37, 468–484; S. Chandrasekhar, Proc. R. Soc. A, 1961, 259, 531–553.
10. V. Vyšín, Proc. Phys. Soc., 1966, 87, 55–60; V. Vyšín, Opt. Commun., 1970, 1, 307–310.
11. I. Vyšín, J. Říha and K. Sváčková, J. Opt. A:Pure Appl. Opt., 2006, 8, 584–593; I. Vyšín and J. Říha, Optik, 2012, 123, 1568–1579; I. Vyšín and J. Říha, Optik, 2016, 127, 1335–1344.
12. Yu. P. Svirko and N. I. Zheludev, Polarisation of Light in Non-Linear Optics, J. Wiley, New York, 1998.
13. T. W. H. Oates, T. Shaykhutdinov, T. Wagner, A. Furchner and K. Hinrichs, Adv. Mater., 2014, 26, 7197–7201; M. S. Davis, W. Zhu, J. Strait, J. K. Lee, H. J. Lezec, S. Blair and A. Agrawal, Nano Lett., 2020, 20, 3656–3662; X. Yin, M. Schäferling, B. Metzger and H. Giessen, Nano Lett., 2013, 13, 6238–6243; H. Kurosawa and S. Inoue, Phys. Rev. A, 2018, 98, 053805; W. Wu, W. Liu, Z. Chun, Y. Ling, J. Ding, X. Wang, L. Huang and H. Li, Opt. Express, 2020, 20, 29496–29512.
14. M. A. Kuntman, E. Kuntman, J. Sancho-Parramon and O. Arteaga, Phys. Rev. B, 2018, 98, 045410; M. S. Davis, W. Zhu, J. K. Lee, H. J. Lezec and A. Agrawal, Sci. Adv., 2019, 5, eaav8262; X. Gui, M. Hentschel, J. Defrance, J. Krauth, T. Weiss and H. Giessen, ACS Photonics, 2019, 6, 3306–3314; L. A. McCarthy, K. W. Smith, S. A. Hosseini, L. Bursi, A. Alabastri, W. S. Chang, P. Nordlander and S. Link, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 16143–16148; N. Karimitari, B. Ai and Y. Zhao, J. Phys. D: Appl. Phys., 2022, 55, 325104.
15. Y. Zhao, A. Galiautdinov and J. Tie, Nanomaterials, 2024, 14, 270.
16. L. D. Landau and E. M. Lifshitz, Electrodynamics of Continuous Media, Pergamon Press, Oxford, 1960.
17. S. F. Boys, Proc. R. Soc. London, 1934, 144, 655–675.
18. E. U. Condon, Rev. Mod. Phys., 1937, 9, 432–457.
19. F. I. Fedorov, Sov. Phys. Uspekhi, 1973, 15, 849–851.
20. W. Pauli, Pauli's Lectures on Physics. Volume 2. Optics and the Theory of Electrons, Dover, New York, 2000.
21. B. V. Bokut’ and A. N. Serdyukov, Sov. Phys. JETP, 1972, 34, 962–964.
22. R. Ossikovski and O. Arteaga, Opt. Lett., 2017, 42, 3690–3693.
23. C. Sturm, V. Zviagin and M. Grundmann, Opt. Lett., 2019, 44, 1351–1354.
24. H. A. Lorentz, The Theory of Electrons, Dover, New York, 2011.
25. L. Rosenfeld, Z. Phys., 1928, 52, 161–174; H. Nakano and H. Kimura, J. Phys. Soc. Jpn., 1969, 27, 519–535.
26. E. B. Graham and R. E. Raab, in Advances in Complex Electromagnetic Materials, ed. A. Priou, et al., Kluwer, 1997.
27. D. J. Caldwell and H. J. Eyring, Theory of Optical Activity, Wiley, New York, 1971.
28. R. E. Raab and O. L. De Lange, Multipole Theory in Electromagnetism, Clarendon Press, Oxford, 2005.
29. F. N. H. Robinson, Macroscopic Electromagnetism, Pergamon Press, Oxford, 1973.
30. L. Rosenfeld, Theory of Electrons, Dover, New York, 1965.
31. D. Franta, J. Appl. Phys., 2020, 127, 223101; D. Franta and J. Vohánka, J. Opt. Soc. Am. B, 2021, 38, 553–561; D. Franta, J. Vohánka and B. Hroncová, J. Opt. Soc. Am. B, 2023, 40, 2928–2941; B. Hroncová, D. Franta, J. Dvořák and D. Pavliňák, J. Opt. Soc. Am. B, 2023, 40, 3209–3220.

---