Analytical and Numerical Analysis of Linear and Nonlinear Properties of an rf-SQUID Based Metasurface

TL;DR

This paper develops a comprehensive model for rf-SQUID based metasurfaces interacting with electromagnetic waves, deriving analytical solutions in the linear regime and numerically exploring nonlinear effects like bistability.

Contribution

It introduces a coupled differential equation model for rf-SQUID metasurfaces, providing analytical solutions and numerical methods to study nonlinear electromagnetic interactions.

Findings

Analytical expressions for reflection, transmission, and absorption in the linear regime.

Numerical demonstration of bistability in nonlinear response.

Convergence and error analysis of the numerical scheme.

Abstract

We derive a model to describe the interaction of an rf-SQUID (radio frequency superconducting quantum interference device) based metasurface with free space electromagnetic waves. The electromagnetic fields are described on the base of Maxwell's equations. For the rf-SQUID metasurface we rely on an equivalent circuit model. After a detailed derivation, we show that the problem that is described by a system of coupled differential equations is wellposed and, therefore, has a unique solution. In the small amplitude limit, we provide analytical expressions for reflection, transmission, and absorption depending on the frequency. To investigate the nonlinear regime, we numerically solve the system of coupled differential equations using a finite element scheme with transparent boundary conditions and the Crank-Nicolson method. We also provide a rigorous error analysis that shows convergence of the scheme at the expected rates. The simulation results for the adiabatic increase of either the field's amplitude or its frequency show that the metasurface's response in the nonlinear interaction regime exhibits bistable behavior both in transmission and reflection.

Equations

$$\vec{\nabla }\times \vec{E}(\vec{r},t)$$

$$\vec{\nabla }\times \vec{H}(\vec{r},t)$$

$$\vec{\nabla }\cdot \vec{D}(\vec{r},t)$$

$$\vec{\nabla }\cdot \vec{B}(\vec{r},t)$$

$$\vec{D}(\vec{r},t)$$

$$\vec{B}(\vec{r},t)$$

$${\partial }_t^2\vec{H}(\vec{r},t)+c^2\vec{\nabla }\times \vec{\nabla }\times \vec{H}(\vec{r},t)=c^2\vec{\nabla }\times \vec{j}(\vec{r},t),$$

$${\vec{H}}^{\Sigma }(z)$$

$$+{\vec{H}}^{+}(z)\Theta (z-a),$$

$${\vec{H}}^{-}(-a)={\vec{H}}^0(-a),{\vec{H}}^0(a)={\vec{H}}^{+}(a).$$

$${\partial }_z{H}_y^{\Sigma }$$

$$+\left(H_y^{+}-H_y^0\right)\delta (z-a)$$

$$={\left({\partial }_z{H}_y\right)}^{\Sigma },$$

$${\partial }_z{\left({\partial }_z{H}_y\right)}^{\Sigma }$$

$$+\left({\partial }_z{H}_y^{+}-{\partial }_z{H}_y^0\right)\delta (z-a).$$

$$\vec{\nabla }\times \vec{\nabla }\times {\vec{H}}^{\Sigma }$$

$$=-{\partial }_z^2{H}_y{\hat{e}}_y+\left({\partial }_z{H}_y^{-}-{\partial }_z{H}_y^{+}\right)\delta (z){\hat{e}}_y.$$

$$\vec{j}(\vec{r},t)=\frac{I(t)}{2bA_{\rm c}}\left(\begin{array}[]{c}z\\ 0\\ -x\end{array}\right)\Theta(R_{a}-\rho)\Theta(\rho-R_{i}),$$

$${\left(\vec{\nabla }\times \vec{j}(\vec{r},t)\right)}_x={\left(\vec{\nabla }\times \vec{j}(\vec{r},t)\right)}_z=0.$$

$$\vec{\nabla }\times \vec{j}(\vec{r},t)$$

$$+\frac{I(t)}{2b{A}_{\mathrm{c}}}\rho \left[\delta (\rho -R_i)-\delta (R_a-\rho )\right]{\hat{e}}_y.$$

$$\vec{\nabla }\times \vec{j}(\vec{r},t)$$

$$+\frac{I(t)}{2b{A}_{\mathrm{c}}}\rho \left[\delta (b+\rho )-\delta (b-\rho )\right]{\hat{e}}_y.$$

$$\underset{b\to 0}{\lim }\delta (b+\rho )=\underset{b\to 0}{\lim }\delta (b-\rho )$$

$$\underset{b\to 0}{\lim }\frac{\Theta (b+\rho )\Theta (b-\rho )}{2b}$$

$$\underset{b\to 0}{\lim }\vec{\nabla }\times \vec{j}(\vec{r},t)=\frac{2I(t)}{A_{\mathrm{c}}}\delta (\rho ){\hat{e}}_y.$$

$${\partial }_t^2\vec{H}-c^2{\partial }_z^2\vec{H}=\left({\partial }_z{\vec{H}}^{+}-{\partial }_z{\vec{H}}^{-}+\frac{2I}{A_{\mathrm{c}}}\vec{n}\right)c^2\delta (z),$$

$$\left({\partial }_t^2-c^2{\partial }_z^2\right)\vec{H}(\vec{r},t)=0.$$

$${\partial }_z{\vec{H}}^{+}(0,t)-{\partial }_z{\vec{H}}^{-}(0,t)=-\frac{2I(t)}{A_{\mathrm{c}}}\vec{n},$$

