Cloaking and anamorphism for light and mass diffusion

TL;DR

This paper explores how transformation optics can be used to create cloaking effects for diffusive light and mass, demonstrating numerical methods and proposing new designs for external cloaks with spatially varying parameters.

Contribution

It extends transformation optics to diffusive regimes, introduces a non-linear transform to address non-physical parameters, and explores cloaking in time domain and non-Euclidean geometries.

Findings

Infinite objects of different shapes scatter diffusive light identically.

A design for an external diffusion cloak with spatially varying parameters is proposed.

Transformation methods can be applied to time-dependent diffusion equations.

Abstract

We first review classical results on cloaking and mirage effects for electromagnetic waves. We then show that transformation optics allows the masking of objects or produces mirages in diffusive regimes. In order to achieve this, we consider the equation for diffusive photon density in transformed coordinates, which is valid for diffusive light in scattering media. More precisely, generalizing transformations for star domains introduced in [Diatta and Guenneau, J. Opt. 13, 024012, 2011] for matter waves, we numerically demonstrate that infinite conducting objects of different shapes scatter diffusive light in exactly the same way. We also propose a design of external light-diffusion cloak with spatially varying sign-shifting parameters that hides a finite size scatterer outside the cloak. We next analyse non-physical parameter in the transformed Fick's equation derived in [Guenneau and Puvirajesinghe, R. Soc. Interface 10, 20130106, 2013], and propose to use a non-linear transform that overcomes this problem. We finally investigate other form invariant transformed diffusion-like equations in the time domain, and touch upon conformal mappings and non-Euclidean cloaking applied to diffusion processes.

Equations

$$\begin{array}[]{ll}&\displaystyle{\frac{1}{\nu({\bf x})}\frac{\partial}{\partial t}L({\bf x},{\bf n},t)+{\bf n}\cdot\nabla\left(L({\bf x},{\bf n},t)\right)+(\mu_{a}+\mu_{s}^{\prime})L({\bf x},{\bf n},t)}\\ &=\displaystyle{S({\bf x},{\bf n},t)+\mu_{s}^{\prime}\int dn^{\prime}L({\bf x},{\bf n}^{\prime},t)f({\bf n},{\bf n}^{\prime})}\;,\end{array}$$

$$\Phi (\mathbf{x},t)=\int dnL(\mathbf{x},\mathbf{n},t).$$

$$L(\mathbf{x},\mathbf{n},t)=\sum _{l=0}^N\sum _{m=-l}^l{\varphi }_{l,m}(\mathbf{x},t)Y_{l,m}(\mathbf{n}),$$

$$S(\mathbf{x},\mathbf{n},t)=\sum _{l=0}^N\sum _{m=-l}^l{s}_{l,m}(\mathbf{x},t)Y_{l,m}(\mathbf{n}).$$

$$f(\mathbf{n},{\mathbf{n}}^{\prime })=\sum _{l=0}^N\sum _{m=-l}^l{g}_l{Y}_{l,m}^{\ast }({\mathbf{n}}^{\prime })Y_{l,m}(\mathbf{n}).$$

$$-\nabla \cdot \left(D(\mathbf{x})\nabla \Phi (\mathbf{x},t)\right)+{\mu }_a(\mathbf{x})\Phi (\mathbf{x},t)+\frac{1}{\nu (\mathbf{x})}\frac{\partial \Phi (\mathbf{x},t)}{\partial t}-S_0=$$

$$-\frac{3D(\mathbf{x})}{\nu (\mathbf{x})}\left[{\mu }_a(\mathbf{x})\frac{\partial \Phi (\mathbf{x},t)}{\partial t}+\frac{1}{\nu (\mathbf{x})}\frac{{\partial }^2\Phi (\mathbf{x},t)}{\partial t^2}\right]+\frac{3D(\mathbf{x})}{\nu (\mathbf{x})}\frac{\partial S_0}{\partial t}-3D(\mathbf{x})\nabla \cdot S_1,$$

$$-\nabla \cdot \left(D(\mathbf{x})\nabla \Phi (\mathbf{x},t)\right)+\frac{1}{\nu (\mathbf{x})}\frac{\partial \Phi (\mathbf{x},t)}{\partial t}+{\mu }_a(\mathbf{x})\Phi (\mathbf{x},t)=S_0.$$

$$-\nabla \cdot \left(D(\mathbf{x})\nabla \Phi (\mathbf{x})\right)-\frac{i\omega }{\nu (\mathbf{x})}\Phi (\mathbf{x})+{\mu }_a(\mathbf{x})\Phi (\mathbf{x})=S_0.$$

$$\left(\begin{array}[]{c}dx\\ dy\\ dz\\ \end{array}\right)=\mathbf{J}_{xx^{\prime}}\left(\begin{array}[]{c}dx^{\prime}\\ dy^{\prime}\\ dz^{\prime}\\ \end{array}\right)\;,\hbox{with}\;\mathbf{J}_{xx^{\prime}}=\frac{\partial(x,y,z)}{\partial(x^{\prime},y^{\prime},z^{\prime})}\;.$$

$$\underline{\underline{D^{\prime }}}=D{\mathbf{T}}_T^{-1},{\nu }^{\prime }=T_{zz}\nu ,{\mu }_a^{\prime }=T_{zz}^{-1}{\mu }_a,S^{\prime }=T_{zz}^{-1}S,$$

$$-\nabla \cdot \left(\underline{\underline{D^{\prime }}}({\mathbf{x}}^{\prime })\nabla \Phi ({\mathbf{x}}^{\prime })\right)-\frac{i\omega }{{\nu }^{\prime }({\mathbf{x}}^{\prime })}\Phi ({\mathbf{x}}^{\prime })+{\mu }_a^{\prime }({\mathbf{x}}^{\prime })\Phi ({\mathbf{x}}^{\prime })=S^{\prime }({\mathbf{x}}^{\prime })\Phi ({\mathbf{x}}^{\prime }),$$

$$\begin{array}[]{lll}&\displaystyle{\int_{\Omega}d{\bf x}\left(\frac{-i\omega+\mu_{a}\nu}{\nu}\Phi\phi\right)+\int_{\Omega}d{\bf x}\left(D\nabla\Phi\cdot\nabla{\phi}\right)}\\ &\displaystyle{-\int_{\partial\Omega}ds\left(D\nabla\Phi\cdot{\bf n}{\phi}\right)+<S,\phi>}=0\;,\end{array}$$

$$\begin{array}[]{ll}&\nabla=(\partial/\partial x,\partial/\partial y,\partial/\partial z)^{T}=\displaystyle{\frac{\partial(x^{\prime},y^{\prime},z^{\prime})}{\partial(x,y,z)}(\partial/\partial x^{\prime},\partial/\partial y^{\prime},\partial/\partial z^{\prime})^{T}}\\ &={\bf J}_{x^{\prime}x}\nabla^{\prime}={\bf J}_{xx^{\prime}}^{-1}\nabla^{\prime}\end{array}$$

$$\begin{array}[]{ll}&\displaystyle{\int_{\Omega}d{\bf x^{\prime}}\left(\frac{-i\omega+\mu_{a}\nu}{\nu}\Phi\phi\hbox{det}({\bf J}_{xx^{\prime}})\right)}\\ &+\displaystyle{\int_{\Omega}d{\bf x^{\prime}}\left\{\left({\bf J}_{xx^{\prime}}^{-1}\nabla^{\prime}{\phi}\cdot D{\bf J}_{xx^{\prime}}^{-1}\nabla^{\prime}\Phi\right)\hbox{det}({\bf J}_{xx^{\prime}})\right\}}\\ &-\displaystyle{\int_{\partial\Omega}ds^{\prime}\left(D{\bf J}_{xx^{\prime}}^{-1}\nabla^{\prime}\Phi\cdot{\bf n}{\phi}\right)\hbox{det}({\bf J}_{xx^{\prime}})}\\ &+<\hbox{det}({\bf J}_{xx^{\prime}})S,\phi>=0\;.\end{array}$$

$$\begin{array}[]{ll}&\displaystyle{\int_{\Omega}d{\bf x^{\prime}}\left(\frac{-i\omega+\mu_{a}\nu}{\nu}\Phi\phi\hbox{det}({\bf J}_{xx^{\prime}})\right)}\\ &-\displaystyle{\int_{\Omega}d{\bf x^{\prime}}\left\{\nabla^{\prime}\cdot\left({\bf J}_{xx^{\prime}}^{-1}D{\bf J}_{xx^{\prime}}^{-T}\nabla^{\prime}{\Phi}\right)\phi\hbox{det}({\bf J}_{xx^{\prime}})\right\}}\\ &+<\hbox{det}({\bf J}_{xx^{\prime}})S,\phi>=0\;,\end{array}$$

$$\begin{array}[]{ll}&\underline{\underline{D^{\prime}}}=D\mathbf{T}_{T}^{-1}\;,\;\nu^{\prime}=\hbox{det}({\bf J}_{xx^{\prime}})^{-1}\nu={T}_{zz}\nu\;,\\ &\;\mu_{a}^{\prime}=\hbox{det}({\bf J}_{xx^{\prime}})\mu_{a}={T}_{zz}^{-1}\mu_{a}\;,\;S^{\prime}={T}_{zz}^{-1}S\;,\end{array}$$

$$\rho (\mathbf{x})c(\mathbf{x})\frac{\partial u}{\partial t}=\nabla \cdot (\underline{\underline{\kappa }}(\mathbf{x})\nabla u),$$

$$\displaystyle\underline{\underline{\kappa}}=\left(\begin{array}[]{ccc}\kappa_{11}&\kappa_{12}\cr\kappa_{21}&\kappa_{22}\end{array}\right)$$

$$\displaystyle\underline{\underline{\varepsilon}}({x,y})=\left(\begin{array}[]{ccc}\varepsilon_{11}&\varepsilon_{12}&0\cr\varepsilon_{21}&\varepsilon_{22}&0\cr 0&0&\varepsilon_{33}\end{array}\right)\;,\;\underline{\underline{\mu}}({x,y})=\left(\begin{array}[]{ccc}\mu_{11}&\mu_{12}&0\cr\mu_{21}&\mu_{22}&0\cr 0&0&\mu_{33}\end{array}\right)\;,$$

$$\nabla \times \left({\underline{\underline{\epsilon }}}^{-1}(x,y)\nabla \times {\mathbf{H}}_l\right)=-\frac{{\partial }^2}{\partial t^2}\left(\underline{\underline{\mu }}(x,y){\mathbf{H}}_l\right)$$

$$\nabla \times \left({\underline{\underline{\mu }}}^{-1}(x,y)\nabla \times {\mathbf{E}}_l\right)=-\frac{{\partial }^2}{\partial t^2}\left(\underline{\underline{\epsilon }}(x,y){\mathbf{E}}_l\right)$$

