Cosmological Birefringence and the Geometric Phase of Photons

TL;DR

This paper demonstrates that cosmological birefringence results from an adiabatic geometric phase in photons caused by axion interactions, and proposes an analogue system using bi-isotropic media for simulation.

Contribution

It reveals the geometric phase origin of cosmological birefringence and establishes an analogue system for experimental simulation.

Findings

Cosmological birefringence linked to geometric phase in photon states.

Equivalence between axion electrodynamics and bi-isotropic media.

Potential for laboratory simulation of cosmological effects.

Abstract

Regarding axion electrodynamics in the background flat FRW universe, we show that cosmological birefringence arises from an adiabatic noncyclic geometric phase that appears in the quantum state of photons because of their interaction with the axion field. We also show that the axion electrodynamics is equivalent to standard electrodynamics in time-dependent bi-isotropic magnetoelectric Tellegen media, which serves as an analogue system that can simulate cosmological birefringence.

Equations

$$\mathcal{L}=-\sqrt{-g}\left(\frac{1}{4}{F}_{\mu \nu }{F}^{\mu \nu }+\frac{\beta }{2M}\varphi F_{\mu \nu }{\tilde{F}}^{\mu \nu }\right)$$

$$d{s}^2=a^2(\eta )\left(-d{\eta }^2+d{\mathbf{x}}^2\right)$$

$$\mathcal{L}=\frac{1}{2}\left({\mathbf{E}}^2-{\mathbf{B}}^2\right)+\frac{2\beta }{M}\varphi \mathbf{E}\cdot \mathbf{B}$$

$$\mathbf{A}(\eta ,\mathbf{x})=\int \frac{d^{3}k}{(2\pi {)}^3}\sum _{\lambda }{\mathbf{e}}_{\mathbf{k}}^{(\lambda )}{\chi }_{\mathbf{k}}^{(\lambda )}(\eta )e^{i\mathbf{k}\cdot \mathbf{x}}$$

$$\mathcal{L}_{\bm{k}}^{(\lambda)}=\frac{1}{2}\biggl{(}\dot{\bm{\chi}}_{\bm{k}}^{(\lambda)T}\dot{\bm{\chi}}_{\bm{k}}^{(\lambda)}-\omega_{k}^{2}\bm{\chi}_{\bm{k}}^{(\lambda)T}\bm{\chi}_{\bm{k}}^{(\lambda)}-\frac{4\beta\lambda\omega_{k}}{M}\phi\dot{\bm{\chi}}_{\bm{k}}^{(\lambda)T}\bm{\chi}_{\bm{k}}^{(\lambda)}\biggr{)}$$

$${\mathcal{H}}_{\mathbf{k}}^{(\lambda )}={\mathbf{\pi }}_{\mathbf{k}}^{(\lambda )T}{\dot{\mathbf{\chi }}}_{\mathbf{k}}^{(\lambda )}-{\mathcal{L}}_{\mathbf{k}}^{(\lambda )},{\mathbf{\pi }}_{\mathbf{k}}^{(\lambda )T}=\partial {\mathcal{L}}_{\mathbf{k}}^{(\lambda )}/\partial {\dot{\mathbf{\chi }}}_{\mathbf{k}}^{(\lambda )}=\left({\pi }_{\mathbf{k}R}^{(\lambda )}{\pi }_{\mathbf{k}I}^{(\lambda )}\right)$$

$${\hat{\mathcal{H}}}_{\mathbf{k}}^{(\lambda )}=\frac{1}{2}\left[{\hat{\mathbf{\pi }}}_{\mathbf{k}}^{(\lambda )T}{\hat{\mathbf{\pi }}}_{\mathbf{k}}^{(\lambda )}+{\omega }_k^2\left(1+\frac{4{\beta }^2}{M^2}{\varphi }^2\right){\hat{\mathbf{\chi }}}_{\mathbf{k}}^{(\lambda )T}{\hat{\mathbf{\chi }}}_{\mathbf{k}}^{(\lambda )}+\frac{2\beta \lambda {\omega }_k}{M}\varphi \left({\hat{\mathbf{\pi }}}_{\mathbf{k}}^{(\lambda )T}{\hat{\mathbf{\chi }}}_{\mathbf{k}}^{(\lambda )}+{\hat{\mathbf{\chi }}}_{\mathbf{k}}^{(\lambda )T}{\hat{\mathbf{\pi }}}_{\mathbf{k}}^{(\lambda )}\right)\right]$$