$$I_{\mathrm{ring}}(t)$$

$$=C{\partial }_{t}U(t)+\frac{U(t)}{R}+I_{\mathrm{JJ}}(t).$$

$$U_{\mathrm{JJ}}(t)=-{\partial }_t{\Phi }_{\mathrm{ring}}(t).$$

$$\phi (t)=\frac{1}{\hslash }{\int }_{\mathrm{JJ}}\vec{p}(t)\cdot d\vec{l}=\frac{1}{\hslash }{\int }_2^1\vec{p}(t)\cdot d\vec{l}$$

$$I_{\mathrm{JJ}}(t)$$

$$U_{\mathrm{JJ}}(t)$$

$$I_{\mathrm{ring}}(t)$$

$$\frac{\hslash }{2e}{\partial }_t\phi (t)$$

$$\vec{p}(\vec{r},t)=\hslash \vec{\nabla }\varphi (\vec{r},t)=2m_e{\vec{v}}_{\mathrm{q}}(\vec{r},t)+2e\vec{A}(\vec{r},t),$$

Full text

Introduction

In the last years, researchers spent tremendous efforts in understanding and developing electrodynamic metamaterials that operate at different frequencies from the GHz range up to the visible Anlage (2011); Jung et al. (2014); Maas et al. (2013); Cao et al. (2014); Enkrich et al. (2005). Metamaterials consist of unit cells that are mostly periodically arranged in space. These artificially structured materials are designed to offer control on the propagation of electromagnetic fields inaccessible with natural materials Turpin et al. (2014). For that, one relies frequently on tiny structures inside a host medium to form the unit cells: the metaatoms. Metaatoms shall assure a strong interaction of the electromagnetic field with matter. Therefore, resonances are often exploited. Moreover, controlling the scattering properties of the individual metaatom is key to tailor the emerging material properties. For a long time, a magnetic response had been looked after but many more properties can be tailored. The metaatoms themselves can be described by purely classical means, e.g., within the context of electrodynamics itself if they are made from ordinary materials such as dielectrics or metals Kurter et al. (2011), but also by quantum mechanical means if required. That would hold when the metaatom consists of, e.g., a flux qubit as an artifical two-level system P. Macha et al. (2014).

A referential example for a metaatom with a strong magnetic response is the split ring resonator (SRR) Singh et al. (2009); Niesler and Peter (2012); Ricci and Anlage (2006); Ricci et al. (2007). An SRR is a metal ring acting as an inductance with a small gap forming a capacitance, i.e., an LC-circuit. In a natural way, determined by its geometry and material, the SRR has a resonance frequency. However, the downside of using resonant structures made from ordinary metals is (a) a spurious intrinsic absorption that lowers the quality factors and with that the achievable dispersion in the effective properties of the actual metamaterial and (b) their limitation to a fixed resonance frequency upon fabrication Anlage (2011).

Both aspects can be mitigated while relying on superconducting materials in the design of metaatoms. First of all, superconductors do not suffer from dissipation Ricci et al. (2005) because they carry current that is not subject to Ohmic resistance due to the bosonic character of their charge carriers Mangin and Kahn (2017). That requires, however, an operational frequency corresponding to an energy that is smaller than the binding energy of the Cooper pairs. This restricts the use of superconductor based metamaterials to the GHz or at most the lower THz frequency range. But superconductors also solve the second aforementioned problem as their properties sensitively depend on the environmental temperature Chen et al. (2010) and magnetic fields they are exposed to Jung et al. (2013); Trepanier et al. (2013, 2017); Ricci et al. (2005). Thus, external parameters have an impact on the intrinsic resonance properties of the metaatoms.

A further option to tune metamaterials is by exploiting nonlinear effects in the interaction of the electromagnetic wave with the metamaterial. A well understood nonlinear element in the field of superconductivity is the Josephson junction (JJ) Stewart (1968). It introduces both nonlinearity into the system and makes use of the low-loss properties of superconducting charge transport. In 2007, it was proposed to put Josephson junctions into the gap of an SRR made of superconducting material and to use these devices as metaatoms Lazarides et al. (2007). Such structures are called rf-SQUID ring resonators (radio frequency Superconducting QUantum Interference Device). They are already well investigated in the context of transmission line theory Butz et al. (2013); Jung (2014); Jung et al. (2013). Additionally, rf-SQUID rings provide a tunable intrinsic inductivity via an externally applied magnetic field Jung et al. (2013). Hence, an rf-SQUID is a natural and promising candidate as a building block to create novel, efficient, and tunable metamaterials.

Besides metamaterials as volumetric matter, is has been appreciated that comparable control over electromagnetic fields can be offered by metasurfaces, i.e., thin films made from a monolayer of metaatoms. Then, it is not refraction and diffraction in the bulk media that shall be controlled but rather reflection and transmission from an array Dolling et al. (2006); Odit et al. (2016); V. Shalaev (2007); Walther et al. (2012). In that context, a question of scientific importance concerns the proper description of rf-SQUID based metasurfaces and the exploration of their linear and nonlinear properties. The present contribution develops a comprehensive theoretical framework for that purpose and explores linear and nonlinear properties.