$$\displaystyle{\bf M}=\left(\begin{array}[]{ccc}m_{11}&m&0\cr m&m_{22}&0\cr 0&0&m_{33}\end{array}\right)=\left(\begin{array}[]{ccc}{\bf M}_{T}&0\cr 0&m_{33}\end{array}\right)\;.$$

$$\nabla\times\Biggl{(}{\bf M}\nabla\times\Bigl{(}u(x,y)\mathbf{e}_{z}\Bigr{)}\Biggr{)}=-\nabla\cdot\Biggl{(}{\bf M}_{T}^{-1}\hbox{det}({\bf M}_{T})\nabla u(x,y)\Biggr{)}\mathbf{e}_{z}\;.$$

$$\begin{array}[]{ll}\nabla\times\Biggl{(}{\bf M}\nabla\times\Bigl{(}u(x,y)\mathbf{e}_{z}\Bigr{)}\Biggr{)}&=\displaystyle{-\Biggl{(}\frac{\partial}{\partial x}\Bigl{(}m_{22}\frac{\partial u}{\partial x}-m\frac{\partial u}{\partial y}\Bigr{)}}\\ &\displaystyle{+\frac{\partial}{\partial y}\Bigl{(}m_{11}\frac{\partial u}{\partial y}-m\frac{\partial u}{\partial x}\Bigr{)}\Biggr{)}\mathbf{e}_{z}}\;.\end{array}$$

$$\displaystyle{\bf M}^{\prime}=\left(\begin{array}[]{ccc}m^{\prime}_{11}&m^{\prime}_{12}\cr m^{\prime}_{21}&m^{\prime}_{22}\end{array}\right)\;.$$

$$\nabla\times\Bigl{(}{\bf M}\nabla\times\Bigl{(}u(x,y)\mathbf{e}_{z}\Bigr{)}\Biggr{)}=-\nabla\cdot\Biggl{(}{\bf M}^{\prime}\nabla u\Biggr{)}\mathbf{e}_{z}\;,$$

$$m_{11}^{\prime }\frac{\partial u}{\partial x}+m_{12}^{\prime }\frac{\partial u}{\partial y}=m_{22}\frac{\partial u}{\partial x}-m\frac{\partial u}{\partial y}$$

$$m_{21}^{\prime }\frac{\partial u}{\partial x}+m_{22}^{\prime }\frac{\partial u}{\partial y}=-m\frac{\partial u}{\partial x}+m_{11}\frac{\partial u}{\partial y},$$

$$\displaystyle\nabla\cdot\left[\frac{1}{\varepsilon_{11}\varepsilon_{22}-\varepsilon_{12}\varepsilon_{21}}\left(\begin{array}[]{ccc}\varepsilon_{11}&\varepsilon_{12}\cr\varepsilon_{21}&\varepsilon_{22}\end{array}\right)\nabla H_{z}\right]=\frac{\partial^{2}}{\partial t^{2}}\left(\mu_{33}H_{z}\right)\;,$$

$$\displaystyle\nabla\cdot\left[\frac{1}{\mu_{11}\mu_{22}-\mu_{12}\mu_{21}}\left(\begin{array}[]{ccc}\mu_{11}&\mu_{12}\cr\mu_{21}&\mu_{22}\end{array}\right)\nabla E_{z}\right]=\frac{\partial^{2}}{\partial t^{2}}\left(\varepsilon_{33}E_{z}\right)\;.$$

Full text

1 Introduction

Following the 2006 proposal of an invisibility cloak via geometric transforms in the covariant Maxwell’s equations [1], and their conformal counterpart in specific regimes [2], enhanced control of electromagnetic and plasmonic waves has been experimentally validated in a number of fascinating studies [3, 4, 5, 6, 7, 8, 9, 10], that also include acoustic [11, 12, 13], hydrodynamic [14, 15, 16, 17, 18, 19, 20, 21], elastodynamic [22, 23, 24, 25, 26, 27, 28, 29, 30] and even seismic [31, 32] waves, see e.g. for a book devoted to these topics [33]. Importantly, transformation optics [1] and conformal optics [2] have been reconciled in the proposal of invisibility carpets [34], that make use of quasi-conformal mappings that relax the constraint on anisotropy [35, 36].

Interestingly, control of wave trajectories can be extended to matter waves [37]. Looking at the mathematical origin of cloaking, which comes from the concept of Dirichlet-to-Neumann map in inverse problems with anisotropic media [38], some of us proposed some mimetism for matter waves [39], which prompted further studies in control of diffusion phenomena that can be achieved through coordinate transformations in the Fourier [40] and Fick’s [41] equations for heat and mass respectively. The latter work has served as an inspiration for the group of Martin Wegener, see Fig. 1, who exploited the Fick’s diffusion equation for multiple light scattering usages, whereby light is governed by ballistic laws [42]: cylindrical and spherical invisibility cloaks made of thin shells of polydimethylsiloxane doped with melamine-resin microparticles surrounding a diffusively reflecting hollow core were experimentally shown to achieve good cloaking performance in a water-based diffusive medium throughout the entire visible spectrum and for all illumination conditions and incident polarizations of light. However, these authors noted that the transformed Fick’s equation as proposed in [41] does not fully retain its form under Pendry’s transform [1], since a heterogeneous isotropic factor appears in front of the time derivative (a problem partially solved in [41] using reduced parameters that preserve diffusion trajectories, but induce some impedance mismatch at the cloak’s boundaries). This issue remained an Achilles heel in the control of mass diffusion through geometric transform until now. Interestingly, scattering cancellation techniques for mantle cloaking in diffusion processes do not suffer from this pitfall [43, 44, 45] as these are based on a totally different mechanism not underpinned by coordinate changes. This cloaking route reminiscent of neutral inclusions, a concept introduced by Kerner 60 years ago [46], has been followed in a number of works in the theory of composites in the quasi-static limit [47, 48, 49]. It can also explain the success of the aforementioned cloaking experiment in a water-paint based environment in [42], see Fig. 1, which has been extended to solid media invoking the concept of neutral inclusion, and no longer transformed Fick’s equation, in [50].

In this topical review paper, we focus our analysis on cloaking of diffusive light in scattering media via geometric transforms, which is governed by the diffusive photon density wave (DPDW) equation. The transformed DPDW equation involves an anisotropic spatially varying diffusivity, as well as other heterogeneous but isotropic parameters. We shall first give a physical interpretation to all the terms in transformed DPDW. We will then consider numerical illustrations of anamorphism and external cloaking in DPDW equation. This part is actually an extension of [45] where DPDW equation was first used to reduce scattering of diffuse light, however without resorting to coordinate changes. We will then consider a non-linear transform that makes Fick’s equation totally form invariant in certain cylindrical geometries. We will finally come back to the interesting question of variation of the field inside the invisibility region in the time domain, since the linear transform which we use for Fick’s equation avoids pitfall of a vanishing factor in front of the time derivative in transformed equations for the diffusion of mass, heat, matter and light. The concluding remarks are followed by perspectives on conformal mappings for diffusion processes and other parabolic equations of interest for geometric transforms.

2 Transformed governing equations for diffusive light waves

Let us set the scene of diffusion of light in scattering media.

2.1 Introduction to the equation for propagation of diffuse photon density waves DPDW

The transport equation for photons in scattering and absorbing media is given by [51, 52, 53]

[TABLE]

with $L({\bf x},{\bf n},t)$ denoting the radiance at position ${\bf x}$, in direction ${\bf n}$, at time $t$ (with units of $Wm^{-2}sr^{-1}$). $f({\bf n},{\bf n}^{\prime})$ denotes the probability of scattering to ${\bf n}^{\prime}$ from direction ${\bf n}$. $\nu$ is the speed of light in the medium and $\mu_{s}^{\prime}$, $\mu_{a}$ are the scattering and absorption coefficients, respectively. $S$ is the source term. The photon fluence is given by:

[TABLE]

In order to get solutions of this equation for diverse geometries, one needs to make the so-called $P_{N}$ approximation [54, 55]. It consists in expanding the radiance $L$ and source $S$ in terms of spherical harmonics $Y_{l,m}$, and truncating at order $l=N$, i.e.

[TABLE]

and

[TABLE]

The phase function can also be expanded as (assuming that the scattering amplitude is only dependent on the change in direction of the photon [54]):

[TABLE]

The approximation $P_{1}$ (i.e. setting $N=1$) is actually, quite appropriate when $\mu_{a}\ll\mu_{s}^{\prime}$ (absorption much smaller than scattering) and $f$ is not too anisotropic [53, 54, 55]. Inserting Eqs. (3)-(4) and the development of $f$ in Eq. (1) yields, under the $P_{1}$ approximation

[TABLE]

with $D$ the modified diffusivity, i.e. $D=1/(3\mu_{s})$ and $\mu_{s}$ the modified scattering coefficient $\mu_{s}=\mu_{s}^{\prime}(1-g_{1})$, where $g_{1}$ depends on the anisotropy of $f$, see (5), and often has a value close to $0.99$. In the frequency domain, i.e. assuming a time dependence of $\Phi\propto e^{-i\omega t}$ and assuming further that the scattering frequency is much larger than $\omega$, i.e. $\nu\mu_{s}\gg\omega$, one can finally get (after ignoring the dipole term of the source) the photon diffusion equation by setting the right side of Eq. (6) equal to zero, so-called DPDW equation, that we will utilize throughout this paper, namely:

[TABLE]

This equation for propagation of diffuse photon density waves with heterogeneous isotropic diffusivity, and heterogeneous absorption coefficient, in the time-harmonic regime, takes thus the final form

[TABLE]

where $\Phi$ represents the distribution of fluence field at each point ${\bf x}=(x,y,z)$ in Cartesian coordinates (for more details on cylindrical DPDW cloaks in the frequency domain, see [45]).

Moreover, let us recall that $D$ is the modified diffusivity, i.e. $D=1/(3\mu_{s})$, with $\mu_{s}$ the modified scattering coefficient and $\nu$ is the speed of light as defined in (1).

Let us now consider a map (a pull-back) from a co-ordinate system $\{x^{\prime},y^{\prime},z^{\prime}\}$ to the co-ordinate system $\{x,y,z\}$ given by the transformation characterized by $x(x^{\prime},y^{\prime},z^{\prime})$, $y(x^{\prime},y^{\prime},z^{\prime})$ and $z(x^{\prime},y^{\prime},z^{\prime})$. This change of co-ordinates is characterized by the transformation of the differentials through the Jacobian:

[TABLE]