$${\hat{\mathcal{H}}}_{\mathbf{k}}^{(\lambda )}=\left({\widehat{\tilde{\mathbf{a}}}}_{\mathbf{k}}^{(\lambda )T}{\hat{\mathbf{a}}}_{\mathbf{k}}^{(\lambda )}+1\right){\omega }_k$$

$${\hat{\mathbf{a}}}_{\mathbf{k}}^{(\lambda )}(\eta )=\frac{1}{\sqrt{2{\omega }_k}}\left[{\omega }_k{\hat{\mathbf{\chi }}}_{\mathbf{k}}^{(\lambda )}+i\left({\hat{\mathbf{\pi }}}_{\mathbf{k}}^{(\lambda )}+\frac{2\beta \lambda {\omega }_k}{M}\varphi {\hat{\mathbf{\chi }}}_{\mathbf{k}}^{(\lambda )}\right)\right]$$

$$\left[{\hat{a}}_{\mathbf{k}i}^{(\lambda )},{\hat{a}}_{\mathbf{k}j}^{(\lambda )†}\right]$$

$$({\hat{a}}_{\mathbf{k}1}^{(\lambda )†}{\hat{a}}_{\mathbf{k}1}^{(\lambda )}+{\hat{a}}_{\mathbf{k}2}^{(\lambda )†}{\hat{a}}_{\mathbf{k}2}^{(\lambda )})|n_{\mathbf{k}}^{(\lambda )},\eta ⟩$$

$${\hat{\mathcal{H}}}_{\mathbf{k}}^{(\lambda )}|n_{\mathbf{k}}^{(\lambda )},\eta ⟩=\left(n_{\mathbf{k}}^{(\lambda )}+1\right){\omega }_k|n_{\mathbf{k}}^{(\lambda )},\eta ⟩.$$

$${\hat{I}}_{\mathbf{k}}^{(\lambda )}={\widehat{\tilde{\mathbf{b}}}}_{\mathbf{k}}^{(\lambda )T}{\hat{\mathbf{b}}}_{\mathbf{k}}^{(\lambda )}+1$$

$$\hat{\bm{b}}_{\bm{k}}^{(\lambda)}(\eta)=\frac{1}{\sqrt{2}}\biggl{\{}\frac{1}{\rho_{k}^{(\lambda)}}\hat{\bm{\chi}}_{\bm{k}}^{(\lambda)}+i\biggl{[}\rho_{k}^{(\lambda)}\left(\hat{\bm{\pi}}_{\bm{k}}^{(\lambda)}+\frac{2\beta\lambda\omega_{k}}{M}\phi\hat{\bm{\chi}}_{\bm{k}}^{(\lambda)}\right)-\dot{\rho}_{k}^{(\lambda)}\hat{\bm{\chi}}_{\bm{k}}^{(\lambda)}\biggl{]}\biggr{\}}$$

$${\ddot{\rho }}_k^{(\lambda )}+{\omega }_k^2\left(1-\frac{2\beta \lambda }{M{\omega }_k}\dot{\varphi }\right){\rho }_k^{(\lambda )}=({\rho }_k^{(\lambda )}{)}^{-3}.$$

$${\hat{I}}_{\mathbf{k}}^{(\lambda )}|{\bar{n}}_{\mathbf{k}}^{(\lambda )},\eta ⟩=\left({\bar{n}}_{\mathbf{k}}^{(\lambda )}+1\right)|{\bar{n}}_{\mathbf{k}}^{(\lambda )},\eta ⟩.$$

$$\frac{d}{d\eta }{\gamma }_{{\bar{n}}_{\mathbf{k}}}^{(\lambda )}(\eta )=⟨{\bar{n}}_{\mathbf{k}},\eta |i{\partial }_{\eta }-{\hat{\mathcal{H}}}_{\mathbf{k}}^{(\lambda )}|{\bar{n}}_{\mathbf{k}},\eta ⟩.$$

$${\gamma }_{{\bar{n}}_{\mathbf{k}}}^{(\lambda )}({\eta }_0)=-\left({\bar{n}}_{\mathbf{k}}^{(\lambda )}+1\right){\int }_0^{{\eta }_0}d\eta ({\rho }_k^{(\lambda )}{)}^{-2}.$$

$$({\rho }_k^{(\lambda )}{)}^{-2}={\omega }_k-\frac{\beta \lambda }{M}ϵ\frac{d\varphi }{d\tau }+O({ϵ}^2)\approx {\omega }_k-\frac{\beta \lambda }{M}\dot{\varphi }.$$