We start by developing an interaction model of free space electromagnetic waves with a thin film loaded with rf-SQUIDs. In the spirit of the metasurface, we assume the spatial extent of the metaatom $L_{\rm MA}$ to be much smaller than the operational wavelength $\lambda$ of the incident waveOHara et al. (2008), i.e., $L_{\rm MA}\ll\lambda$. Also the thickness $d$ of the metasurface is considered to be much smaller than $\lambda$, i.e., $d\ll\lambda$. Hence, the interaction region can be regarded as infinitesimally thin in the propagation direction of the waves (Caputo et al., 2012). On the one hand, we will describe the dynamics of the system by the continuity of the magnetic field and a jump discontinuity of its first derivative with respect to space, derived from Maxwell’s equations Pfeiffer and Grbic (2013); Selvanayagam and Eleftheriades (2013). On the other hand, we will use circuit theory and macroscopic quantum effects to describe the inner dynamics of the current and voltage drop inside the rf-SQUID. From these considerations we derive in Section I two coupled differential equations that describe (a) the propagation of the incident field coupled to the rf-SQUIDs and (b) the temporal evolution of the internal dynamics of the rf-SQUID metasurface. The wellposedness of our system of equations is proven in Section II. We take this as a justification for the reliability of our approach. The optical response of the metasurface in the linear regime is discussed analytically in Section III. Selected properties of the optical response of the metasurface in the nonlinear regime are discussed numerically in Section IV. For these simulations, we outline an efficient scheme and discuss details of the spatial and temporal discretization of the governing differential equations. This discussion also contains a rigorous error analysis showing error estimates for both discretizations. Finally, we conclude on our work in Section V.

I Derivation of the model

The derivation of a model that describes the interaction of an rf-SQUID ring film with an electromagnetic wave will take into account Maxwell’s theory of electrodynamics and circuit theory to express the dynamics in the rings. For the latter we rely on the resistively and capacitively shunted junction model (RCSJ model) of the Josephson junction (JJ) in the rf-SQUID ring Kleiner et al. (2004); Caputo et al. (2012, 2015). Moreover, we use models of macroscopic quantum effects, such as the Josephson effects and flux quantization. These different aspects are documented in the following subsections. The final purpose of this section is to derive a set of coupled differential equations [cf. (31)] that describe in a self-consistent manner the evolution of the electromagnetic field and the internal dynamics in the rf-SQUID ring film.

I.1 Electrodynamics - Maxwell’s equations

To describe the interaction of an rf-SQUID ring film with electromagnetic fields, we start with Maxwell’s equations describing the evolution of electromagnetic fields in time Jackson (1999),

[TABLE]

We set the polarization and the magnetization of the film’s host material to zero, since for simplicity, we consider the film to be located in vacuum, such that

[TABLE]

Differentiating (1a) with respect to time and applying the $\operatorname{curl}$ operator to (1b) together with (2) yields

[TABLE]

where $c$ is the speed of light in vacuum. This is the governing wave equation that we have to solve to express the dynamics of the electromagnetic field.

As illustrated in Fig. 1, we assume that the film comprising the rf-SQUID rings has a thickness of $d=2a$. Without loss of generality it is located around $z=0$ inside the $x$-$y$-plane, such that $z\in[-a,a]$. This thickness shall be much smaller than the wavelength of the incident light, i.e., $d\ll\lambda$. The orientation of the rings can be arbitrary but we bias our description towards the assumption that the strongest interaction is observed when the rings are upright in the film and the normal vector of the rf-SQUID rings points in $y$-direction. We consider normally incident light which renders our model to be translationally invariant in $x$-$y$-direction, thus $\vec{H}(\vec{r},t)=\vec{H}(z,t)$. Moreover, we assume linear polarization for the magnetic field in $y$-direction. This assures a strong coupling of the magnetic field to the ring at their preferential orientation.

We start with the evaluation of the left-hand side of (3) and have a look at the double $\operatorname{curl}$ of the linearly polarized magnetic field $\vec{H}(\vec{r},t)=H_{y}(z,t)\hat{e}_{y}$. It needs a special treatment since the magnetic field $\vec{H}$ is not differentiable twice with respect to space. We make a piecewise ansatz in the three different regions of space (to the left, to the right, and inside the film) and introduce the notation

[TABLE]

where $\vec{H}$ is forced to be continuous, i.e.,

[TABLE]

We compute

[TABLE]

using the chain rule $\partial_{z}\Theta(f(z))=\delta(f(z))\partial_{z}f(z)$ and (5). Following the same arguments as before, we arrive at

[TABLE]

Note that the differences in the brackets do not vanish in general. However, performing the limit $a\rightarrow 0$, we can simplify the expression further and get

[TABLE]

For the evaluation of the right-hand side of (3), consider a current density created by a current flowing within a superconducting metal ring. We parametrize the current density in the plane $y=0$ that fully contains the enclosed area of the ring. We call the ring’s cross sectional area $A_{\rm c}=\pi(R_{\rm a}-R_{\rm i})^{2}/4$, where $R_{i}=R-b$ and $R_{a}=R+b$ are the inner and outer radius of the ring, respectively. Therefore, one can parametrize the current density’s motion as

[TABLE]

where $\rho=\sqrt{x^{2}+z^{2}}$ and $2b$ is the ring’s thickness. Due to the cylindrical symmetry, at the position of the origin, where the center of the ring is placed and the interaction with the electromagnetic field takes place, we see that

[TABLE]

The $\operatorname{curl}$ of the current density in (9) is therefore given by

[TABLE]