From a geometric point of view, the matrix $\mathbf{T}\!=\!\mathbf{J}_{xx^{\prime}}^{T}\mathbf{J}_{xx^{\prime}}/\det(\mathbf{J}_{xx^{\prime}})$ allows for a description of the metric tensor inside the cloak. The only method of operation for the transformed coordinates is to replace the diffusivity (homogeneous and isotropic) and potential by equivalent ones. The effective diffusivity becomes heterogeneous and anisotropic, while the absorption coefficient and the speed of light deserve a new expression. Their properties are given by

[TABLE]

where $\mathbf{T}_{T}^{-1}$ stands for the upper diagonal part of the inverse of $\mathbf{T}$ and $T_{zz}$ is the third diagonal entry of ${\bf T}$.

The transformed equation associated with (8) reads

[TABLE]

where importantly the source term $S^{\prime}$ is a hallmark of the geometric transform. Indeed, if one places a source in the transformed medium, it seems to radiate from a shifted location, what can be considered as a mirage effect [56].

An elegant way to derive (11) is to multiply (8) by a test function that is an infinitely differentiable function which vanishes on the boundary $\partial\Omega$ of a domain $\Omega$, and to further integrate by parts over $\Omega$, which leads to 111We have used the fact that $\int_{\Omega}d{\bf x}\nabla\cdot\left((D\nabla\Phi){\phi}\right)=\int_{\Omega}d{\bf x}(D\nabla\Phi)\cdot\left(\nabla{\phi}\right)+\int_{\Omega}d{\bf x}{\phi}\nabla\cdot\left(D\nabla\Phi\right)$ and we invoked the divergence theorem in the left hand side of the equation, which esnures that this term is equal to $\int_{\partial\Omega}ds\left(D\nabla\Phi\cdot{\bf n}{\phi}\right)$, with ${\bf n}$ the unit outward normal to the boundary $\partial\Omega$.:

[TABLE]

where ${\bf n}$ is the unit outward normal to the boundary $\partial\Omega$ of the integration domain $\Omega$. Moreover, $<,>$ denotes the duality product between the space of Distributions (these are generalized functions such as Dirac, Heaviside functions) and the space of test functions (i.e. infinitely differentiable functions with compact support on $\Omega$) i.e. a pairing in which one integrates a distribution against a test function.

We now apply to (2.1) the coordinate change ${\bf x}=(x,y,z)\rightarrow{\bf x^{\prime}}=(x^{\prime},y^{\prime},z^{\prime})$ (a push-forward) and noting that

[TABLE]

where $\nabla^{\prime}$ is the gradient in the new coordinates, we end up with

[TABLE]

Upon integration by parts, and noting that ${\bf J}_{xx^{\prime}}^{-1}\nabla^{\prime}{\phi}\cdot D{\bf J}_{xx^{\prime}}^{-1}\nabla^{\prime}\Phi={(\nabla^{\prime}{\phi})}^{T}{\bf J}_{xx^{\prime}}^{-T}D{\bf J}_{xx^{\prime}}^{-1}\nabla^{\prime}\Phi$ we obtain:

[TABLE]

which is the integral form of (11). This lays the foundation of transformation light diffusion 222The boundary integral disappears between (2.1) and (2.1), which should not be a surprise since this second integration by parts is just the opposite of the one performed before coordinate change, see (2.1), and in any case one could simply consider test functions $\phi$ that vanish on the boundary of $\Omega$..

Note that in the case of cylindrical domains, when the geometric transform only affects the transverse coordinates i.e. the transformation is characterized by $x(x^{\prime},y^{\prime})$, $y(x^{\prime},y^{\prime})$, $z(z^{\prime})$ one has:

[TABLE]

where $\mathbf{T}_{T}^{-1}$ stands for the upper diagonal part of the inverse of $\mathbf{T}$ and $T_{zz}$ is the third diagonal entry of ${\bf T}$. This is left as an exercise for the zealous reader, who can otherwise find the derivation in section 4.5. on external cloaking.

We would like to also stress that we shall use pull-back transformations in the sequel for historical reasons: cloaking was first proposed for Maxwell’s equations, which behave appropriately under pull-back geometric transforms, unlike push-forward [56], and this makes possible designs of arbitrarily shaped cloaks [57, 58, 59]. Let us also take this opportunity to give proper credit to the seminal work on masking by Teixeira [60] and on anamorphism by Nicolet, Zolla and Geuzaine [61], that served as an inspiration for other research groups [62, 63, 64, 65]. In the former, the formalism of differential geometry is used to demystify cloaking and masking, while the latter shows that if an object is placed in a transformed medium, it will scatter like another object with deformed boundary, and different electromagnetic parameters, see Fig. 2, which is characterized by the pull-back transform. Nicolet et al. illustrated their theory with an object placed inside a heterogeneous anisotropic shell, which is a singular cloak. In what follows, we shall pursue an alternative route to anamorphism, which amounts to placing an object in the invisibility region of a non-singular cloak that scatters like a different object, see [39] for the case of matter-waves.

3 Transformation Optics for star domains

In this section we discuss a generalization of Pendry’s linear transformation [1]. It is applied to the general framework where it sends a domain to another domain, and is nowhere singular. The corresponding derived electromagnetic material is hence nowhere singular. Moreover, this allows the creation of mimesis whereby an object dressed with a device acquires the same electromagnetic behavior as another one of our choice. The case where the polynomial is of degree $1$ is discussed in [39].

3.1 Governing equation for heat

We consider the anisotropic diffusion equation in a two-dimensional domain without heat source

[TABLE]

where $u$ represents the distribution of temperature evolving with time $t>0$, at each point ${\bf x}=(x,y)$ in the domain. Moreover, the thermal conductivity $\underline{\underline{\kappa}}$ ($W.m^{-1}.K^{-1}$ i.e. watt per meter kelvin in SI units) has the form

[TABLE]

with $\kappa_{12}=\kappa_{21}$. Also, $\rho$ is the density ($kg.m^{-3}$ i.e. kilogram per cubic meter in SI units) and $c$ the specific heat (or thermal) capacity ($J.K^{-1}.kg^{-1}$ i.e. joule per kilogram kelvin in SI units).

3.2 Governing equations for transverse electromagnetic waves

Let us now consider the Maxwell’s equations in a cylindrical domain without source. Let us also assume that the medium is described by anisotropic permittivity $\underline{\underline{\varepsilon}}$ and permeability $\underline{\underline{\mu}}$ that have the form

[TABLE]

with $\varepsilon_{12}=\varepsilon_{21}$ and $\mu_{12}=\mu_{21}$.

In transverse electric polarisation, the Maxwell operator writes as

[TABLE]

where ${\bf H}_{l}=H_{z}(x,y,t){\bf e}_{z}$ and in transverse magnetic polarisation it writes as

[TABLE]

where ${\bf E}_{l}=E_{z}(x,y,t){\bf e}_{z}$.

We would like to show that one can recast (23) and (24) in the form of scalar wave equations. For this, we need the following result:

Property: Let ${\bf M}$ be a real symmetric matrix defined as follows

[TABLE]

Then we have

[TABLE]

Indeed, we note that

[TABLE]

Furthermore, let ${\bf M}^{\prime}$ be defined as

[TABLE]

From (3.2), we have

[TABLE]

if and only if

[TABLE]

[TABLE]

which is true if ${\bf M}^{\prime}={\bf M}_{T}^{-1}\hbox{det}({\bf M}_{T})$.

Using the above property, from (23), we derive that

[TABLE]

and from (24), we derive that

[TABLE]

We note that for isotropic permittivity and permeability, (37) and (40) reduce to the usual expressions.

If we were to solve the heat equation (12) in anisotropic media using some electromagnetic code, (37) and (40) might be useful for instance in a commercial finite element package like COMSOL.

3.3 Construction of a star cloak via non-linear transform

The non-linear transformation (44) we propose here readily extends to any star domain in $\mathbb{R}^{p},$ that is, any domain with a vantage point from which all points, including the boundary’s, are within line-of-sight. The transformation preserves every line passing through that fixed vantage point. In more details, consider three bounded star domains $D_{j}$ in $\mathbb{R}^{p},~{}p=2,3$ with piecewise smooth boundaries ${\mathbf{S}}_{j}:=\partial D_{j},$ all sharing the same chosen vantage point, $j=0,1,2,$ and further suppose they have the following configuration: $D_{2}$ contains $D_{1}$ which in turn, contains $D_{0}.$ We would like to make the domain $D_{2}$ acquire the same response to any electromagnetic waves as the different domain $D_{0}.$ More precisely, we seek to fill the region $D_{2}\setminus D_{1}$ with an electromagnetic medium whose properties (permittivity and permeability) are given by the tensor ${\mathbf{T}}^{-1}$ in Equation (54) below, so that $D_{2}$ appears to electromagnetic waves as if it were $D_{0}$ with an infinitely conducting boundary ${\mathbf{S}}_{0}:=\partial D_{0}.$

The transformation (44) is the identity outside $D_{2}$, that is, in $\mathbb{R}^{p}\setminus D_{2}$. It sends the region $D_{2}\setminus D_{0}$ to the region $D_{2}\setminus D_{1}$, in such a way that the boundary ${\mathbf{S}}_{2}:=\partial D_{2}$ of $D_{2}$ stays point-wise fixed, while ${\mathbf{S}}_{0}$ is mapped to ${\mathbf{S}}_{1}:=\partial D_{1}.$ We also set an infinitely conducting boundary condition on the inner boundary ${\mathbf{S}}_{0}$ of the region $D_{2}\setminus D_{1}.$ Any type of defect could be concealed inside $D_{1}$ and $D_{2}$ will yet still have the same electromagnetic response as $D_{0}$ with an infinitely conducting boundary ${\mathbf{S}}_{0}.$ In practice, we may divide the domains $D_{j}$ into subdomains, the part of whose boundaries lying inside ${\mathbf{S}}_{j}$ is a smooth hypersurface.

Here is a detailed description of the non-linear transformation and the subsequent construction of the tensor ${\mathbf{T}^{-1}}$. Consider a point ${\bf\underline{x}}$ of $D_{1}\setminus D_{0}$ with ${\bf\underline{x}}=(x^{1},x^{2},...)$ relative to a system of coordinates centered at the chosen vantage point 0. The line passing through x and 0 meets the boundaries ${\mathbf{S}}_{0},$ ${\mathbf{S}}_{1},$ ${\mathbf{S}}_{2}$ at the unique points ${\bf\underline{x}_{0}}=(x_{0},y_{0},...),$ ${\bf\underline{x}_{1}}=(x_{1},y_{1},...)$, ${\bf\underline{x}_{2}}=(x_{2},y_{2},...)$ respectively. Inside $D_{2}\setminus D_{0}$, the transformation reads:

[TABLE]

where $n$ is a strictly positive integer (note that the linear case whereby $n=1$ was studied in [39]).