$${\gamma }_{n_{\mathbf{k}}}^{(\lambda )}({\eta }_0)=\left(n_{\mathbf{k}}^{(\lambda )}+1\right)\left(\frac{\beta \lambda }{M}\Delta \varphi -{\omega }_k{\eta }_0\right),\Delta \varphi ({\eta }_0)=\varphi ({\eta }_0)-\varphi (0).$$

$${\Gamma }_{n_{\mathbf{k}}}^{\text{geom}}({\eta }_0)=\frac{\beta }{M}\left(n_{\mathbf{k}}^{(+)}-n_{\mathbf{k}}^{(-)}\right)\Delta \varphi ({\eta }_0).$$

$${\nabla }_{\mu }{F}^{\mu \nu }+\xi {\nabla }_{\mu }\left(\varphi {\tilde{F}}^{\mu \nu }\right)$$

$${\nabla }_{\mu }{\tilde{F}}^{\mu \nu }$$

$$E_i=F_{i0},D_i=\sqrt{-g}\left(F^{0i}+\xi \varphi {\tilde{F}}^{0i}\right)$$

$$B_i=\frac{1}{2}{ϵ}_{ijk}{F}_{jk},H_i=\frac{1}{2}{ϵ}_{ijk}\sqrt{-g}\left(F^{jk}+\xi \varphi {\tilde{F}}^{jk}\right).$$

$$\nabla \cdot \mathbf{D}=\nabla \cdot \mathbf{B}=0,-{\partial }_0\mathbf{D}+\nabla \times \mathbf{H}={\partial }_0\mathbf{B}+\nabla \times \mathbf{E}=0$$

$$\mathbf{D}$$

$$\mathbf{H}$$

$$\mathbf{D}$$

$$\mathbf{B}$$

Full text

axion electrodynamics, cosmological birefringence, geometric phase

I Introduction

Dark matter, which neither emits nor scatters electromagnetic (EM) radiation, is considered to be one of the main cosmic constituents of the universe, comprising some 23% of its energy. The components of dark matter are unknown, but there are several plausible candidates. One interesting class is furnished by massive pseudo-Goldstone bosons (axions and axion-like particles) Kolb1990 ; Marsh2016 ; Leon2017 ; pseudo-scalar particles first proposed in order to resolve the strong CP problem in quantum chromodynamics Peccei1977 . A distinguishing feature of the axion field is that it couples to EM field, notably the cosmic microwave background (CMB) radiation. As a subject of growing interest, axion-electrodynamics, i.e., the standard electrodynamics modified by photon-axion interaction, provides a theoretical framework for possible violation of parity and Lorentz invariance Carroll1990 .

CMB radiation furnishes the best testing ground for axion-electrodynamics. Its polarization arises from the asymmetric Thomson scattering in the epoch of recombination, yielding the longest available look-back time as compared to any other cosmological source. Current measures and constraints on the polarization pattern of CMB anisotropies produce an upper limit on the linear polarization rotation angle of the order of few degrees (see, e.g., Wu2009 ; Chiang2010 ; Jarosik2011 ). The EM-axion field interaction affects the polarization of EM waves propagating over cosmological distances Carroll1990 ; Carroll1991 ; Harari1992 ; Finelli2009 . In this connection, Carroll et al. Carroll1990 ; Carroll1991 were first to show that different polarizations of EM wave propagate in the same direction with different phases, causing a rotation of the polarization plane (cosmological birefringence). Comparison of the polarization rotation (which depends on the axion field features) with the CMB polarization observations, can give experimental confirmation of the existence (or otherwise) of axions and their physical properties.

Because of the dependence on polarization of optical and transverse acoustical wave propagation in inhomogeneous media, there occurs similar polarization rotation (the Rytov rotation), which is established to be a manifestation of the Berry phase of the quanta (photons/phonons) of the quantized wave fields Torabi2012 ; Torabi2008 ; Mehrafarin2009 . In this work, we show that the cosmological birefringence likewise arises from an adiabatic noncyclic geometric phase that appears in the quntum state of the photons because of their interaction with the slowly varying pseudo-scalar field. This phase, which is gauge and reparametrization invariant, is associated with the image of the path traversed by the state vector of photons in the projective Hilbert space. A noncyclic geometric phase also has been shown to characterize the mixing phenomenon of photons and axions in the presence of strong magnetic fields Capolupo .