After expressing $R_{a}$ and $R_{i}$ through $R$ and shrinking the ring $(R\rightarrow 0)$, such that it is contained inside the thin film, (11) reads

[TABLE]

In the limit of a vanishing thickness of the ring, such that $b\rightarrow 0$, we notice, that

[TABLE]

Therefore, only the fist term in (12) remains and we are left with

[TABLE]

The Dirac distribution in (14) confines the curl of the current density to $z=0$. Since the model for the entire film is assumed to be translationally invariant in $x$-direction, we omit the confinement to $x=0$ here. This accounts for the existence of other rings along the $x$-direction inside the film. Please note, that the limit in (14) is performed in such a way, that the current density inside the ring remains constant. For an arbitrarily oriented ring’s area with unit normal vector $\vec{n}$, from (3), (8) and (14), we can summarize

[TABLE]

linking the current’s motion in the ring and the hereby radiated magnetic field $\vec{H}$. Please note, that for $z\neq 0$, we deal with a free wave equation for the magnetic field in the negative and positive half space, respectively,

[TABLE]

Additionally, from (15), we obtain a jump condition for the first spatial derivative of the magnetic field at the position of the film $z=0$,

[TABLE]

illustrating that $\vec{H}$ is not differentiable twice at $z=0$. The latter two equations are the most important equations describing the evolution of the field coupled to a film that carries a current. In the next subsection we elaborate on the details of the rf-SQUID ring to express the current in the film that is driven by an external field.

I.2 Circuit theory - Kirchhoff’s rules

We apply the RCSJ model of a Josephson junction to the rf-SQUID ring Kleiner et al. (2004); Clarke and Braginski (2004). It states, that a JJ can be replaced in circuit diagrams by the junction itself (JJ), a shunted capacitor (C), and a shunted resistor (R). Additionally, the ring’s loop is taken into account as an inductance connected in series, see Fig. 2(b). Kirchhoff’s nodal rule yields

[TABLE]

Kirchhoff’s mesh rule yields $\sum_{\rm ring}U(t)=-\partial_{t}\Phi_{\rm ring}(t)$, where $\Phi_{\rm ring}(t)$ denotes the flux penetrating the ring’s enclosed surface acting as the electromotive force. As the voltage drop across all of the three shunted elements is the same, we pick the voltage drop across the JJ for convenience and write

[TABLE]

I.3 Macroscopic quantum effects

To link the current through a Josephson junction and the voltage drop across it to the phase difference $\varphi$, we use the gauge invariant definition

[TABLE]

of the phase difference of the superconducting wave functions on either side of the JJ Josephson (1962). The integration path $2\rightarrow 1$ in (20) refers to Fig. 2(a). For $\varphi(t)$, Josephson’s equations

[TABLE]

hold. Hence, from (18) and (19), we obtain

[TABLE]

Yet another macroscopic quantum effect has to be taken into account, namely the flux quantization in a superconducting loop. This effect occurs when considering a superconducting bulk material device containing a hole Doll and Näbauer (1961). We state the general expression of the momentum of a Cooper pair of mass $m=2m_{e}$ and charge $q=2e$ inside a superconductor Mangin and Kahn (2017)

[TABLE]

where $\vec{v}_{\rm q}$ denotes the velocity of the Cooper pairs and $\vec{A}$ is the magnetic vector potential, that obeys $\vec{B}=\vec{\nabla}\times\vec{A}$. The term in the middle of (23) is generated when applying the momentum operator to the general expression of a condensate’s wave function in real space, i.e., $\Psi(\vec{r},t)=\sqrt{n}e^{i\phi(\vec{r},t)}$ with position-independent density distribution $n$. We integrate (23) along a closed loop around the superconducting ring

[TABLE]

The integration path is chosen such that the distance from the path to the surface of the ring is everywhere larger than the London penetration depth $\lambda_{\rm L}$ of the electromagnetic field into the ring. Then, the integration path only coincides with a current carrying region inside the JJ and $\vec{v}_{q}(\vec{r},t)\cdot d\vec{l}=0$ holds elsewhere inside the ring. Thus, the current has to be integrated only across the JJ. Please note, that even if the magnetic field penetrates deep into the superconductor and there is no such current-free region, we only have to correct the geometric inductance $L$ in (27) by a kinetic term, i.e., $L\rightarrow L^{\prime}=L+L_{\rm kinetic}$. By Stokes’ theorem, we obtain

[TABLE]

where $\Phi_{0}=h/2e$ is the flux quantum, $A_{\rm r}=\pi R_{\rm i}^{2}$ is the enclosed area of the ring, and $m$ an integer. Furthermore, $\Phi(t)$ is the externally applied magnetic flux via a magnetic field and in the first term of the right hand side we used (20). We can see that the total flux $\Phi_{\rm ring}(t)$, which penetrates the ring, consists of the externally applied flux $\Phi(t)=\Phi_{\rm ext}(t)$ and an additional term that describes a screening current $I_{\rm ring}$ in the ring. Its role is to force the enclosed flux onto an integer multiple value of the flux quantum $\Phi_{0}$, i.e.,

[TABLE]

where the current is determined by (21a). Equation (26) can be reformulated to

[TABLE]

In the case of the interaction of an electromagnetic wave with the ring, the external flux penetrating the ring’s enclosed surface is provided by the magnetic component of the wave, i.e.,

[TABLE]