For the practical implementation, we actually need the inverse ${\bf\underline{x}^{\prime}}\mapsto{\bf\underline{x}}$ of (44) which is the general setting that reads $x^{i}=x_{0}^{i}+\alpha_{i}~{}({(x^{\prime i})}^{n}-{(x^{i}_{1})}^{n}),~{}\rm{where}~{}\alpha_{i}=\frac{{(x_{2}^{i}-x^{i}_{0})}^{n}}{{(x_{2}^{i})}^{n}-{(x^{i}_{1})}^{n}}.$ In the 3-space with coordinates $(x^{1},x^{2},x^{3})=(x,y,z)$ we can write it as

[TABLE]

The boundary points ${\mathbf{x}_{0}}=(x_{0},y_{0},...),$ ${\mathbf{x}_{1}}=(x_{1},y_{1},...)$, ${\mathbf{x}_{2}}=(x_{2},y_{2},...)$ are functions of $x,y,\cdots$ and hence of $x^{\prime},y^{\prime},\cdots$ We will discribe them explicitly in Sections 3.4, 4.3 , for some particular cases.

The tensor ${\bf T}^{-1}$ below gives the expression of the electromagnetic properties of the transformed material inside $D_{2}\setminus D_{1}.$ Indeed, denote the permittivity and the permeability of the vacum by $\varepsilon_{0}$ and $\mu_{0}$ respectively. Then one derives the permittivity and the permeability of the new electromagnetic transformed material obtained using transformation (53), via the formula 333If we replace vaccuum by an anisotropic media, one has $\underline{\underline{\varepsilon}}^{\prime}(\mathbf{x}^{\prime})={\bf J}_{xx^{\prime}}^{-1}{\underline{\underline{\varepsilon}}(\mathbf{x})}{\bf J}_{xx^{\prime}}^{-T}\hbox{det}({\bf J}_{xx^{\prime}})$ and $\underline{\underline{\mu}}^{\prime}(\mathbf{x}^{\prime})={\bf J}_{xx^{\prime}}^{-1}{\underline{\underline{\mu}}(\mathbf{x})}{\bf J}_{xx^{\prime}}^{-T}\hbox{det}({\bf J}_{xx^{\prime}})$ instead of (49).

[TABLE]

where $\epsilon_{0}\mu_{0}=c^{-2}$, with $c$ the speed of light in vacuum.

In this article, we focus on light diffusion phenomena, but the above permittivity and permeability tensors (the square root of their product can be encapsulated within an anisotropic spatially varying refractive index) simply translate into diffusivity, scattering and absorption coefficients. Therefore and we consider the cases where the regions $D_{j}$ are cylinders over some plane curves (triangular, square, elliptic, sun flower-like cylinders, etc.) Thus, we consider the transformation mapping the region enclosed between the cylinders ${\bf S}_{0}$ and ${\bf S}_{2}$ into the space between ${\bf S}_{1}$ and ${\bf S}_{2}$ as in Fig. 3, whose inverse is

[TABLE]

The tensor ${\bf T}^{-1}$ is thus given by

[TABLE]

with

[TABLE]

where the $a_{ij}$ are

[TABLE]

and finally the partial derivatives are as follows

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Now after having derived the general formulas for mimesis, we supply the expressions of the boundary points $(x_{0},y_{0}),$ $(x_{1},y_{1})$, $(x_{2},y_{2}),$ for some particular cases of regions $D_{0},$ $D_{1},$ $D_{2}$ needed for the numerical implementations. We finally turn our attention to the numerical simulations. The explicit illustrations we have supplied to exemplify the work within this article have boundaries whose horizontal plane sections are parts of ellipses (sunflower-like petal, cross, circle) or lines (square, hexagram), see Fig. 7.

3.4 Formulas for piecewise linear boundaries

In this section, we supply the explicit expressions of the coordinates $(x_{0},y_{0}),$ $(x_{1},y_{1})$, $(x_{2},y_{2})$ of the boundary points used above, when the boundary of the star domain is piecewise linear. If a piece of the boundary of a star domain is a line of the form $y=ax+b,$ then clearly the line through the origin and a point $(x^{\prime},y^{\prime})$ intersects this piece of boundary at $(x_{i},y_{i})=(\frac{bx^{\prime}}{y^{\prime}-ax^{\prime}},\frac{by^{\prime}}{y^{\prime}-ax^{\prime}})$ and hence $\frac{\partial x_{i}}{\partial x^{\prime}}=\frac{by^{\prime}}{(y^{\prime}-ax^{\prime})^{2}}$,   $\frac{\partial x_{i}}{\partial y^{\prime}}=-\frac{bx^{\prime}}{(y^{\prime}-ax^{\prime})^{2}}$,   $\frac{\partial y_{i}}{\partial x^{\prime}}=\frac{aby^{\prime}}{(y^{\prime}-ax^{\prime})^{2}}$,   $\frac{\partial y_{i}}{\partial y^{\prime}}=-\frac{abx^{\prime}}{(y^{\prime}-ax^{\prime})^{2}}$. Of course, in the case where this piece of boundary is a vertical segment with equation $x=c,$ then the above intersection is at $(c,c\frac{y^{\prime}}{x^{\prime}}).$

3.5 Formulas for piecewise elliptic boundaries

We suppose here that a piece of at least one of our boundary curves is part of a nontrivial ellipse $\mathcal{E}_{i}$ with equation of the form $(x-a)^{2}/c_{i}^{2}+(y-b)^{2}/d_{i}^{2}$ where $(a,b)$ is our ventage point. Of course a line $(x(t),y(t))=(a,b)+t(x^{\prime}-a,y^{\prime}-b)$ passing through $(a,b)$ and a different point $(x^{\prime},y^{\prime})$ intersects $\mathcal{E}_{i}$ at two distinct points. For the construction, we need the point $(x_{i},y_{i})$ of that intersection which is nearer to $(x^{\prime},y^{\prime})$ in the sense that $x_{i}-a$ and $y_{i}-b$ have the same sign as $x^{\prime}-a$ and $y^{\prime}-b,$ respectively. We have

[TABLE]

So we have

[TABLE]

and

[TABLE]

We can then apply formulas (53)-(61) to build any generalized cloaks involving boundaries of elliptic types, where the center of the ellipse is the ventage point $(a,b)$. This will be used in Fig. 7.

4 Illustrative numerical results

In this section, we apply our theory of non-linear transforms for cloaking and anamorphism of star shaped domains to illustrative cases. We start with the simple case of carpet cloaks for diffusive light. However, before we start looking at the design of cloaks, we would like to explain how we can rigorously model a diffusion problem, which is necessarily set in unbounded domain with a bounded computational domain.

4.1 Perfectly Matched Layers through geometric transforms

To deal with the open problem, a judicious choice of coordinate transformation allows the finite element modeling of the infinite exterior domain [66]. Without loss of generality, we consider two disks $D(O,A)$ and $D(O,B)$ of center $O=(0,0)$ and radii $A$ and $B>A$ strictly including $\Omega$, we define a corona $C=D(O,B)\setminus\overline{D(O,A)}$, but the mapping can be applied to star domains. Let $(x,y)$ be a point in $\mathbb{R}^{2}\setminus\overline{D(O,A)}$ (the infinite outer domain) and $(X,Y)$ be a point in $C$, the transformation is then given by:

[TABLE]

where $R$ denotes the Euclidean norm $\sqrt{X^{2}+Y^{2}}$. This transformation may be viewed as a mapping of the finite corona $C$ with the non orthogonal coordinate system $(X,Y)$ to the infinite domain with the cartesian coordinate system $(x,y)$.

This way, the finite element discretization appears as a chained map from the reference space to the transformed space and from the transformed space to the physical space. It is also worth noting that taking Dirichlet or Neumann boundary conditions at a finite distance (without geometric transformation) would change the physics of the problem, by enforcing some values of the field at a finite distance. Importantly, since we deal with a diffusion equation which is such that the field’s amplitude is evanescent away from the source, there is no need to consider a complex valued mapping like in [67] for perfectly matched layers (PMLs), which were originally introduced by Berenger [68] to absorb electromagnetic waves propagating in all directions, without reflection. We finally point out that PMLs for elastodynamic waves require some advanced mathematical formalism of scattered field theory [69]. We implement our PMLs for diffusive light in fig. 4, with a computational domain consisting of a square of sidelength $4$m, surrounded by PMLs of width $0.5$m. The outer boundary of the computational domain has homogeneous Neumann data everywhere, except for the bottom, which has Dirichlet data.

4.2 Anamorphism with carpet cloak for virtual source

We now consider a light diffusion carpet cloak of height $2$m, and width $2$m. The inner boundary of the carpet cloak is infinite conducting (transformed Neumann boundary conditions). We report these results in Fig. 4 for an angular finite frequency of $\omega=\pi/50$ rad.s*-1*, which is appropriate for the carpet dimensions (one can easily scale down the size of carpet, and thus get high frequencies that would be in the visible spectrum as we deal with a linear problem). Some homogeneous Dirichlet boundary conditions are set on the ground plane, the inner boundary of the carpet and the infinite conducting obstacle in panels (a), (b). However, in panels (c), (d), we set the source on the inner boundary of the carpet and the object with a non-homogeneous Dirichlet data (with a normalized data of $1$), and one can see that it is virtually impossible for an external observer to deduce the shape of the source of diffusive light in (c), as it looks like it is a bare plane diffusive pseudo-wave. Let us emphasize that we used our specially designed PMLs as detailed in section 4.1.

4.3 Anamorphism with circular cloaks: Squaring the circle

In this section, we make a circle have the same (light diffusion) signature as an imaginary small square lying inside its enclosed region and sharing the same centre. In particular, its appearance to an external observer will look like that of a square. As above, the transformation will map the region enclosed between the small square and the inner circle (cloaking surface) into the circular annulus bounded by the inside and outside circles, where the sides of the square are mapped to the inner circle and the outer one stays fixed point wise. To do so, we again use the diagonals of the square to part those regions into sectors. Indeed, the diagonals provide a natural triangulation by dividing the region enclosed inside the square into four sectors. The natural continuation of such a triangulation gives the needed one. In each sector, we apply the same formulas as above, where $(x_{0},y_{0})$ are obtained from the small square as in Section 3.4 and for both $(x_{1},y_{1})$and $(x_{2},y_{2})$, we use the same formulas for ellipses in Section 3.5. For the numerical simulation, we consider a pseudo-wave angular frequency $\omega=\pi/10$ rad.s*-1* in a computational square domain of sidelength $4$m with Dirichlet data on top and bottom, and Neumann data on left and right sides, respectively. We first consider a small square of side $L_{0}=0.4$m and two circles radii $R_{1}=0.2$m, $R_{2}=0.4$m centred at $(a,b)=(0,0)$, see Fig. 5 for a result of numerical simulations for the circular cloak (a) mimicking the square (b). We next consider in Fig. 5 a circular cloak with same boundaries (c) that now mimicks a rectangle of sidelengths $L_{0}=0.4$m and $l_{0}=0.2$m (d). The tiny, yet visible, differences in the fields in Fig. 5, can be made more pronounced by considering a smaller computational domain of sidelength $2$m. In Fig. 6 (a), we therefore consider the same circular cloak surrounding the circular scatterer as in Fig. 5 (a), which we compare to the scattering by the same square scatterer as in Fig. 5 (b), see Fig. 6 (b). These two panels (a) and (b) in Fig. 6 look exactly the same, however panel (c) of Fig. 6, which considers the scattering by an infinite conducting circular obstacle of radius $0.2$m, just like the one surrounded by the cloak in panel (a), is markedly different.