Berry, in his seminal work Berry1984 , showed that in a cyclic adiabatic quantum evolution, the wave function ends up in the original state up to a phase factor that is the sum of a dynamical phase and an additional geometric phase, known as the Berry phase. This phase has been generalized to the case of nonadiabatic Aharonov1987 , nonunitary, and noncyclic evolutions Samuel1988 . Also, using a quantum kinematic approach and the Pancharatnam criterion Pancharatnam1956 , Mukunda and Simon showed that a gauge and reparametrization invariant geometric phase can be naturally associated with any smooth open curve traversed by unit vectors in Hilbert space Mukunda1993 . Since its discovery, examples of geometric phase have appeared in many different areas of physics Chruscinski2004 ; Baggio2017 . It is now well known that geometric phases arise in systems with time-dependent Hamiltonians and non-stationary quantum states Zeng1995 .

Generalized time-dependent harmonic oscillator (TDHO) is a ubiquitous time-dependent quantum systems, which is apt for the application of the dynamical invariant method introduced by Lewis and Riesenfeld (LR) Lewis1969 . In the present work, we show that axion electrodynamics in flat Friedmann-Robertson-Walker (FRW) universe, when quantized, can be described by a collection of generalized TDHO Hamiltonians for the photons. We use the LR invariants to find the phase of the quantum state of photons interacting with the axion field. By analyzing the phase factor, we then establish that cosmological birefringence arises from the geometric part of the phase.

Finally, we present the axion electrodynamics from another perspective. We propose a similar system where the relevant features of photon-axion mixing in expanding spacetime can be reproduced analogously. This approach, known as analogue gravity Barcelo2011 , is the study of phenomena in curved spacetimes by means of analogue systems whose kinematics are governed by the same mathematical formalisms. With analogue models, it is possible to simulate predicted phenomena that have not yet been experimentally observed. They, therefore, provide a powerful tool for testing cosmological models. The analogue model that we propose to use for simulating cosmological birefringence is light traveling in time-dependent Tellegen medium Tellegen1948 . We show that the propagation of light in time-dependent bi-isotropic magnetoelectric media shares the same behavior as in axion electrodynamics in flat FRW universe. The results obtained from this simulation can, thus, provide a substitution for direct experimental observation.

The organization of this paper is as follows. In section II, we construct the Lagrangian of axion electrodynamics in flat FRW universe and derive the resulting generalized TDHO Hamiltonian for the photons via canonical quantiztion. In section III, we obtain the phase of the photon quantum state by means of the LR invariant method. We then analyze this phase and derive the polarization rotation from the geometric part of the phase. Finally, in section IV, we present the analogue system.

II Canonical quantization of axion electrodynamics in flat FRW spacetime

The Lagrangian density for axion electrodynamics in curved spacetime is ($\hbar=c=1$)

[TABLE]

where $M$ (mass scale for broken global symmetry in Peccei-Quinn model Peccei1977 ) has mass dimension 1, $\beta$ is a dimensionless coupling, $\phi$ is the pseudo-scalar axion field, $F^{\mu\nu}=\nabla^{\mu}A^{\nu}-\nabla^{\nu}A^{\mu}$ is the EM field tensor, and $\tilde{F}^{\mu\nu}=\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}F_{\rho\sigma}$ is its dual. $\epsilon^{\mu\nu\rho\sigma}$ is the complete antisymmetric tensor, which is related to the absolutely antisymmetric Levi-Civita symbol $[\mu\nu\rho\sigma]$ through the metric according to: $\epsilon_{\mu\nu\rho\sigma}=\sqrt{-g}\,[\mu\nu\rho\sigma]$, $\epsilon^{\mu\nu\rho\sigma}=-(\sqrt{-g})^{-1}[\mu\nu\rho\sigma]$. For a spatially flat FRW universe the line element reads:

[TABLE]

where $a$ is the scale factor and $\eta$ is the conformal time. In the Coulomb gauge $A^{\mu}=(0,\bm{A})$ (the vector potential $\bm{A}$ satisfying $\nabla\cdot\bm{A}=0$), (1) becomes

[TABLE]

where $\bm{E}=-\partial_{\eta}\bm{A}$ and $\bm{B}=\nabla\times\bm{A}$ are the electric and magnetic fields, respectively. We write the Fourier transform of the vector potential as

[TABLE]