For the sake of brevity and comprehensible readability, we drop the spatial and temporal dependencies from now on, whenever the situation is unambiguous. Hence, from (22a), using (22b), (27), and (28), we arrive at

[TABLE]

Equation (29) is a nonlinear oscillator $\varphi(t)$, that is driven by the magnetic field vector $\vec{H}(z=0,t)$ at the position of the film. Equation (17) describes the back-action of the current in the film on the magnetic field via the jump condition of its first spatial derivative. We know that accelerated charges send out radiation, such that the current can be regarded as the source of the electromagnetic field $\vec{H}(\vec{r},t)$. Equations (17) and (29) constitute the central equations of the interaction model.

I.4 Normalization of the model

To investigate their mathematical structure, we boil (17) and (29) down to dimensionless equations by introducing

[TABLE]

We use both dimensionless time $\tau:=\tilde{\omega}t$ and space $\xi:=z/\tilde{\lambda}$ variables, where $\tilde{\omega}$ defines a characteristic time scale of the oscillator and $\tilde{\lambda}$ a characteristic length scale of the system. Inserting the above relations, (17) and (29) transform to the dimensionless expressions

[TABLE]

For the sake of simplicity, we will now rename the spatial and temporal coordinates back to the original ones and write $\tau\rightarrow t$ and $\xi\rightarrow z$, both still being dimensionless, i.e.,

[TABLE]

After the physical model of the dynamics in the system has been derived, we now discuss the wellposedness of the problem from a mathematical point of view. This offers a clear indication that the derived system of equations is reasonable.

II Wellposedness

We now show that (31) together with initial conditions

[TABLE]

has a unique solution ${\vec{h}}:[0,T]\to H^{1}(\mathbb{R})^{3}\cap H^{2}(\mathbb{R}\!\setminus\!\{0\})^{3}$ and $\varphi:[0,T]\to\mathbb{R}$, where $H^{k}(\mathbb{R})$ denotes the Sobolev space of order $k\in\mathbb{N}$.

We prove wellposedness of (31) with (32) using Ref. Leibold, 2017. To keep notation short, we introduce the spaces

[TABLE]

equipped with the respective standard norms. Using the short notation $u=(\vec{u}_{h},u_{\varphi})=({\vec{h}},\varphi)$, $u_{0}=({\vec{h}}_{0},\varphi_{0})$ and $u_{\texttt{t},0}=(\vec{h}_{\texttt{t},0},\varphi_{\texttt{t},0})$, we derive the weak form

[TABLE]

where $m:X\times X\to\mathbb{R},a,b:V\times V\to\mathbb{R}$ and $f:V\to X$ are defined by

[TABLE]

Since $\kappa>0$, $m$ is an inner product for $X$. Moreover $b$ is positive semidefinite and continuous with

[TABLE]

Furthermore, by Gauss’s theorem $a$ is continuous. It is also symmetric and satisfies a Garding inequality, i.e.,

[TABLE]

for all $c_{G}>\max\{1,\kappa^{-1}\}$. Finally the right-hand side $f$ is Lipschitz continuous with constant $\kappa\beta$, i.e.,

[TABLE]

Therefore, Theorem 3.3 in Ref. Leibold, 2017 yields the existence of a unique solution $u\in C^{2}(0,T;X)\cap C^{1}(0,T;V)$ of (31) with (32) for initial values u_{0}\in\bigl{(}H^{2}(\mathbb{R}^{3}\!\setminus\!S)\times\mathbb{R}\bigr{)}\cap V satisfying the jump condition (31b) and $u_{\texttt{t},0}\in V$ (even if $u_{0}$ is a bistable point of the potential).

We want to emphasize that an analogous proof holds true without the assumption of the film being translationally invariant with respect to the $x$-$y$-plane, i.e., ${\vec{h}}$ and $\varphi$ being also functions of $x$ and $y$. This situation is relevant if the metasurface shall show a position dependent response to encode further functionalities.

III Analytical treatment in linear approximation

First, we consider a special case of (31) and calculate reflection and transmission from the film in the linear regime, where we assume that the effective interaction potential can be accurately described by a parabola. As the one-dimensional case suffices to describe the reflection and transmission coefficients of the film, we continue to assume that the film is placed in the $x$-$y$-plane. We further assume the electromagnetic field and the film to be translationally invariant with respect to the $x$-$y$-plane and calculate the reflection and transmission coefficient of an incident plane wave that is $y$-polarized with its magnetic field.

We make a small-amplitude ansatz to linearize the differential equations

[TABLE]

where $\vec{h}_{s}$ and $\varphi_{s}$ are the static components of the magnetic field and the phase difference. For $\varphi_{s}=0$, we obtain

[TABLE]

where we used $\omega_{0}=\sqrt{1+\beta}$. Analogously to the calculations in Ref. Caputo et al., 2015, we proceed with a time-harmonic ansatz to investigate the model in frequency space

[TABLE]

Equation (35) yields

[TABLE]

We can plug $\phi(\omega)$ from (37) into the jump condition (31b) and as a result we eliminate one equation, i.e.,

[TABLE]

where we defined

[TABLE]

III.1 Reflection and transmission coefficients

For the spatial dependence of the magnetic field, we also make a harmonic ansatz in the two half-spaces $z<0$ and $z>0$ and impose continuity of the magnetic field at the position of the film $z=0$, see Fig. 3, i.e.,

[TABLE]