4.4 Star shaped and sunflower cloaks

We now consider sunflower and star shaped light diffusion cloaks. We apply formulas (53)-(61) to build the sunflower cloaks involving boundaries of elliptic types, where the center of each ellipse is a ventage point $(a,b)$ located on a radius starting from $(0,0)$. The sunflower cloaks have 8 petals, which are deduced from a rotation through an angle $\pi/4$. See Fig. 7 (a),(b) for further details on geometry and a numerical simulation of the diffusion light scattering of such cloaks at an angular frequency $\omega=\pi/10$ rad.s*-1*. On the other hand, the design of star shaped cloaks requires an adapted triangulation of the corresponding region. That is, a triangulation that takes into account the singularities at the vertices of the boundary of the region. So that, each vertex belongs to the edge of some triangle. To the resulting triangles, one applies the same maps as in Sect. 3. See Fig. 7 (c), (d) for further details on geometry of star shaped cloak and result of numerical simulations at an angular frequency $\omega=\pi/20$ rad.s*-1*, i.e. a quasi-static case.

4.5 External cloaking for diffusive light

External cloaking is rooted in a 1994 paper [70] where a core-shell resonant system induced anomalous resonances that were later on use in the first external cloaking theory for small objects (dipoles) located in the close neighborhood of a cloak with a negatively refracting shell [71]. However, this route did not work for finite size objects and finite frequencies [72]. One possible way to extend external cloaking to finite obstacles and beyond the quasi-static limit is to make use of complemntary media, but this requires a priori knowledge of the obstacle to hide, as the cloak consists of its complementary medium [73] (let us note in passing that complementary media for optical space cancellation were originally introduced in [74]). An alternative route to external cloaking makes use of a space folding that leads to spatially varying negtaively refracting shell [75]. We propose here to adapt [75] to the case of diffusion processes. Let us consider the ’space-folding’ transform

[TABLE]

which leads to the transformed diffusivity and determinant of Jacobian

[TABLE]

where ${\bf I}$ is the $2\times 2$ identity matrix.

Let us now detail the computation of the diffusion matrix (the reader versed in transformation optics techniques can skip this part, but we find it might be useful for newcomers in this field). We start with the computation of the Jacobian matrix of the compound transformation $(x,y,z)\rightarrow(r,\theta,z)\rightarrow(r^{\prime},\theta^{\prime},z^{\prime})\rightarrow(x^{\prime},y^{\prime},z^{\prime})$ which can be expressed as:

[TABLE]

where ${\bf J}_{xr}=\partial(x,y,z)/\partial(r,\theta,z)=\partial(r\cos\theta,r\sin\theta,z)/\partial(r,\theta,z)={\bf R}(\theta)\hbox{diag}(1,r,1)$, ${\bf J}_{rr^{\prime}}=\partial(r,\theta,z)/\partial(r^{\prime},\theta^{\prime},z^{\prime})=\partial(g(r^{\prime}),\theta,z)/\partial(r^{\prime},\theta^{\prime},z^{\prime})$ and ${\bf J}_{r^{\prime}x^{\prime}}=\partial(r^{\prime},\theta^{\prime},z^{\prime})/\partial(x^{\prime},y^{\prime},z^{\prime})={\bf J}_{x^{\prime}r^{\prime}}^{-1}={[{\bf R}(\theta)\hbox{diag}(1,r^{\prime},1)]}^{-1}=\hbox{diag}(1,1/r^{\prime},1){\bf R}(-\theta^{\prime})$ since the rotation matrix ${\bf R}(\theta^{\prime})$ is such that ${\bf R}(\theta^{\prime})^{-1}={\bf R}(-\theta^{\prime})$.

We then compute the transformation matrix:

[TABLE]

where we have used the fact that the rotation matrix satisfies ${\bf R}(\theta)^{-1}={\bf R}(\theta)^{T}$ and $\theta^{\prime}=\theta$.

We now note that from (69) one has:

[TABLE]

Results of finite element computations are shown in Fig. 8 for radii $r_{c}=0.5$m, $r_{s}=0.8$m and a finite angular frequency of $\omega=\pi/10$ rad.s*-1* in panel (a) and $\pi/2$ rad.s*-1* in panel (b), where we note that we used our specially designed PMLs as detailed in section 4.1.

Such a design of external cloak for light of mass diffusion might seem far-fetched, but we note that propagation of photon density waves in random amplifying media studied in [76] with DPDW makes a theory consistent with physics of negatively refracting media [76, 77], notably versus space folding and lensing [78] that is aking to apparent negative thermal diffusivity in [79].

5 A non-linear transform for form invariant Fick’s equation of mass diffusion

Let us now assume that $\nu=1$ and $\mu_{a}=0$ in (8), this leads us to Fick’s equation, where $D$ is a chemical diffusivity in units of $m^{3}.s^{-1}$ and $\Phi$ a concentration. In [41], where the transformed Fick’s equation counterpart of (11) with $\nu=1$ and $\mu_{a}=0$ was studied in the time-domain, it was noted that the parameter $\nu^{\prime}=T_{zz}\nu=\rm{det}({\bf J})$ has no physical meaning for $\nu^{\prime}\neq 1$.

In order to overcome this pitfall, we now consider a geometric transform for cloaking which is such that det(J) is a constant. This can be done assuming that $n=2$ in (44), and further considering some boundaries parameterized by an azimuthal angle $0<\theta\leq 2\pi$:

[TABLE]

whose Jacobian is such that $\rm{det}({\bf J})(\theta)=1/\sqrt{\frac{r_{2}^{2}(\theta)-r_{1}^{2}(\theta)}{{(r_{2}(\theta)-r_{0}(\theta))}^{2}}}$.

If we further consider some radially symmetric (circular) boundaries, we can see that $\rm{det}({\bf J})$ becomes a constant. Interestingly, if we take $r_{0}=0$, we end up with the nonlinear counterpart of Pendry’s formula

[TABLE]

which has Jacobian

[TABLE]

which is such that $\rm{det}({\bf J}_{xx^{\prime}})={(1-r_{1}^{2}/r_{2}^{2})}^{-1/2}$.

Provided that $\rm{det}({\bf J}_{xx^{\prime}})$ is constant (what requires circular boundaries), it is possible to recast (11) as

[TABLE]

with $\underline{\underline{D^{\prime\prime}}}({\bf x^{\prime}})=D\rm{det}({\bf J}^{-1}){\bf T}_{T}^{-1}=D{\bf J}_{xx^{\prime}}^{-1}{\bf J}_{xx^{\prime}}^{-T}$. When $\mu_{a}=0$ and $\nu=1$, Fick’s equation is form invariant under (86).

We note that (86) also solves the issue of reduced parameters that brought a slight discrepancy in Fourier’s equation in [40], which were used for a layered thermal cloak via homogenization. Unfortunately, the form invariant Fick’s equation under (86) holds only for circular boundaries, and breaks down in the spherical case.

Let us now detail the computation of the diffusion matrix. We start with the computation of the Jacobian matrix of the compound transformation $(x,y)\rightarrow(r,\theta)\rightarrow(r^{\prime},\theta^{\prime})\rightarrow(x^{\prime},y^{\prime})$ which can be expressed as 444Noting that $r^{\prime}=\sqrt{\alpha r^{2}+r_{1}^{2}}$ hence $r=\alpha^{-1/2}\sqrt{{r^{\prime}}^{2}-r_{1}^{2}}$, we can see that (5), simplifies into ${\bf J}_{xx^{\prime}}=\alpha^{-1/2}{\bf R}(\theta)\hbox{diag}(r^{\prime}/\sqrt{{r^{\prime}}^{2}-r_{1}^{2}},\sqrt{{r^{\prime}}^{2}-r_{1}^{2}}/r^{\prime}){\bf R}(-\theta^{\prime})$.:

[TABLE]

where ${\bf J}_{xr}=\partial(x,y)/\partial(r,\theta)$, ${\bf J}_{rr^{\prime}}=\partial(r,\theta)/\partial(r^{\prime},\theta^{\prime})$ and ${\bf J}_{r^{\prime}x^{\prime}}=\partial(r^{\prime},\theta^{\prime})/\partial(x^{\prime},y^{\prime})$; Moreover, $\alpha=(r_{2}^{2}-r_{1}^{2})/r_{2}^{2}$ and we note that $\hbox{det}({\bf J}_{xx^{\prime}})=\alpha^{-1}$ as the rotation matrices ${\bf R}(\theta)$ and ${\bf R}(\theta^{\prime})$ are unimodular.

We then compute the transformation matrix:

[TABLE]

where we have used the fact that the rotation matrix satisfies ${\bf R}(\theta)^{-1}={\bf R}(\theta)^{T}={\bf R}(-\theta)$ and $\theta^{\prime}=\theta$.

We now observe that from (86) $r^{2}=(r^{\prime 2}-r_{1}^{2})/\alpha$, hence the transformed diffusivity inside the circular coating of the cloak can be expressed as

[TABLE]

where the eigenvalues of the diagonal matrix (principal values of conductivity) are from (5)

[TABLE]

with $r_{1}$ and $r_{2}$ the inner and the outer radii of the cloak.

In order to eliminate $T_{zz}=\rm{det}({\bf J})$ in the factor of $i\omega\Phi$ as in (88), we multiply the eigenvalues by $\rm{det}({\bf J})$ and find

[TABLE]