where $\lambda=\pm$ runs over the two circular polarization states, $\bm{e}_{\bm{k}}^{(\lambda)}$ represent the corresponding polarization vectors, and $\bm{e}_{\bm{k}}^{(\lambda)\star}=\bm{e}_{-\bm{k}}^{(\lambda)},\,\chi_{\bm{k}}^{(\lambda)\star}=\chi_{-\bm{k}}^{(\lambda)}$. Substituting (4) in the Lagrangian $L=\int\mathcal{L}d^{3}x$ and forming the row matrix $\bm{\chi}_{\bm{k}}^{(\lambda)T}=\left(\chi_{\bm{k}R}^{(\lambda)}\quad\chi_{\bm{k}I}^{(\lambda)}\right)$ from the real and imaginary parts of $\chi_{\bm{k}}^{(\lambda)}$, the Lagrangian for each Fourier mode and polarization state becomes ($L=\int\sum_{\lambda}\mathcal{L}_{\bm{k}}^{(\lambda)}(2\pi)^{-3}d^{3}k$)

[TABLE]

where dot denotes differentiation with respect to conformal time, and $\omega_{k}=|\bm{k}|$. The corresponding Hamiltonian is given by

[TABLE]

$\bm{\pi}_{\bm{k}}^{(\lambda)}$ being the conjugate momenta. Upon canonical quantization (operators represented by “hat”), we have

[TABLE]

which is the Hamiltonian of a generalized TDHO of frequency $\omega_{k}$.

This Hamiltonian can be written as

[TABLE]

where the two-component column matrix

[TABLE]

and $\hat{\tilde{\bm{a}}}_{\bm{k}}^{(\lambda)T}=\left(\hat{a}_{\bm{k}1}^{(\lambda){\dagger}}\quad\hat{a}_{\bm{k}2}^{(\lambda){\dagger}}\right)$ are the standard particle annihilation and creation matrix operators, respectively. They satisfy ($i,j=1,2$)

[TABLE]

$\ket{n_{\bm{k}}^{(\lambda)},\eta}$ being the Hamiltonian eigenstate given by

[TABLE]

It follows that Hamiltonian (7) describes any number $n_{\bm{k}}^{(\lambda)}$ of photons of a given wavelength and polarization, interacting with the axion field. In the next section, we use the dynamical invariant method Lewis1969 to obtain the phase of their quantum state.

III Geometric Phase and Cosmological Birefringence

The invariant operator, which by definition is a constant of motion satisfying the von-Neumann equation, has been derived for the generalized TDHO Hamiltonian in Gao1990 ; Gao1991 . For Hamiltonian (7), it takes the form

[TABLE]

where

[TABLE]

is the lowering and $\hat{\tilde{\bm{b}}}_{\bm{k}}^{(\lambda)T}=\left(\hat{b}_{\bm{k}1}^{(\lambda){\dagger}}\quad\hat{b}_{\bm{k}2}^{(\lambda){\dagger}}\right)$ is the raising matrix operator, and $\rho_{k}^{(\lambda)}(\eta)$ is an auxiliary variable that satisfies

[TABLE]

The components of the raising and lowering operators obey standard relations similar to (II) in terms of the eigenstates of the invariant operator, given by

[TABLE]

A photon state can be written as the superposition of the eigenstates of the invariant operator with phase factors (LR phase) that follow from

[TABLE]

This phase has been calculated for a generalized TDHO, with the result Gao1990

[TABLE]

In our case, $\rho_{k}^{(\lambda)}$ obeys equation (12), $\eta=0$ denotes the time when CMB photons left the last scattering surface and $\eta_{0}$ is the present time.

In the adiabatic limit of slow time variation, the auxiliary variable attains a definite form. In the auxiliary equation (12), we introduce the adiabatic parameter $\epsilon$ ($\ll 1$) and write $\tau=\epsilon\eta$ to get

[TABLE]

Adiabatic approximation, therefore, holds for $\beta\dot{\phi}/M\omega_{k}\ll 1$, i.e. for sufficiently slowly varying axion field. It follows also from differentiation of (16) that in the adiabatic limit $\dot{\rho}_{k}^{(\lambda)}\approx 0$. Hence, $\hat{\mathscr{H}}_{\bm{k}}^{(\lambda)}\approx\omega_{\bm{k}}\hat{I}_{\bm{k}}^{(\lambda)}$, implying that the eigenstates of the invariant operator coincide with those of the Hamiltonian in this limit. Thus, $\bar{n}_{\bm{k}}^{(\lambda)}$ corresponds to the number of photons of a given wavelength and polarization, $n_{\bm{k}}^{(\lambda)}$.