Note that due to the normalization of the model and the propagation direction of the electromagnetic wave along the $z$-axis, it holds $\omega(\vec{k})=k_{\rm z}$. From (38) and using the ansatz in (40), we find

[TABLE]

Using (41) as well, we obtain

[TABLE]

Our goal is to express both the reflected wave $\vec{f}_{\rm ref}$ and the transmitted wave $\vec{f}_{\rm trans}$ through the incoming wave $\vec{f}_{\rm in}$ only. On that account, we project (43b) onto $\vec{n}$ and obtain

[TABLE]

As desired, using (44), we can write (43a) and (43b) as functions of the incoming field amplitude $\vec{f}_{\rm in}$ only. We obtain

[TABLE]

For the following, we make assumptions concerning the geometry of the problem. Assuming, that $\theta$ is the inclination angle of the ring’s normal vector $\vec{n}$ with respect to the incoming field amplitude $\vec{f}_{\rm in}$, i.e., $\vec{f}_{\rm in}\cdot\vec{n}=|\vec{f}_{\rm in}|\cos\theta$, we find the reflection coefficient $R(\omega,\theta)$ and the transmission coefficient $T(\omega,\theta)$ according to

[TABLE]

where we defined the ”common denominator” as

[TABLE]

Since the absorption function $A(\omega,\theta)$ has to fulfill the energy conservation relation $A+R+T=1$, we find by comparison

[TABLE]

For $\theta=0$, Fig. 4 shows the reflection and transmission coefficients, as well as the absorption function for the two cases $\beta=0$ (a) and $\beta=1$ (b). Both cases show the resonance at the resonance frequency $\omega_{0}=\sqrt{1+\beta}$ of the rings. However, for a notably large current through the Josephson junction $(\beta>0)$, in the DC limit the film acts as a reflector. This can be explained taking into account self-induction of the rf-SQUID ring loop for small frequencies.

IV Simulations in the Nonlinear Interaction Regime

Up to this point, the effects have been computed analytically in the linear interaction regime after having performed linear approximations in (34). We next consider nonlinear effects by numerical simulations. As soon as the amplitude of the incoming magnetic field $\vec{h}_{\rm inc}$ exceeds a critical value, the trigonometric expressions in the equations of motion of our model can no longer be replaced by the linear term of their Taylor expansion.

IV.1 Wellposedness on a bounded domain

In order to introduce a spatial discretization, we restrict the computational domain to a bounded subdomain $\Omega_{\ell}\coloneqq[\ell_{1},\ell_{2}]\subset\mathbb{R}$ for some $\ell_{1}<0<\ell_{2}$. This leads to the following simplified model problem, where the magnetic field ${\vec{h}}$ and the phase $\varphi$ satisfy

[TABLE]

The wellposedness of the reduced system is shown analogously to the approach for the general setting. First we introduce the spaces

[TABLE]

equipped with the standard norms. The bilinear forms $m:X_{1}\times X_{1}\to\mathbb{R},a:V_{1}\times V_{1}\to\mathbb{R}$ and the right-hand side $f:V_{1}\to X_{1}$ are defined as before, but with $\Omega_{\ell}$ instead of $\mathbb{R}$. Only the bilinear form $b:V_{1}\times V_{1}\to\mathbb{R}$ changes significantly:

[TABLE]

but as all bilinear forms and the right-hand side have the same properties as before, Theorem 3.3 in Ref. Leibold, 2017 again yields the existence of a unique solution $u\in C^{2}(0,T;X_{1})\cap C^{1}(0,T;V_{1})$ of (50).

IV.2 Space discretization

We discretize in space using finite elements on a grid $\widehat{\Omega}_{\ell}$ of $\Omega_{\ell}$. In order to resolve the jump condition (50c) correctly, we require [math] to be a grid point. We denote the maximal length of the intervals in $\widehat{\Omega}_{\ell}$ by $h$. We further introduce the space $\mathcal{P}_{k}(\widehat{\Omega}_{\ell})$, consisting of piecewise polynomials of degree at most $k\in\mathbb{N}$ in every interval in $\widehat{\Omega}_{\ell}$, and the space $\widehat{V}\coloneqq\mathcal{P}_{k}(\widehat{\Omega}_{\ell})^{3}\times\mathbb{R}$.

[TABLE]

where the initial values $\widehat{u}_{0}$ and $\widehat{u}_{\texttt{t},0}$ discrete versions of their continuous counterparts. The discrete bilinear forms $\widehat{m},\widehat{a}:\widehat{V}\times\widehat{V}\to\mathbb{R}$ are approximations of $m$ and $a$, where the integrals are replaced by a quadrature rule of order at least $k^{2}$. Therefore, the discrete bilinear forms coincide with their continuous counterparts on $\widehat{V}\times\widehat{V}$. Hence, they satisfy the same assumptions and we get from Theorem 3.6 in Ref. Leibold, 2017 the following semi-discrete error estimate.

Theorem (semi-discrete error estimate): For the exact solution $u=(\vec{u}_{h},u_{\varphi})^{T}\in C^{2}(0,T;W^{k,2}(\Omega_{\ell})\times\mathbb{R})\cap C^{1}(0,T;W^{k+1,2}(\Omega_{\ell})\times\mathbb{R})$ of the continuous problem and the discrete solution $\widehat{u}\in C^{2}(0,T;\widehat{V})$ of (52), the following estimate holds for all $t\in[0,T]$

[TABLE]

IV.3 Full discretization