This should be compared with eigenvalues of diffusivity and determinant of Jacobian obtained from Pendry’s transform when $n=1$ 555Noting that for Pendry’s transform $r^{\prime}=\beta r+r_{1}$, with $\beta=(r_{2}-r_{1})/r_{2}$ hence $r=\beta^{-1}({r^{\prime}}-r_{1})$, we can see that ${\bf J}_{xx^{\prime}}={\bf R}(\theta)\hbox{diag}(1,r)\hbox{diag}(1/\beta,1)\hbox{diag}(1,1/r^{\prime}){\bf R}(-\theta^{\prime})={\bf R}(\theta)\hbox{diag}(1/\beta),r/r^{\prime}){\bf R}(-\theta^{\prime})$, and det(${\bf J}_{xx^{\prime}}$)$=r/(\beta r^{\prime})=(r^{\prime}-r_{1})/(\beta^{2}r^{\prime})$ and ${\bf T}^{-1}={\bf R}(-\theta^{\prime})^{-1}\hbox{diag}(\beta^{2},{(r^{\prime}/r)}^{2}){\bf R}(-\theta^{\prime})^{-T}(r/(\beta r^{\prime}))={\bf R}(-\theta^{\prime})^{-1}\hbox{diag}((r\beta)/r^{\prime},r^{\prime}/(r\beta)){\bf R}(\theta)^{-1}={\bf R}(-\theta^{\prime})^{-1}\hbox{diag}((r^{\prime}-r_{1})/r^{\prime},r^{\prime}/(r^{\prime}-r_{1})){\bf R}(\theta)^{-1}$.:

[TABLE]

Note that $(\rm{det}{\bf J})$ vanishes for $r^{\prime}=r_{1}$, in which case the frequency dependence is discarded. Physically, replacing the frequency by a time derivative, this means it should take an infinite amount of time for mass to cross the inner boundary of the cloak. This led to counterintuitive results for the transformed heat equation in [40], where it was observed that heat inside the invisibility region is constant at each time step, and similar results hold for the transformed Fick’s equation in [41]. We would like to investigate what happens inside this invisibility region with the nonlinear transform that avoids the vanishing determinant.

6 Non-linear transform for form invariant diffusion equations in time domain

In this section, we draw interesting consequences on the application of non-linear transform (86) to various diffusion equations in the time domain.

6.1 Form invariant Fick, Fourier and Schrödinger equations

Provided that $\rm{det}({\bf J})$ is constant (what requires circular boundaries), the transformed light diffusion equation in time domain writes as

[TABLE]

with $\underline{\underline{D^{\prime\prime}}}({\bf x^{\prime}})=D\rm{det}({\bf J}^{-1}){\bf T}_{T}^{-1}=D{\bf J}_{xx^{\prime}}^{-1}{\bf J}_{xx^{\prime}}^{-T}$. When $\mu_{a}=0$ and $\nu=1$, the time dependent Fick’s equation is form invariant under (86), which was not the case in [41].

The transformed Fick’s equation in [41] has served as an inspiration for the work on light diffusion cloaks by Wegener’s group, and using our geometric transform (86), it should be now possible to achieve well behaved cloaks in the transient domain, since there is no longer the issue of the heterogeneous determinant as a factor of the time derivative. Moreover, Fourier’s equation for heat also takes a nicer form, when using our transform:

[TABLE]

as now there is no more issue with vanishing values of product of density by specific heat at the inner boundary of the cloak. Therefore, the thermodynamic cloak design experimentally validated in [80] might be further improved.

Finally, one should note that there is a strong analogy between (6.2) and the time-dependent Schrödinger equation in the context of quantum mechanics. Some time-independent quantum cloaks for matter waves have been proposed by the groups of Zhang and Greenleaf using coordinate transforms [37], but they suffered from an inherent singularity due to the vanishing determinant on the inner boundary of the cloak, which is overcome with the non-linear transform using the transformed Schrödinger equation

[TABLE]

where $\hbox{det}({\bf J})$ can now be safely removed from the equation, which has not been studied thus far in the cloaking literature. Here, what plays the role of the transformed conductivity in (6.2) is $\underline{\underline{m^{\prime}}}^{-1}={\bf J}^{-1}m^{-1}({\bf x}^{\prime}){\bf J}^{-T}\hbox{det}({\bf J})$, which is a transformed (anisotropic heterogeneous) mass. Note that when $\hbox{det}({\bf J})$ is not a constant, $\hbox{det}({\bf J})V$ is a transformed (heterogeneous) potential. Moreover, $\hbox{det}({\bf J})$ in the left-hand side is a time-scaling parameter, $\hbar$ is the Planck constant and $\Psi$ is the wave function. In order to facilitate the physical interpretation of (95), one can now divide throughout by $\hbox{det}({\bf J})$, since it is a constant thanks to (86).

6.2 Homogenization route to anisotropic diffusivity

We have considerably improved the cloak’s design, since we do not have to handle the heterogeneous determinant like in [40, 41, 81]. Nonetheless, we still face the obstacle of an anisotropic diffusivity inside the cloak. This problem can be handled using techniques of homogenization see e.g. chapter 1 in [82, 83, 84] for a detailed derivation in the case of the conductivity and wave equations. The book of Milton on the theory of composites also contains many results on homogenization of layered media, that can be of interest to the reader [48]. Let us consider that the field is solution of the following DPDW equation (outside the source)

[TABLE]

inside the heterogeneous isotropic cloak $\Omega_{f}$, where $D({\bf x},{\bf y})$, $\mu_{a}({\bf x},{\bf y})$ and $\nu({\bf x},{\bf y})$ are piecewise continuous functions of period $1$ in ${\bf y}$. When $0<\eta\ll 1$, the field $\Phi_{\eta}({\bf x})$ oscillates rapidly in the structured cloak $\Omega_{f}$, with oscillations on the order of $\eta$. To filter these variations, we consider an asymptotic expansion of $\Phi_{\eta}$ solution of (6.2) in terms of a macroscopic (or slow) variable ${\bf x}=(r,\theta)$ and a microscopic (or fast) variable ${\bf x}_{\eta}=(\frac{r}{\eta},\theta)$.

The homogenization of DPDW equation can be derived by considering the ansatz

[TABLE]

where for every ${\bf x}\in\Omega_{f},\;u^{(i)}({\bf x},\cdot)$ is 1-periodic along $r$. Note that we evenly divide $\Omega_{f}$ ($R_{1}\leq r\leq R_{2}$, $0\leq\theta<2\pi$) into a large number of thin curved layers of radial length $(R_{2}-R_{1})/\eta$, but the time variable $t$ in (6.2) is not rescaled.

Rescaling the differential operator in Eq. (6.2) accordingly as $\nabla=\nabla_{\bf x}+{\bf e}_{r}\displaystyle{\frac{1}{\eta}\frac{\partial}{\partial r}}$, and collecting the terms sitting in front of the same powers of $\eta$, we obtain

[TABLE]

which is the homogenized DPDW equation where $\Phi_{0}$ is the leading order term in the ansatz (96). Here, $<f>({\bf x})=\int_{Y}f({\bf x},{\bf y})\,d{\bf y}$ is the arithmetic mean over a unit cell $Y$ along the radial axis and $\underline{\underline{D_{hom}}}$ is a homogenized rank-2 diagonal tensor, which has the physical dimensions of a homogenized anisotropic diffusivity $\underline{\underline{D_{hom}}}={\rm Diag}(D_{r},D_{\theta})$ given by

[TABLE]

We note that if the cloak consists of an alternation of two homogeneous isotropic layers of thicknesses $d_{A}$ and $d_{B}$ and diffusivities $D_{A}$, $D_{B}$, light speed $\nu_{A}$, $\nu_{B}$ and absorption coefficients $\mu_{a,A}$, $\mu_{a,B}$, we obtain the following effective parameters

[TABLE]

which are respectively described by one harmonic and three arithmetic means, where $\eta=d_{B}/d_{A}$ is the ratio of thicknesses for layers $A$ and $B$ and $d_{A}+d_{B}=1$.

6.3 Time dependence of field inside the invisibility region

We consider a multilayered cloak with radii $r_{1}=150$ mm and $r_{2}=300$mm, with $20$ layers of diffusivity varying between $0.25$ and $1680$ (in units of $10^{-5}m^{2}.s^{-1}$) and compare its diffusion against a shell with same radii and a thin inner layer of diffusivity $0.25$ (graphene could be used) and a thick outer layer of diffusivity $1680$ (still in units of $10^{-5}m^{2}.s^{-1}$). We show in Figures 9-12 the spatial distribution of mass at representative time steps (this is similar for light, heat etc., except that one has to apply a scaling factor on time due to the very diffusion speed for chemical species, light, heat etc.). We start with the case of long times, when we reach the stationary regime in Figure 9, where we note in panels (a) and (b) that the field is uniform inside the invisibility region of the cloak, unlike for the two-layer shell. We plot the profile of mass concentration in panels (b) and (d) along the x-axis highlighted by the white segment in panels (a) and (b). However, when we look at short times, we can see quite interesting antagonistic behaviours in the cloak and the shell, see Figures 11-10: field increases with time almost uniformly in the invisibility region in panels (a) and (b), whereas it varies non monotonically in panels (b) and (c), and eventually reaches a uniform gradient. Let us explain this phenomenon.

If we assume that the field (mass, or photon density or temperature) on the inner disc $D$ is zero, then the weak maximum principle ensures us that the temperature is zero everywhere inside the disc $D$. However, if the temperature on the boundary $\partial D$ of $D$ this theorem only states that $\max_{D}u=\max_{\partial D}u$. In order to demonstrate that the field inside the disc is a constant for any given time $t$, we need to apply the mean value property, which states that [85]:

Theorem (Nirenberg): Let $u\in C_{1}^{2}(D\times[0,T))$ (i.e. continuous with first and second space derivatives continuous and first time derivative continuous) solve the heat equation (94) then:

[TABLE]

for each $E_{\lambda}(x,y,t)\subset D$. Here, $E_{\lambda}(x,y,t)$ is the so-called heat ball which is a super-level set of the fundamental solution of the heat equation defined as $E_{\lambda}(x,y,t)=\{(v,w,s)\in\mathbb{R}^{2}\times[0,\infty)\mid s\leq t\;,\;\Phi(x-v,y-w,t-s)\geq\lambda\}$ and $\Phi(x,y,t)={(4t\pi)}^{-1}\exp(-(x^{2}+y^{2})/t)$.

As a corollary of this theorem, if D is connected and there exists $(x_{0},y_{0},t_{0})\in D\times[0,T)$ such that $u(x_{0},y_{0},t_{0})=\max_{\overline{D\times[0,T)}}u$, then by picking $r$ small enough so that $E_{\lambda}((x_{0},y_{0}),t_{0})\subset D\times[0,T)$, and using the mean value property, we conclude that $u$ is constant inside $E_{\lambda}(x_{0},y_{0},t_{0})$.

Next for any $(u_{0},v_{0},s_{0})\in D\times[0,T)$ such that the line segment connecting $(x_{0},y_{0})$ to $(v_{0},w_{0})$ is in $D$, we can show that $u(v_{0},w_{0},s_{0})=u(x_{0},y_{0},t_{0})$ whenever $s_{0}<t_{0}$ by covering the line segment connecting $(v_{0},w_{0},s_{0})$ and $(x_{0},y_{0},t_{0})$ with the heat balls.