In adiabatic evolution, photons initially in an eigenstate of the Hamiltonian will evolve into the instantaneous eigenstate after a time $\eta_{0}$, acquiring the phase $\gamma_{n_{\bm{k}}}^{(\lambda)}(\eta_{0})$. Substitution of (16) into (15) yields

[TABLE]

The second term in above, which is immediately seen from (10) to come from the second term in (14), is the usual dynamical phase acquired along open curve $\mathcal{C}$ traversed by the state vector. For any arbitrary quantum evolution, the geometric phase is the difference between the total phase and the dynamical phase Mukunda1993 . The total phase is the relative phase of the end point of $\mathcal{C}$ with respect to its starting point, as determined by the Pancharatnam criterion Pancharatnam1956 , and the dynamical phase is a locally additive functional of $\mathcal{C}$. Their difference is a gauge and reparametrization invariant phase that is only a functional of the image of $\mathcal{C}$ in the projective Hilbert space.

Thus, for a general photon state that consists of $n_{\bm{k}}^{(+)}$ ($n_{\bm{k}}^{(-)}$) photons of positive (negative) helicity, the adiabatic noncyclic geometric phase acquired during its journey from the last scattering surface to the present epoch is given by

[TABLE]

Therefore, a photon state with net helicity will exhibit the geometric phase shift (18) as it travels through the background slowly varying axion field. This is our main result. For a linearly polarized photon, being the superposition of two opposite circularly polarized photons, this will have the effect of a rotation of the polarization plane by $\beta\Delta\phi/M$, which corresponds to the classical result of Carroll1990 ; Carroll1991 for cosmological birefringence.

IV The Analogue System

Let us consider the possibility of viewing axion electrodynamics in curved spacetime as Maxwell electrodynamics in an effective medium in flat spacetime, the properties of which medium are determined by the metric of the original curved spacetime (see, e.g.Skrotskii1957 ; Mashhoon1973 ).

The axion electrodynamics field equations obtained from Lagrangian (1) are

[TABLE]

where $\xi=2\beta/M$. Consider a background flat spacetime with Cartesian coordinates, in which the electric and magnetic fields are defined using the decompositions $F_{\mu\nu}\rightarrow(\bm{E},\bm{B})$, and $\sqrt{-g}(F^{\mu\nu}+\xi\phi\tilde{F}^{\mu\nu})\rightarrow(-\bm{D},\bm{H})$. That is, define

[TABLE]

In terms of these quantities, field equations (IV) take the form

[TABLE]

which are the source-free Maxwell equations in a material medium in flat spacetime. For flat FRW metric with cosmic time $dt=ad\eta$, the constitutive relations obtained from (20), thus, read:

[TABLE]

Rearranging these equations in the form of standard magnetoelectric equations, we have

[TABLE]

where $\varepsilon=a(1+\xi^{2}\phi^{2})$ is the electric permittivity, $\mu=a$ is the magnetic permeability, and $\alpha=-a\xi\phi$ is the magnetoelectric coefficient. Equations (IV) are the characteristic equations of bi-isotropic magnetoelectric media. Maxwell electrodynamics in this time-dependent medium presents the analogue system for Axion electrodynamics.

In 1948, Tellegen Tellegen1948 suggested that an assembly of randomly distributed parallel or antiparallel electric-magnetic dipole twins can furnish a new type of magnetoelectric material that is non-reciprocal and non-chiral. Since then, the very existence of the Tellegen media (now also known as axion media) in nature has been the subject of theoretical debate (see, e.g., Lakhtakia1994 ; Sihvola1995 ; Tretyakov2003 ). However, there is experimental evidence that the natural composite chromium sesquioxide $\text{Cr}_{2}\text{O}_{3}$ can be modelled by a uniaxial Tellegen medium Hehl2008 . Moreover, the realization of an artificial isotropic Tellegen medium has been demonstrated by Gosh et al. Gosh2007 ; Gosh2008 . In Gosh2008 , the authors have shown that by combining permanent electric and magnetic moments in particles, it is possible to construct a cross-correlation between the electric and magnetic properties of matter. These magnetoelectric particles interact with each other and their surroundings, and the nature of these interactions will determine the magnetoelectric response. By making this response decadently oscillatory (e.g. by applying an a.c. field with appropriate frequency) like the axion field Kolb1990 ; Marsh2016 , one can construct the time-dependent medium defined by the constitutive equations (IV) and study EM wave propagation through it. This provides an analogue system for the simulation of cosmological birefringence.