We use the Crank-Nicolson scheme for the time discretization of (52). First, we define $\widehat{A},\widehat{B}:\widehat{V}\to\widehat{V}$ via

[TABLE]

for all $\widehat{u},\widehat{v}\in\widehat{V}$. The Crank-Nicolson scheme with time step $\tau$ is then given by

[TABLE]

From Corollary 3.7 in Ref. Leibold, 2017 we get the following result.

Theorem (fully discrete error estimate): For the exact solution $u=(\vec{u}_{h},u_{\varphi})^{T}\in C^{4}(0,T;X_{1})\cap C^{3}(0,T;V_{1})$ and the numerical approximation $\widehat{u}^{n}\in\widehat{V}$ obtained by the Crank-Nicolson scheme (54), the following error estimate holds for all $t_{n}=n\tau\leq T$

[TABLE]

IV.4 Validation

In Fig. 5(a), the error estimate (53) is numerically confirmed for $\alpha=0,\beta=3,\kappa=1$ and $k=1$. For the initial values, we chose a Gaussian-modulated sinusoidal pulse of the form

[TABLE]

for $z\in[-\frac{1}{2},\frac{1}{32}]=\widehat{\Omega}_{\ell}$. Since the exact solution is unknown, we computed a reference solution on a finer grid. As predicted, we see linear convergence in the spatial resolution for the error measured in the energy norm.

Since it is not possible to have any inflow with the boundary conditions (50d), the support of the initial values yields a lower bound for the size of the spatial domain. However, for the reflection and transmission coefficients it is sufficient to know the magnetic field at any pair of points $\varepsilon_{1}<0<\varepsilon_{2}$ arbitrarily close to the film, i.e., $-\varepsilon_{1},\varepsilon_{2}\ll 1$. Therefore, we adapt the boundary condition at $\ell_{1}$ in order to allow an incident wave $\vec{h}_{\textnormal{in}}:[0,T]\to\mathbb{R}$ entering the computational domain $\Omega_{\ell}$ from the left side. Using d’Alembert’s formula, $\vec{h}_{\textnormal{in}}$ is uniquely defined by the initial values ${\vec{h}}_{0}$ and $\vec{h}_{\texttt{t},0}$. So in the following, we replace (50d) by the new boundary conditions

[TABLE]

Although not covered by our analysis, numerical experiments also show first order convergence, as can be seen in Fig. 5(b).

The benefit of these boundary conditions is the possibility to drastically reduce the computational domain. In fact, the choice $-\ell_{1}=\ell_{2}=h$ means that the grid $\widehat{\Omega}_{\ell}$ contains only the two intervals $(-h,0)$ and $(0,h)$. Therefore, the numerical effort for the spatial discretization is completely independent of the spatial resolution.

As the Crank-Nicolson scheme is unconditionally stable, one can even keep the number of time steps constant. Therefore, there is no dependency between spatial resolution and the computational effort. So the maximal computable spatial resolution is solely restricted by the machine epsilon, as the condition number of the resolvent $\begin{pmatrix}I&-\frac{\tau}{2}I\\ \frac{\tau}{2}\widehat{A}&I+\frac{\tau}{2}\widehat{B}\end{pmatrix}$ is growing proportionally to $h^{-2}$.

IV.5 Simulation results

To investigate the nonlinear effects, we increase the amplitude of the incoming light’s magnetic field component and observe the excited phase difference $\varphi$ as well as the reflected and transmitted field amplitude. We apply two qualitatively different types of sources in the simulation setup.

One option is to sweep the amplitude $h_{\rm inc}$ of the incoming magnetic field at a fixed driving frequency $\omega(t)=\omega_{\rm D}$. The corresponding magnetic field for $t\in[0,T]$ has the form

[TABLE]

where $h_{\rm inc,min}$ is still part of the linear interaction regime, but $h_{\rm inc,max}$ is not. 2. 2.

Another way to observe nonlinear effects is to sweep the frequency of the incoming light at fixed amplitude. In the linear interaction regime, the system provides its maximum amplitude response of $\varphi$ at resonance frequency $\omega_{\rm D}=\omega_{0}$. This is not necessarily the case, when we go to the nonlinear interaction regime. When we perform a frequency sweep of the incoming light at fixed amplitude, the incoming magnetic field will be of the form

[TABLE]

We will choose the setting in such a way, that at some time, the incoming field is in resonance with the structure, i.e., $\omega_{\rm min}<\omega_{0}<\omega_{\rm max}$ holds.

Figure 6(a) shows simulation results of the first kind, applying a source term according to (57) to an rf-SQUID with parameters $\alpha=0.1$ and $\beta=1.5$. One can see the amplitude of the stationary state oscillation of $\varphi$, belonging to different incident plane wave amplitudes. The blue triangles pointing to the right indicate dynamic parameter sweep simulation results from small amplitudes upwards towards larger ones. Vice versa, the red triangles pointing to the left indicate dynamic sweep simulation results from large amplitudes downwards towards smaller ones. One can observe, that in a certain range of amplitudes the curves do not coincide. The hysteresis loop that occurs for $\varphi$ can also be observed in the reflected and transmitted field amplitudes, see Fig. 6(b)-(c). In a certain bistable region, the amplitudes of the reflected and transmitted waves, respectively, do not only depend on the amplitude of the incident plane wave, but also on the direction this amplitude value has been approached from.