Finally, since $D$ is connected, any $(v_{0},w_{0})$ can be connected from $(x_{0},y_{0})$ via finitely many line segments. And therefore we have $u(v,w,s)=u(x_{0},y_{0},t_{0})$ for all $(v,w)\in D$, $s<t_{0}$.

We have therefore shown that if the invisibility region $D$ is connected, which is clearly the case for a disc, but also for the invisibility region shown in Figures 9-12, and if there exists $(x_{0},y_{0})\in D$ such that $u(x_{0},y_{0})=\max_{\overline{D\times[0,T)}}u$, then we are ensured that $u\equiv u(x_{0},y_{0})$ in $D\times[0,T)$. Bearing in mind that the value of the field on the boundary of the invisibility region $D$ in Figures 9-12 is the solution of the transformed equation (6.2), we note that it is an appromation of the field solution of the diffusion equation (8) at the point $r=0$ in (86) mapped onto $r^{\prime}\in\partial D$. We therefore deduce that the field in the circular invisibility region in Figures 9-12 is a constant field equal to the value of the field at the center point of the disc $D$ before geometric transform. In the present case, the value of $u$ solution of (8) at the center point is $T/2$, hence after the geometric transform the value of $u$ solution of (6.2) on $\partial D$ equals $T/2$, and so is the value of $u_{0}$ solution of its approximate homogenized equation (6.2). Applying the above corollary of the mean value theorem, we deduce that $u\equiv T/2$ in the disc $D$. Importantly, when we consider a different point as the origin of the geometric transform, the maximum value inside $D$ can be smaller or larger than $T/2$, but it cannot exceed $T$. Indeed, if we assume that the value of the field at the point taken as the origin $r=0$ of the geometric transform is $T$, maximum of the field solution of (8), then the above discussion leads to $u\equiv T$ inside the invisibility region.

7 Conclusion and perspectives

In this review article, we have proposed some models of generalized cloaks that create some optical illusion for light in diffusive media: infinite conducting obstacles dressed with these cloaks scatter light like other infinite conducting obstacles. In these cloaks, an electric wire could in fact be hiding a larger object near it. In fact, any object could mimic the signature of another.

For instance, we design a cylindrical cloak so that a circular obstacle behaves like a square obstacle, thereby bringing about one of the oldest enigma of ancient time : squaring the circle! The ordinary cloaks then come as a particular case, whereby objects appear as an infinitely small infinite conducting object (of vanishing scattering cross-section) and hence become invisible. Such generalized cloaks are described by nonsingular permittivity and permeability, even at the cloak’s inner surface. We have proposed and discussed some interesting applications. The results within this paper can be extended to cloaking and anamorphism for electromagnetic, acoustic, hydrodynamic and elastodynamic waves propagating in anisotropic heterogeneous media.

Regarding the fabrication of the light diffusion cloak, there is of course the route work [42] but we note the possibility to use instead a self-assembly technique for colloidal nanoparticles like in [87]. For the mass diffusion cloak, it has been shown in [88] that certain membranes of graphene oxyde allow for control of drug diffusion, we thus believe GO could be used as a building block of a biocloak (see Fig. 13 for an artistic view) that would also consist of an alternation of peptides as proposed in [41]. Indeed, a cloak whose innermost layer is made of graphene would meet the requirement of high anisotropy of transformed diffusivity in (90).

Let us conclude this review by some perspectives on diffusion media mimicking celestial mechanics [89, 90, 91, 92, 93, 94], conformal mapping for diffusion processes [2] and some governing equations in other fields of physics.

7.1 Black hole analogues for heat, light and other diffusion phenomena

Concept of black hole was first discussed by Laplace in the framework of the Newton mechanics, whose event horizon radius is the same as the Schwarzschild’s solution of the Einstein’s vacuum field equations. If all those objects having such an event horizon radius but different gravitational fields are called as black holes, then one can simulate certain properties of the black holes using electromagnetic fields and metamaterials due to the similar propagation behaviours of electromagnetic waves in curved space and in inhomogeneous metamaterials. In a recent theoretical work by Narimanov and Kildishev [89], with an independent proposal by Cheng and Cui [90], an optical black hole has been proposed based on metamaterials, in which the theoretical analysis and numerical simulations showed that all electromagnetic waves hitting it are trapped and absorbed. The Schwarzschild metric of a black hole in cylindrical coordinates $(r,\theta)$ is

[TABLE]

and in spherical coordinates $(r,\theta,\Theta)$ it is

[TABLE]

where $L$ is the event horizon of the black hole.

If we are interested in control of heat, transformation thermodynamics [40], we get that a medium with anisotropic heterogeneous conductivity and heterogeneous denisty and heat capacity

[TABLE]

where $r=\sqrt{x^{2}+y^{2}}$, should be the analogue of a black hole for heat.

We obtain the following transformed conductivity and product of density and heat capacity in cylindrical coordinates

[TABLE]

and in spherical coordinates

[TABLE]

where $L$ is the event horizon.

It is of course possible to translate this in the language of transformation light diffusion, or mass diffusion, and this is left as an exercise. We further note that other designs of optical black holes [92, 93, 94] and wormholes [95, 96] can be applied to diffusion process.

7.2 Conformal mappings for light, heat, mass and other diffusion physics

Inspired by the seminal paper of Leonhardt on conformal optics [2], we consider a homogeneous isotropic medium within which the following governing equation for a diffusion process holds:

[TABLE]

for $u(x,y)$ in an infinite transverse plane $-\infty<x,y<\infty$. Note that for instance (8) can be encapsulated in (103), simply by taking $\Omega=D^{-1}(i\omega\nu^{-1}-\mu_{a})$. In practice, the metamaterial will have a finite extension and the diffusion direction of pseudo waves may be slightly tilted without causing an appreciable difference to the ideal case (note that unlike for [2], $\Omega$ lies in the complex plane, not on the positive real line).

To describe the spatial coordinates in the propagation plane, complex numbers $z=x+iy$ are used with the partial derivatives $\partial_{x}=\partial_{z}+\partial_{z}^{\star}$ and $\partial_{y}=i\partial_{z}-i\partial_{z}^{\star}$, where the asterisk symbolizes complex conjugation. In this way, the Laplace operator takes the form $\partial^{2}_{x}+\partial^{2}_{y}=4\partial^{\star}_{z}\partial_{z}$ so that (103) is expressed in complex coordinates as follows

[TABLE]

In the case of gradually varying diffusivity, that is $\Omega:=\Omega(z)$, one needs to assume that $\Omega$ must not vary by much over the scale of a diffusion length.

Suppose we introduce new coordinates $w$ described by an analytic function $w(z)$ that does not depend on $z^{\star}$. Such functions define conformal maps [97] that preserve the angles between the coordinate lines. Because

[TABLE]

we obtain in $w$ space a diffusion equation with the transformed index profile $\Omega^{\prime}(w)$ that is related to the original one as

[TABLE]

Suppose that the diffusive medium is designed such that $\Omega(z)$ is the modulus of an analytic function $g(z)$. The integral of $g(z)$ defines a map $w(z)$ to new coordinates where, according to (105), the transformed index $\Omega^{\prime}$ is unity. Consequently, in $w$ coordinates, the diffusion is indistinguishable from empty space where diffusion flux follows straight lines. The medium performs a diffusion conformal mapping to homogeneous space. If $w(z)$ approaches $z$ for $w\to\infty$, all incident pseudo-waves appear at infinity as if they have diffused through empty space, regardless of what has happened in the medium. However, as a consequence of the Riemann Mapping Theorem [97], nontrivial $w$ coordinates occupy Riemann sheets with several $\infty$, one on each sheet.

Following [2], we consider the simple map

[TABLE]

that is realized by the pseudo refractive index profile $\Omega=\mid 1-a^{2}/z^{2}\mid$. The constant $a$ characterizes the spatial extension of the medium. The function (107) maps the exterior of a circle of radius $a$ on the $z$ plane onto a given Riemann sheet and the interior onto another. Pseudo-waves traveling on the exterior $w$ sheet may well be passing through the branch cut between the two branch points $\pm 2a$. In continuing their diffusion, the pseudo-waves approach $\infty$ on the interior $w$ sheet. From the point of view of an observer on the physical $z$ plane, they cross the circle of radius $a$ and approach the singularity of the pseudo refractive index at the origin. For general $w(z)$, only one $\infty$ on the Riemann structure in $w$ space corresponds to the true $\infty$ of physical $z$ space and the others to singularities of $w(z)$. Instead of traversing space, pseudo-waves may cross the branch cut to another Riemann sheet where they approach $\infty$. For an observer in physical space, the rays seem to be attracted towards some singularities of the pseudo refractive index. Instead of becoming invisible, the medium casts a shadow that is as wide as the apparent size of the branch cut. Nevertheless, the diffusion on Riemann sheets could serve as a mathematical tool for developing the design of metamaterials that shape the flow of diffusion processes.

We illustrate numerically this conformal cloak for diffusion processes in Fig. 14. Results of finite element computations are shown in Fig. 14 for a conformal cloak given by (107) for $a=0.25$m and a finite angular frequency of $\omega=\pi/10$ rad.s*-1* in panel (a) where we note that we used our specially designed PMLs as detailed in section 4.1, but here the background medium is spatially varying, and so are PMLs. In panel (b), the object illuminated by the same pseudo wave in a homogeneous diffusive medium acts as a concentrator, the streamlines are pinched inside it (so the flux therein is large).

7.3 Non-Eudlidean cloaking for diffusion processes

In the sequel, we make use of the correspondence between the anisotropic diffusion equation for heat (12) and the anisotropic equations (40) and (40) for transverse electromagnetic waves. In this paper, we focus on a 2D version of the non-Euclidean invisibility cloak that has been described in the previous sections in flat space. We perform the earlier described geometric transformation from 2D virtual space to the plane $xy$ of physical space, and then, we add an extra dimension $z$, assuming that the material parameters only depend upon $xy$ coordinates, thereby generating a cylindrical 3D physical space. We consider transverse electromagnetic waves, which split in $s$ and $p$ polarizations. For calculating the line elements in this 3D physical space, we then have to add the term $dz^{2}$ to the line elements corresponding to the 2D situation. To describe the Euclidean part of the cloak, we use bipolar cylindrical coordinates $(\tau,\sigma,z)$ in physical space that are related to the Cartesian $(x,y,z)$ as

[TABLE]

Here, the parameter $a$ characterizes the size of the cloak. Analogous relations hold in virtual space where $x$, $y$, and $\sigma$ are replaced by $x^{\prime}$, $y^{\prime}$, and $\sigma^{\prime}$, respectively. The mapping between the two spaces is given in the bipolar cylindrical coordinates by the following relation between $\sigma$ and $\sigma^{\prime}$:

[TABLE]

and

[TABLE]

and the coordinates $\tau$, $z$ coincide in the two spaces (the mapping only involves inplane coordinates). The line element $ds^{2}=d{x^{\prime}}^{2}+d{y^{\prime}}^{2}+d{z^{\prime}}^{2}$ in virtual space is equal to

[TABLE]

The tensors of permittivity and permeability for this region have been calculated by Tyc et al. in [101] using the same general method as in [100], see also the supplemental material [102]. A short, but insightful, review of [100] has been written by experts in generalized cloaking [103] and can serve as an entry door in this field. In [101], the following expressions can be found:

[TABLE]

which can be translated in terms of anisotropic diffusivity and product of density and heat capacity in the anisotropic governing equation for heat, see Eq. (12):

[TABLE]

The mapping between the non-Euclidean part of virtual space and the corresponding part of physical space requires more work. Here, we will not reproduce all the corresponding formulas of [102], but give some important steps, as was done in [101]. One starts with a sphere of radius $r=(4/\pi)a$ (see Eq. (S43) in [102]) parameterized by spherical coordinates— the latitudinal angle $\theta$ and the longitudinal angle that is denoted by $\sigma^{\prime}$. This is not yet the sphere of virtual space, but is related to it by a suitable Möbius transformation (see Eq. (S40) in [102]) represented on the sphere via stereographic projection (see Eq. (S37)-(S39) in [102]). Then, the coordinates $(\sigma^{\prime},\theta)$ are mapped to the bipolar coordinates $(\sigma,\tau)$ of physical space as follows. The mapping between $\sigma^{\prime}$ and $\sigma$ is given by (110), where $\sigma^{\prime}$ is set to one of the intervals $[-2\pi,-\pi]$ or $[\pi,2\pi]$ by adding or subtracting $2\pi$. The corresponding intervals of $\sigma$ are then $[-\pi,-3\pi/4]$ and $[3\pi/4,\pi]$. The mapping between $\theta$ and $\tau$ can be derived by combining Eqs. (S41) and (S46) in [102]. After some calculations, this yields

[TABLE]

The line element of the non-Euclidean part of virtual space is

[TABLE]

This line element is slightly different from the line element in Eq. (S47) of [102], as noted in [101], where details can be found.

Tyc et al. simplify Eq. (115) in a more amenable form in [101] that leads to the following entries for the tensors of permittivity and permeability

[TABLE]

Regarding the heat equation, this translates as follows in terms of anisotropic diffusivity and product of density by heat capacity

[TABLE]

One can easily adapt these parameters to the DPDW and mass diffusion equations and this is left as a simple exercise.

Interestingly, Smerlak derived the Fokker-Plank equation in the general curved space-time framework back in 2012 [99]. Bearing this in mind, we would like to end this review article by some perspectives on non-Eudlidean cloaking for diffusion processes in curved spaces.

7.4 Other transformation based governing equations

There is a vast literature on diffusion equations, that span chemical, medical, physical and engineering sciences, and that also appear in economics and image processing. There is therefore room for further investigations of transformation based diffusion equations applied to societal problems where anisotropy might play a more prominent role. Since we have a particular interest in mass diffusion applied to biocloaks, we first consider a diffusion model with a mathematical setup by Haidar [104] for a neutron wave propagation akin to heat pseudo-waves, motivated by modulated neutron source experiments with applications in Boron and Gadolinium capture for cancer treatment [105]. Let $\phi$ denote the neutron flux distribution in a Boron and/or Gadolinium loaded cancerous tissue of thin slab domain $\Omega$, the diffusion equation of a flux of speed $\nu$ thermal neutrons in $\Omega$ is similar to that of diffusive light (7), more precisely

[TABLE]

where $\phi({\bf x},t)=\nu N({\bf x},t)$, with $N$ the neutron density. Moreover, $D$ is the diffusion length of these neutrons in $\Omega$, $\mu_{a}$ is the macroscopic absorption cross section of these neutrons and $D/\mu_{a}=L^{2}$, where $L$ is the diffusion length of these neutrons. In light water, $D=0.164$ cm, $\mu_{a}=0.0198$ cm*-1*, and $L=2.88$ cm. The life time of thermal neutrons is $T={(\nu\mu_{a})}^{-1}=213$ $\mu$s in $H_{2}O$. The source $S_{0}$ is periodic, and boundary conditions on $\partial\Omega$ can be zero flux $\phi({\bf x},t)\mid\partial\Omega=0$ (Dirichlet), modulated neutron flux intensity $[D({\bf x})\nabla\Phi({\bf x},t)]\cdot{\bf n}=\xi S_{0}$ (Neumann) where $\xi$ is a coupling parameter between $\Omega$ and an adjacent domain, which is filled for instance with non cancerous cell.

Another equation of interest is the Fokker-Planck equation, that describes the time evolution of the probability density function of the velocity of a particle under the influence of drag forces and random forces, as in Brownian motion. For an Ito process driven by the standard Wiener process ${\bf W}_{t}$ and described by the stochastic differential equation $d{\bf X}_{t}=\boldsymbol{\mu}({\bf X}_{t},t)dt+\boldsymbol{\sigma}({\bf X}_{t},t)d{\bf W}_{t}$ where ${\bf X}_{t}$ and ${\displaystyle{\boldsymbol{\mu}}(\mathbf{X}_{t},t)}$ are N-dimensional random vectors, ${\displaystyle{\boldsymbol{\sigma}}(\mathbf{X}_{t},t)}{\boldsymbol{\sigma}}(\mathbf{X}_{t},t)$ is an $N{\displaystyle\times}M$ matrix and ${\displaystyle\mathbf{W}_{t}}$ is a $M$-dimensional standard Wiener process, the probability density ${\displaystyle p(\mathbf{x},t)}$ of the random variable ${\displaystyle\mathbf{X}_{t}}$ satisfies the Fokker-Planck equation [98]

[TABLE]

with drift vector ${\displaystyle{\boldsymbol{\mu}}=(\mu_{1},\ldots,\mu_{N})}$ and diffusion tensor $D_{ij}({\bf x},t)=\sum_{k=1}^{M}\sigma_{ik}({\bf x},t)\sigma_{jk}({\bf x},t)$.

Under coordinate change ${\bf x}\longmapsto{\bf x}^{\prime}$, this equation takes the form

[TABLE]

with ${\bf J}=\partial(x_{1},...,x_{N})/\partial(x^{\prime}_{1},...,x^{\prime}_{N})$ the Jacobian, ${\bf J}^{-1}$ its inverse and ${\bf J}^{T}$ its transpose.

Other PDEs of interest for coordinate changes include the Boltzmann transport equation [106, 107, 108], which thanks to Mouhot and Villani [109] is know to have well-behaved solutions (if a system obeying the Boltzmann equation is perturbed, then it will return to equilibrium, rather than diverging to infinity). In a simplified form, it can be written similarly to the Vlasov-Poisson nonlinear equation (also known as Landau equation)

[TABLE]

where the unknown function, $\Phi$ is the density of electrons in a gas in the phase space of all positions and velocities of particles (${\bf x},{\bf v}$), $m$ is the mass of an electron and the mean field electrostatic force ${\bf F}=-e{\bf E}=-e\nabla\Delta^{-1}(4\pi\rho)$, with $e=1.6021766208\times 10^{-19}$ Coulomb and $\rho({\bf x},t)$ the density of charges. All quantities have $N\geq 2$ indices.

Under coordinate change ${\bf x}\longmapsto{\bf x}^{\prime}$, assuming that ${\bf v^{\prime}}({\bf x}^{\prime})={\bf v}({\bf x}^{\prime})$, the Vlasov-Poisson nonlinear equation (120) takes the form

[TABLE]

with ${\bf J}=\partial(x_{1},...,x_{N})/\partial(x^{\prime}_{1},...,x^{\prime}_{N})$.

Finally, we would like to mention that geometric transforms might find also interesting applications in mathematical models for market. For instance, the Black-Scholes equation [110, 111]. Let $s_{i}(t)\in[0,+\infty)$, $i=1,2,...,N$ denote the value of the $i^{th}$ underlying asset at time $t>0$ and $u({\bf s},t)$ denote the price of an option, where ${\bf s}=(s_{1},s_{2},...,s_{N})$. Then the option price follows the following generalized Black–Scholes partial differential equation

[TABLE]

with final condition $u({\bf s},T)=\Lambda({\bf s})$ (payoff function at maturity $T$), where $r$ is the constant riskless interest rate, $\sigma_{i}$ are volatilities of $s_{i}$, and $\rho_{ij}$ are the asset correlations between $s_{i}$ and $s_{j}$.

Under coordinate change ${\bf s}\longmapsto{\bf s}^{\prime}$, this equation takes the form

[TABLE]

with ${\bf J}=\partial(s_{1},...,s_{N})/\partial(s^{\prime}_{1},...,s^{\prime}_{N})$ the Jacobian, ${\bf J}^{-1}$ its inverse and ${\bf J}^{T}$ its transpose. This transformed Black-Scholes equation might find some applications in control of options trading. The transformed Black-Scholes model could help hedge the option by buying and selling the underlying asset in just the right way and, as a consequence, to eliminate risk. This could serve as the basis of hedging strategies such as those engaged in by investment banks and hedge funds, by adding anisotropy in the existing models.

We would like to finalise this review article by some perspectives on image processing. Back in 2009, some of us have applied a simple algorithm for comparative analysis of different datasets related to heparan sulphate oligosaccharides [112]. In this work, it was observed that when a cloud of points displays some preferential direction (anisotropy) for repeated experiments, the quantification of their reproducibility can be made through evaluation of an average Lyapunov exponent characterizing the area-preserving nature of a sequence of effective ellipses [112]. We believe that the evolution of the data sets with time could be encapsulated in a diffusion equation with anisotropic coefficients. It is actually well-known that diffusion of heat smoothes the temperature function $u$, and the noting that $u(x,t)$ is equivalently found by minimizing problems of $L^{2}$ norm of gradients using the identity $\partial u/\partial t=\nabla\cdot(\nabla u/{\|u\|})$, some denoising algorithm can be applied to images thanks to anisotropic diffusion equation [113, 114, 115]. We believe that image processing could benefit from geometric transforms so as to improve smoothing algorithms.

Acknowledgements

A.D. and S.G. acknowledge funding from ERC grant ANAMORPHISM. S.G. T.P. acknowledge funding from AMIDEX foundation of Aix-Marseille University through project BIOCLOAK. M.F. acknowledges funding from the Qatar National Research Fund (QNRF) through a National Priorities Research Program (NPRP) Exceptional grant, NPRP X-107-1- 027.