We further apply a plane wave source according to (59), i.e., we keep the amplitude of the incident plane wave fixed throughout the entire simulation and sweep its driving frequency $\omega_{\rm D}$ over a frequency range, which contains the resonance frequency $\omega_{0}=\sqrt{1+\beta}$. The results are shown in Fig. 7(b)-(d). The resonance frequency is indicated by the vertical black dashed lines. In Fig. 7(a), $h_{\rm inc}=10^{-3}$ is in the linear interaction regime. The other three figures are plots of a simulation done at $h_{\rm inc}=1.2$, when nonlinear effects already play a role. The damping parameter is chosen to be $\alpha=0.1$ and the SQUID parameter is $\beta=2.5$. The sweeps were done first from $\omega_{\rm min}=1.3$ to $\omega_{\rm max}=2.1$ in an increasing way (blue triangles pointing to the right), afterwards vice versa in a decreasing way (red triangles pointing to the left). The response of the system is different in either case, exhibiting the manifestations of the nonlinear terms in the equations of motion. We can observe that the resonance of both the phase difference $\varphi$ and the reflected wave $h_{\rm ref}$ are shifted to smaller frequencies. Compare this observation to the Duffing oscillator, which takes into account the cubic term in the $\sin$-expansion as well Sharma et al. (2012); Wu et al. (2014). Hence, the minimum of the transmitted wave amplitude occurs at smaller frequencies than in the linear case as well. Thus, one can tune the effective resonance frequency of the entire system by increasing the amplitude of the incoming magnetic field amplitude $h_{\rm inc}$ to smaller frequencies.

V Discussion of realistic parameters

In order to retain a physical intuition of the dimensions of the involved quantities, we briefly insert realistic values into the parameters of the normalized system introduced in section I.4. From Ref. Butz et al., 2013 and Ref. Butz, 2015, we take the values presented in the experimental study of a transmission line based rf-SQUID interaction setup. We plug in the geometric inductance $L=83\,\rm pH$, an intrinsic capacitance of the JJ of $C=0.02\,\rm pF$, a critical current of the JJ of $I_{\rm cr}=1.8\,\rm\upmu A$, and a resistance $R=1600\,\rm\Omega$. This results in the characteristic frequency $\tilde{\omega}\approx 2\pi\cdot 124\,\rm GHz$ and a scaling parameter for the spatial coordinates of $\tilde{\lambda}\approx 0.4\,\rm mm$. Due to the large resistance, the damping parameter is $\alpha\approx 0.04$ and the rf-SQUID in this case is nonhysteretic with $\beta\approx 0.45$. We want to point out here, that in this contribution we used a higher value $\alpha=0.1$ (corresponding to a resistance of $R=640\,\rm\Omega$) to investigate the effect of dissipative losses, which otherwise wouldn’t have shown its impact on the dynamics of the system. All other parameters remain unchanged by this replacement of the resistance. We would also like to emphasize that the authors of Ref. Butz et al., 2013 and Ref. Butz, 2015 used an additional shunted capacitance of $C=2\,\rm pF$ to lower the resonance frequency of the rf-SQUID by roughly one order of magnitude to the range of $\omega\approx 2\pi\cdot 10\,\rm GHz$. In this study, we generically used $\kappa=1$. We obtain this value by assuming $R_{\rm i}=10\,\rm\upmu m$, according to Fig. 1, and a cross sectional area of $A_{\rm c}=1600\,\upmu m^{2}$. This corresponds to a torus-shaped ring with quite small aspect ratio. The operational frequency is in the range of $f\approx 100\,\rm GHz$. Consequently, the wavelength $\lambda$ is in the order of $\rm mm$. The above choice of geometry is therefore justified and we meet the condition $d\ll\lambda$, since the wavelength is two orders of magnitude larger than the radius of the rings. The scaling parameter for the magnetic field evaluates to $\vec{h}=1.2\,\frac{\rm m}{\rm A}\cdot\vec{H}$. This means, that a value of $h=1$ corresponds to a magnetic induction of $B=1.05\,\rm\upmu T$.

VI Conclusion

We have derived an interaction model of an rf-SQUID loaded infinitesimally thin film with electromagnetic waves. In a strictly mathematical treatment, we showed that our problem is wellposed. Therefore, a unique solution of the system of coupled differential equations exists.

We have treated the model in the linear small-amplitude interaction regime analytically. In this limit, we derived analytical expressions for the reflection and transmission coefficients of the film as well as for the absorption function.

To investigate nonlinear effects, we proposed a numerical scheme based on the finite element method and the Crank-Nicolson scheme. We further showed rigorous error estimates and presented a numerical scheme based on transparent boundary conditions with inflow, where the computational effort is independent of the spatial resolution. With these methods, we simulated the dynamics in the system numerically and found bistable and hysteretic behavior in the nonlinear interaction regime.

In future work, interaction between the rf-SQUIDs inside the film has to be taken into account. It has been proposed already to assume an interaction via mutual inductance between the rings Hizanidis et al. (2016); Lazarides et al. (2015); Lazarides and Tsironis (2013). Moreover, one has to investigate, if the coupling of the electric component of the wave to the rf-SQUID is relevant in the description of the interaction.

Acknowledgements

The authors gratefully acknowledge financial support by the Deutsche Forschungsgemeinschaft (DFG) through CRC 1173. We also wish to thank Jan Brehm and Alexey Ustinov for many helpful comments.