%---------------------------------------------
% modelisation cryo aimants supra 
% 
%dimanche 29 mars 2020, 19:19:30 (UTC+0200)
%
%---------------------------------------------
\documentclass{article}
\usepackage{epsfig,graphicx}
\usepackage{float}
\usepackage{amsmath}
%,dsfont}
\usepackage{authblk}
\renewcommand{\thepage}{ }

\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\cotan}{\mathrm{cotan}}

\newtheorem{theorem}{Theorem}[section]
\newtheorem{definition}[theorem]{Definition}
\newtheorem{proposition}[theorem]{Proposition}
\parskip .3cm


\title{\bf Modeling cryo-magnets}

\author[1]{J. G. Caputo \thanks{caputo@insa-rouen.fr}}
\author[2]{I. Danaila \thanks{ionut.danaila@univ-rouen.fr}}
\author[2]{C. Tain \thanks{cyril.tain@univ-rouen.fr}}


\affil[1]{Laboratoire de Math\'ematiques, INSA de Rouen Normandie\\ 76801 Saint-Etienne du Rouvray, France.}
\affil[2]{Laboratoire Raphael Salem, Universit\'e de Rouen et CNRS.}


\renewcommand\Authands{ and }
\date{ }

\begin{document}
\maketitle




\date{\ }

\begin{abstract}
Cryomagnets can be made by applying a pulse of magnetic field
to cooled bulk superconductors, resulting in a "trapped field" 
that can reach several teslas. These easily transportable devices
could be used in many applications, in particular magnetic levitation
and motors.  \\
We analyze the theoretical models proposed; these are based on
Maxwell's equations and a constitutive equation which is usually
Bean-Kim. In particular, we show that they cannot yield amplitude
dependant effects and this is contrary to the experimental observations
of "giant flux jumps". A solution to this is to couple directly 
Maxwell's equations to the Ginzburg-Landau free energy and derive a
Maxwell-Ginzburg-Landau Lagrangian for the vector potential
and the order parameter. We discuss this new model and present
1D and 2D simulations.

\end{abstract}

\maketitle

\end{document}

\tableofcontents


\section{Fait et a faire}


\begin{itemize}

\item  Analyse equations constitutives London et Bean {\bf FAIT}
\item  modele dynamique Maxwell-Ginzburg-Landau 1D {\bf FAIT}
\item 2D Maxwell-London code {\bf FAIT} 

\item compute fixed points  , stability {\bf A FAIRE}
\item adapt FDTD 1D code to Maxwell-Ginzburg-Landau 1D {\bf A FAIRE}

\item 2D code Maxwell-Ginzburg-Landau {\bf A FAIRE}

\end{itemize}



\section{The model}

Consider Maxwell's equations \cite{feynman}
\begin{eqnarray}
\label{max1} \nabla \cdot \mathbf{E} = {\rho \over \epsilon_0} ,\\
\label{max2} \nabla \times \mathbf{E} = -\mathbf{B}_t ,\\
\label{max3} \nabla \cdot \mathbf{B}=0,\\
\label{max4} c^2 \nabla \times \mathbf{B}= {\mathbf{J} \over \epsilon_0} + \mathbf{E}_t
\end{eqnarray}
where $\mathbf{E},\mathbf{B}$ are respectively the electric and magnetic fields,
$\rho$ is a static charge density, the indices represent time derivative, $\mathbf{J}$ is a current (moving charge)
density and $c$ is the light velocity. 


We model a superconducting material. Then no
static charges are present so that we can assume $\rho=0$.
We introduce the vector potential $\mathbf{A}$, such that
\be\label{vpot} \mathbf{B} = \nabla \times \mathbf{A}, \ee
Equation (\ref{max3}) is satisfied and equation (\ref{max2}) implies
\be\label{vpott} \mathbf{E} = - \mathbf{A}_t .\ee
We also use the London gauge
\be\label{lgauge} \nabla \mathbf{A}=0.\ee
Then, the last Maxwell equation (\ref{max4}) can be written as
\be\label{At}
\mathbf{A}_{tt} -c^2 \Delta \mathbf{A} ={\mathbf{J} \over \epsilon_0} .  \ee


The main problem now is to give an expression for $\mathbf{J}$.
A first choice is the London approximation, see \cite{vanduzer}
\be\label{J1}
\mathbf{J} = {c^2 \over \lambda^2} \mathbf{A} , \ee
where $\lambda$ is the London characteristic length.
The second expression is the Bean-Kim model \cite{bean}
\be\label{J2}
\mathbf{J} =  \mathbf{J_0} {B \over B+B_0}  \tanh({E \over E_c}),\ee
where $E,B$ are the respective magnitudes of  $\mathbf{E},\mathbf{B}$.
Equations (\ref{At}) and (\ref{J1}) or (\ref{J2}) are our main equations.
We will label this system Maxwell-super.



\section{Existence of a non trivial fixed point}

The experiments on magnetic flux jumps, see for example \cite{weinstein} 
report that there is a treshold in magnetic field above which the
sample switches from a transient to a permanent magnetization.
Fig. \ref{wein} shows the trapped field as a function of the radius
for different applied fields $B$ from 1.44 T to 2.69 T
Notice the jump observed for $B > 1.57 T$. 
\begin{figure}[H]
\centerline{ \epsfig{file=figs/wein.eps,height=4 cm,width=12 cm,angle=0}}
\caption{ Trapped field observed as a function of the radius for
different applied fields $B$ from 1.44 T to 2.69 T. }
\label{wein}
\end{figure}
This indicates that
the system switches from one type of state close to zero to another.
Dynamically, this means that there should exists a non trivial 
fixed point in the coupled system Maxwell-super. 


Equations (\ref{At}) and (\ref{J1}) or (\ref{J2}) do not support
such a non trivial fixed point. This is easy to see for (\ref{J1}).
For (\ref{J2}) , write the system (\ref{At},\ref{J2}) as
\begin{eqnarray}
\mathbf{A}_{t} = \mathbf{C} , \\
\mathbf{A}_{t} = c^2 \Delta \mathbf{A} 
-{1 \over \epsilon_0} \mathbf{J_0} {B \over B+B_0}  \tanh({\parallel \mathbf{A}_{t} \parallel \over E_c}) ,
\end{eqnarray}
where $\mathbf{C}$ is an auxiliary vector and where we used
the relation (\ref{vpot}).
Assume all quantities are scalars.
Looking for a fixed point, we write $A_t=C_t=0$ and find that
the only solution is $A=0$.
The conclusion of this simple analysis is that the constitutive
equations (\ref{J1}) or (\ref{J2}) are not adequate to describe
the effect reported in the experiments of flux jump.


We therefore suggest to use the full Ginzburg-Landau functional
to describe the superconductor. This is the object of the next section.


\section{The Maxwell-Ginzburg-Landau model}


The Ginzburg-Landau free energy can be written as (see \cite{tinkham}
for example)
\be\label{freegl}
F(\mathbf{A}, \psi) = -\alpha |\psi|^2 + {\beta \over 2} |\psi|^4 +
{1 \over 4 m } |(-i \hbar \nabla - {2 e \over c} \mathbf{A} )\psi |^2 . \ee
The $\alpha$ term is usually taken as $\alpha_0(T-T_c)$. We assume
$T < T_c$ so we put $-\alpha$.



We start with a 1D version of the problem to see if it makes sense.
For that, in the spirit of the much simpler study \cite{ckm06},
we write a Lagrangian for the Maxwell equations for the single
component $a$ of the vector potential and the 
Ginzburg-Landau potential 
$${\cal L}(a,a_t,a_x,\psi,\psi_t,\psi_x) = {a_t^2 \over 2 c^2}-{a_x^2 \over 2}
+I(x) [ {i\over 2}(\psi_t \psi^* -\psi \psi^*_t) +\alpha |\psi|^2 
- {\beta \over 2} |\psi|^4$$
\be\label{lag} - \xi^2 |\psi_x|^2 -i \mu a (\psi^*\psi_x -\psi \psi_x^*)
- \sigma a^2 |\psi|^2
- \gamma |\psi|^2 a_t ] ,\ee
where $I(x)$ is the indicator function of the superconductor, 
\be\label{constants}
\mu = {e \hbar \over 2 m},~~\sigma = {e^2 \over m},~~\xi^2 = {\hbar^2 \over 4 m } , \ee
and the $\gamma$ term is the coupling of the electric field $a_t$
to the charge in the superconductor. 
The kinetic term for $\psi$ is of the Schrodinger type because we want
a relaxation like dynamics.


The Euler-Lagrange equations are 
\begin{eqnarray}
\label{ela} \partial_t({\partial {\cal L} \over \partial a_t }) + 
\partial_x({\partial {\cal L} \over \partial a_x }) 
-{\partial {\cal L} \over \partial a } =0, \\
\label{elu} 
\partial_t({\partial {\cal L} \over \partial \psi^*_t }) + 
\partial_x({\partial {\cal L} \over \partial \psi^*_x }) 
-{\partial {\cal L} \over \partial \psi^* } =0 . 
\end{eqnarray}



I will detail calculations for completeness, we can comment them
later. 
$${\partial {\cal L} \over \partial a_t }= {a_t \over c^2} -I \gamma  |\psi|^2$$
$${\partial {\cal L} \over \partial a_x }= -a_x $$
$${\partial {\cal L} \over \partial a }= I[-2a \sigma |\psi|^2 -i \mu  (\psi^*\psi_x -\psi \psi_x^*)] $$
so that we get for $a$
\be\label{ela1} a_{tt} - c^2 a_{xx}= Ic^2[-2 a \sigma |\psi|^2 
-i \mu  (\psi^*\psi_x -\psi \psi_x^*)
+\gamma |\psi|^2_t], \ee


For $\psi$ we have
$${\partial {\cal L} \over \partial \psi^*_t }= -i I {\psi \over 2} , $$
$${\partial {\cal L} \over \partial \psi^*_x}=  I[-\xi^2 \psi_x +i \mu a \psi] , $$
$${\partial {\cal L} \over \partial \psi^* }= 
I[ i {\psi_t \over 2}  + \alpha \psi - \beta |\psi|^2 \psi -i \mu  a \psi_x - \sigma a^2 \psi - \gamma \psi a_t] , $$
so the final equation is
\be \label{elpsi1} I[i \psi_t +\xi^2 \psi_{xx} 
+\alpha \psi - \beta |\psi|^2 \psi -2i \mu  a \psi_x - \sigma a^2 \psi -i\mu a_x \psi- \gamma \psi a_t] -\delta(x)[ -\xi^2 \psi_{x} +i \mu a \psi]=0 . \ee


To make sense of these generalized PDE's let us integrate them
over one of the edges of the domain, say $x=-L$.
The first equation yields 
$$\int_{-L-\epsilon}^{-L+\epsilon} a_{tt} 
-[{a_x  c^2}]_{-L-\epsilon}^{-L+\epsilon}
= c^2\int_{-L}^{-L+\epsilon} (-2 a\sigma |\psi|^2 +\gamma |\psi|^2_t -i \mu  (\psi^*\psi_x -\psi \psi_x^*)) .$$
Taking the limit $\epsilon \to 0$ and assuming bounded variations
of the integrands we get that $a_x$ is continuous at the interface
$x=-L$.
The second equation gives
$$\int_{-L}^{-L+\epsilon} [i \psi_t +\xi ^2\psi_{xx} + \dots]
-(-\xi^2 \psi_{x} + i \mu a \psi)_{-L}=0$$
so that taking the limit $\epsilon \to 0$, we recover the standard
boundary condition \cite{vanduzer}
\be\label{ju}
-\xi^2 \psi_{x} + i \mu a \psi   = 0 , \ee
at $x=-L$.


The final system 
including the coupling conditions at the interfaces $x=\pm L$ is then 
\begin{eqnarray}
\label{ela2} a_{tt} - c^2 a_{xx}= Ic^2[-2 a \sigma |\psi|^2
-i \mu  (\psi^*\psi_x -\psi \psi_x^*)
+\gamma |\psi|^2_t], \\
\label{elpsi2} i \psi_t +\xi^2 \psi_{xx} -2i \mu a \psi_x
+(\alpha - \sigma a^2-i\mu a_x - \gamma a_t) \psi - \beta |\psi|^2 \psi  =0 , \\
\label{elapsi2} -\xi^2 \psi_{x} + i \mu a \psi   = 0,~~~x=\pm L .
\end{eqnarray}


\subsection{Fixed points}

The fixed points of the system of equations
(\ref{ela2}-\ref{elapsi2}) are given by
\begin{eqnarray}
\label{fla2} - c^2 a_{xx}= Ic^2[-2 a \sigma |\psi|^2
-i \mu  (\psi^*\psi_x -\psi \psi_x^*)], \\
\label{flpsi2} \xi^2 \psi_{xx} -2i \mu a \psi_x
+(\alpha - \sigma a^2-i \mu a_x) \psi - \beta |\psi|^2 \psi  =0 , \\
\label{flapsi2} -\xi^2 \psi_{x} + i \mu a \psi   = 0,~~~x=\pm L .
\end{eqnarray}


To simplify the problem, we first look for solutions such that $\psi_x=0$.
The equations reduce to
\begin{eqnarray}
\label{fca}  a_{xx}= I[2 a \sigma |\psi|^2 ], \\
\label{fcpsi} (\alpha - \sigma a^2-i \mu a_x) \psi - \beta |\psi|^2 \psi  =0 , 
\end{eqnarray}
Two solutions exist for $\psi$
\be\label{solpsi}
\psi=0,~~|\psi|^2 = {  \alpha - \sigma a^2 + i \mu a_x \over \beta}  . \ee
since $| \psi | ^2 \in \mathbf{R}$ and $a_x$ is continous at $\pm L$, $a_x$ must be zero on $[-L,L]$ \\
As a result $a$ is a constant on $[-L,L]$\\
From (\ref{solpsi}), since $|\psi|^2$ is real we have $a_x=0$ and 
\be \label{afix} a = \pm \sqrt{\alpha \over \sigma} \ee


{\bf A revoir} \\
Plugging the second expression into (\ref{fca}) we get
$$a_{xx} +  I a {2 \alpha \sigma \over \beta} ({\sigma \over \alpha} a^2 -1)=0 .$$
Recalling the London characteristic length $\lambda$ and introducing 
the amplitude $a_0$ as 
\be\label{london}
{1 \over \lambda^2} = {2 \alpha \sigma \over \beta} ,~~
{1 \over a_0^2} = {\sigma \over \alpha} ,  \ee
we get a simplified equation
\be\label{duffing}
a_{xx} +  I a ({a^2 \over a_0^2} -1)=0 ,\ee
where the unit of length is now $\lambda$. 


{\bf This should be solved using Newton for different values of a0
to see what type of solutions we get}


\section{Numerical solution of the 1D Maxwell-Ginzburg-Landau equation}

The dynamical equations (\ref{ela2}-\ref{elapsi2}) can be solved using 
the Finite Difference Time Domain (FDTD) algorithm in a similar 
way as in \cite{bcg19}. This algorithm was shown to be conditionally
stable.
The algorithm is the following {\bf THIS SHOULD BE ADAPTED}



\begin{eqnarray} 
&& \mu_0\frac{H^{n+\frac{1}{2}}_{j+\frac{1}{2}}-H^{n-\frac{1}{2}}_{j+\frac{1}{2}}}{\Delta t} = - \frac{E^{n}_{j+1}-E^{n}_{j}}{\Delta x} ,  \label{fdtd1} \\
        && \frac{D^{n+1}_{j}-D^{n}_{j}}{\Delta t} = -\frac{H^{n+\frac{1}{2}}_{j+\frac{1}{2}} - H^{n+\frac{1}{2}}_{j-\frac{1}{2}}}{\Delta x},  \label{fdtd2} \\
&& \frac{\tau^{2}}{\epsilon_{0}}\; \frac{P^{n+1}_{j}-2P^{n}_{j}+P^{n-1}_{j}}{\Delta t^{2}} + \tilde{\alpha} P^{n}_{j} + \tilde{\beta} (P^{n}_{j})^3 =E^{n}_{j},  \label{fdtd3} \\
        && E^{n+1}_{j}=\frac{1}{\epsilon_{0}}(D^{n+1}_{j}-P^{n+1}_{j}) ,  \label{fdtd4} \\
&& \text{where} \;\;\; \tilde{\alpha} = {\tau^2 \over \epsilon_0} \alpha, ~~~
\tilde{\beta} = {\tau^2 \over \epsilon_0} \beta.
\end{eqnarray}
An array $ \displaystyle H^{n+\frac{1}{2}}_{j+\frac{1}{2}}$ approximates an exact magnetic field evaluated at $x=(j+\frac{1}{2}) \Delta x$ and $t=n \Delta t$. The other arrays for $P, D, $ and $E$ are interpreted similarly.



The Matlab code is given in the Appendix.



\section{Numerical solution of the 2D Maxwell-London equation}

\subsection{Equation with A: cylindrical coordinates}

We assume $\mathbf{A} = a(z,r) \mathbf{e_{\theta}}$ and
denote in the following $a = a(z,r)$.
From Maxwell's equation and the London equation (\ref{At}) and (\ref{J1}) , we get
$$\mathbf{A}_{tt} - c^2\Delta \mathbf{A} = -I_{\Omega_d}\frac{c^2}{\lambda^2}\mathbf{A}$$
where $I_{\Omega_d}$ est the indicator function of the superconductor.
Using $\mathbf{A} = (0,0,a(z,r))$ in the $(z, r, \theta)$ cylindrical 
coordinates we deduce
$$\Delta\mathbf{A} = (\frac{\partial^2 a}{\partial z^2}+\frac{\partial^2 a}{\partial r^2}+\frac{1}{r}\frac{\partial  a}{\partial r}- \frac{ a}{r^2})\mathbf{e_\theta}.$$
The wave equation for $\mathbf{A}$ yields 
\be\label{att0}a_{tt}-c^2(\frac{\partial^2 a}{\partial z^2}+\frac{\partial^2 a}{\partial r^2}+\frac{1}{r}\frac{\partial a}{\partial r}- \frac{a}{r^2}) =- I_{\Omega_d}\frac{c^2}{\lambda^2} a . \ee


The magnetic field $\mathbf{B}$ is given by 
$$\mathbf{B} = \nabla \cdot \mathbf{A}$$
which gives
\be\label{boa} 
\mathbf{B} = (-\partial_z a, ~ \frac{1}{r}\partial_r(r a_r),~ 0) .\ee


\subsection{The dynamic problem : the London approach}

We are dealing with the following problem
\be\label{att2} a_{tt}-c^2(\frac{\partial^2 a}{\partial z^2}+\frac{\partial^2 a}{\partial r^2}+\frac{1}{r}\frac{\partial  a}{\partial r}- \frac{ a}{r^2}) =- I_{\Omega_d}\frac{c^2}{\lambda^2} a\ee
with the initial conditions for $a$ and $a_t$
\begin{eqnarray} \label{a0}
 a^0(z,r,t=0)={\exp}(-\frac{1}{2}(z-z_0)^2) , \\
\label{da0}
\frac{\partial  a^0}{\partial t}(z,r,t=0) = -c(z-z_0)\text{exp}(-\frac{1}{2}(z-z_0)^2) ,
\end{eqnarray} 
where $ a^0 = a(t=0)$. 
The boundary conditions on the different segments $\Gamma_i, ~i=1,\dots,8$
of $\partial \Omega$ (see Fig. \ref{schema}) are
\be\label{abound}
\frac{\partial  a}{\partial n} = 0 ~~{ \rm on}~~ \Gamma_{7,8,4,9,5,6} \ee
\begin{figure}[H]
\centering
\resizebox{10 cm}{5 cm}{\includegraphics[angle=0]{figs/schema.eps}}\hfill
\caption{The computational domain for the PDE problem 
(\ref{att2} ,\ref{a0},\ref{da0},\ref{abound}).}
\label{schema}
\end{figure}




In order to solve it numerically, we use an order 2 centered finite difference 
scheme in time.
$$ a_{tt} (t=n dt) \approx \frac{ a^{n+1}-2 a^n+ a^{n-1}}{dt^2},$$
where $ a^n= a(t=n dt)$.
We deduce the following variational formulation for $ n\geq 1$\\
{Find }  $a(.,t)\equiv  a^{n+1} \in V$, such that $\forall v \in V$
\be\label{watt}
\int_{\Omega} \nabla  a \nabla v + \frac{1}{c^2dt^2}\int_{\Omega} av + \int_{\Omega}(\frac{ av}{r^2} - \frac{1}{r} \frac{\partial  a}{\partial r}v)+\int_{\Omega}\frac{1}{\lambda^2}I_{\Omega_d} av +  \frac{1}{c^2dt^2}\int_{\Omega}(-2 a^n+ a^{n-1})v=0\ee
with  initial conditions 
 $$a^0={\exp}(-\frac{1}{2}(z-z_0)^2), ~~
 a^1 = {\exp}(-\frac{1}{2}(z-z_0+c dt)^2) $$


The results are displayed in figure \ref{At}
\begin{figure}[H]
\centering
\resizebox{15 cm}{5 cm}{\includegraphics[angle=0]{figs/A_instat10.pdf}
\includegraphics[angle=0]{figs/A_instat20.pdf}}
\resizebox{15 cm}{5 cm}{\includegraphics[angle=0]{figs/A_instat30.pdf}
\includegraphics[angle=0]{figs/A_instat40.pdf}}
%\resizebox{15 cm}{5 cm}{\includegraphics[width = 3in]{figs/ffglut_0000} \includegraphics[width = 3in]{figs/ffglut_0001}}
%\resizebox{15 cm}{5 cm}{\includegraphics[width = 3in]{figs/ffglut_0002} \includegraphics[width = 3in]{figs/ffglut_0003}}
\caption{Amplitude of the vector potential at times $t = 10 dt$ , $t = 20 dt$ (top left and right) and $t=30 dt$ and $t=40dt$ (bottom).}
\label{At}
\end{figure}


The code in Freefem++ is given in the appendix.


%{\bf Acknowledgements}
%The calculations were done at the CRIANN computing center.

\begin{thebibliography}{99}



\bibitem{cgl13} 
J.-G. Caputo, L. Gozzelino, F.Laviano, G.Ghigo, R.Gerbaldo, J.Noudem, Y.Thimont, P.Bernstein, 
"Screening magnetic fields by a superconducting disk: a simple model", 
J. Appl. Phys. 114, 233913 (2013) 

\bibitem{feynman}  R. Feynman, course of physics, Electromagnetism.
\bibitem{tinkham} M. Tinkham, "Introduction to Superconductivity", second edition, Dover, (2004).
\bibitem{vanduzer} T. Van Duzer and C. W. Turner, "Principles of superconductive devices and circuits", Edward Arnold, (1998).
\bibitem{bean} C.P. Bean, Rev. Mod. Phys. 36, 31 (1964).

\bibitem{bcg19} M. Brio, J.G. Caputo, K. Gwirtz, J. Liu and A. Maimistov,
 " Scattering of a short electromagnetic pulse from a Lorentz-Duffing film: theoretical and numerical analysis",
Wave motion, 89, 43-56, (2019).
\\https://arxiv.org/abs/1804.07061

\bibitem{weinstein} R.Weinstein et al. IEEE Transactions on Applied 
Superconductivity 25 , 6601106 (2014)

\bibitem{ckm06} J.-G. Caputo, E. V. Kazantseva, A.I. Maimistov,
''Electromagnetically induced switching of ferroelectric thin films'',
Phys. Rev. B 75, 014113 (2007)
\\http://arxiv.org/abs/nlin.PS/0607001


\end{thebibliography}


\end{document}
\section{Appendix: Freefem code for 2D Maxwell-London} 
\begin{verbatim}
// calculus of the vector potential A for the London model:
// 1/c2*Att - Lap(A) = -Indic*1/lambda2*A
// boundary condition : a gaussian propagating pulse
// a cylindrical shaped superconducter - boundary : [0, Rs] * [-Zs, Zs]
// units : millimeters
// Omega = external boundary = ]0;Rmax[*]-Zmax;Zmax[ 

//===============================================================
// Parametres
//===============================================================
real lambda = 1.6e-4; // en mm
real lambda2 = lambda*lambda;

real Rmax = 20;
real Zmax = 15;
real Rs = 10;
real Zs = 1;
real c = 3e8;  // speed in mm/s
real c2 = c*c;
int nbsegul = 5; // nb of segment per unit length
real h = 1./nbsegul;
real cfl = 0.8;
real dt = cfl*h/c;
real dt2 = dt*dt;
real T = 120*dt;
cout<<"dt = "<<dt<<endl;
cout<<"c2*dt2 = "<<c2*dt2<<endl;

//===============================================================
// mesh 
//===============================================================


	//external boundary
border Gamma1(t=0,Rs){x=t; y=-Zs;label=1;};
border Gamma5(t=0,Rmax){x=t; y=-Zmax;label=5;};
border Gamma6(t=-Zmax,Zmax){ x=Rmax; y=t; label=6;}; 
border Gamma7(t=0,Rmax){ x=Rmax-t; y=Zmax; label=7;} ;
border Gamma8(t=0,Zmax-Zs){ x=0; y=Zmax-t; label=8;} ;
border Gamma4(t=-Zs,Zs){ x=0; y=-t; label=4;} ;
border Gamma9(t=0,-Zmax+Zs){ x=0; y=-Zs+t; label=9;} ;
	//superconductor
border Gamma2(t=-Zs,Zs){x=Rs;y = t; label=2;};
border Gamma3(t=0,Rs){x=Rs-t;y=Zs; label=3;};
	

int nbseg = Rmax*nbsegul;
mesh Th=buildmesh( Gamma1(nbseg*Rs/Rmax)+Gamma5(nbseg)+Gamma6(nbseg*2*Zmax/Rmax)+Gamma7(nbseg)+Gamma8(nbseg*(Zmax-Zs)/Rmax)+Gamma4(nbseg*2*Zs/Rmax)+Gamma2(nbseg*2*Zs/Rmax)+Gamma3(nbseg*Rs/Rmax)+Gamma9(nbseg*(Zmax-Zs)/Rmax));
plot(Th,wait=1, cmm="maillage du domaine");



//===============================================================
// Finite element space
//===============================================================
fespace Vh(Th,P1); 

// subdomain identification
fespace Ph(Th,P0);
	
//labels of subdomains
	Ph labelR=region;
	int label1=labelR(Rs/2,Zs/2);  // superconductor
	int label2=labelR(Rmax,Zmax);  // upper right corner
	cout<<"labels domaines ="<<label1<<"  et "<<label2<<endl;

//domain of superconductor
Ph IdSupra = (region == label1) ? 1 : 0 ; 
plot(IdSupra, cmm="fonction indicatrice du supraconducteur",wait=1,fill=1); 

//=========================
// building matrices
//=========================

//--------- matrix for the system
Vh u,v;
real tgv = 1e30;
varf Vsys(u,v) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v)+IdSupra*u*v/lambda2)
						-int2d(Th)(dx(u)*v/x)
						+int2d(Th)(u*v/x^2)
						+int2d(Th)(1./(c2*dt2)*u*v)
						+on(5,7,u=1);
matrix Asys=Vsys(Vh,Vh,tgv=tgv);
set(Asys,solver = UMFPACK);

//--------- matrix for the right hand side
varf Vmass(u,v)=int2d(Th)(1./(c2*dt2)*u*v)+on(5,7,u=1);
matrix Amass = Vmass(Vh,Vh,tgv = tgv);

//--------- vector for the right hand side
Vh BC;
varf Vbc(A,v) = on(5,7,A=1);
BC[] = Vbc(0,Vh,tgv=1);
BC[] = -BC[];
BC[] += 1; // 1 everywhere except 0 on boundaries 5,7


//===============================================================
// Variationnal formulation
// finite difference scheme of order 2 in time
// variationnal formulation in space
//===============================================================
   Vh A,A0,A1;
   A0 = exp(-0.5*(y-0.5*Zmax)*(y-0.5*Zmax)); 
   //A1 = exp(-0.5*(y-0.5*Zmax+c*dt)*(y-0.5*Zmax+c*dt));	
   A1 = A0-c*dt*(y-0.5*Zmax)*exp(-0.5*(y-0.5*Zmax)*(y-0.5*Zmax));

   real t=dt;
   int it = 1;
   real [int] rhs = A0[];
   real cputime = clock();
   //plot(A0,value = true,cmm="temps="+0+"  A.max="+A0[].max,ps = "A_instat"+0+".eps");
  // plot(A1,value = true,cmm="temps="+t+"  A.max="+A1[].max,ps = "A_instat"+1+".eps");
while(t < T)
{
  t += dt; 
  it += 1; 
  A[] =  2*A1[]-A0[];
  rhs = Amass*A[];
  rhs .*= BC[];
  A[] = Asys^-1*rhs;
  if(it%30==0)
  {
  cout <<"=========== temps="<<t<<endl;
  plot(A,value = true,cmm="temps="+t+"  A.max="+A[].max,ps = "A_instat"+it+".eps",grey=1,wait=true);
  }
  A0[]=A1[];
  A1[]=A[];
}
cout<<"temps (varf+) = "<<clock()-cputime<<endl;
\end{verbatim}


\section{Appendix: FDTD code } 

\begin{verbatim}

%------------------------------------------------ 
%   1D FDTD Lorentz media
%   (Ey, Hz) & Py. Pulse propagates in x direction
%
%   Original Maxwell:
%       d(Dy)/dt = -d(Hz)/dx    (1)
%    mu*d(Hz)/dt = -d(Ey)/dx    (2)
%             Dy =  eps*Ey+P    (3)
%
%   P: Lorentz media
%   (tau^2/eps0) P_tt + (Gamma tau^2/eps0) P_t + alpha P 
%      + beta P^3 = E    (4)
%
%
%   Algorithm
%    a)  advance H, H_t = ..
%    b)  advance P  P_tt = ..
%    c)  advance D  D_t = ..
%    d)  compute E  E = (D-P)/eps
%
%    (4) rewritten as
%      delta P_tt + xi P_t + alpha P + beta P^3 = E (5)
%    discretized as
%    delta (P_n+1 + P_n-1 -2P_n)/dt^2 + xi (P_n+1 - P_n-1)/(2dt) 
%     + alpha P_n + beta P_n^3 = E_n
%     in layer   jl < n < jr
% 
%------------------------------------------------
function fdtd1d_Lorentz
clear all;

display('Program start ...')

%Tot_steps = 40000;                  % Total number of time steps
Tot_steps = 20000;                  % Total number of time steps
ns =400 ;                            % Nb of snapshots
nsamp = Tot_steps/ns;
%%pwid = 14e-15;%/(2*sqrt(log(2)));   % Pulse width
pwid = 14e-15;%/(2*sqrt(log(2)));   % Pulse width
pcen = pwid*3;                      % Pulse center
Start_time = pcen*2;                


eps0= 8.854e-12;        % permittivity in vacuum
mu0 = pi*4e-7;          % permeability in vacuum    
c   = 299792458;        % light speed in vacuum c = 299792458 meter/second
Z0  = sqrt(mu0/eps0);   

% parameters dielectric slab 
alpha   = 2.2e10;
tau     = 3.e-16;  % tau from Prof Caputo's paper (tau =1.5e-17  blow up)
Gamma   = 0.;      % no damping
beta    = 0;
%beta    = 1.e29 ;    % nonlinearity  beta P^3
% width dielectric slab 
wp      = 5.e-6;


%%lambda      = 2*pi*tau*c/sqrt(eps0*alpha);
lambda      = 1.4e-6;     % wavelength
omega       = 2*pi*c/lambda;% angle frequency  
freq        = c/lambda;     % 2GHz pulse
period      = lambda/c;     % period
wavenumber  = 2*pi/lambda;  % wavenumber

% Computational Domain Lx < x < Ux, grid size Nx
Lx=0e-8;    
Ux=100e-6;
%Nx=10000;   % number of grid points 
Nx=2500;   % number of grid points 

dx=(Ux-Lx)/Nx;
dt=0.25*dx/c; 


% output files fab 
fab = fopen('eab','w');
ja = round(Nx/5);
jb = round(4*Nx/5);

% width dielectric slab 
jl = Nx/2-round(wp/(2*dx));
jr = Nx/2+round(wp/(2*dx));


xx = linspace(Lx,Ux,Nx+1); xx = xx'; % xx(i) is the x-coordinates of the i-th grid point

PML_size = 100;
sg_max = 32./(PML_size*dx*Z0);

Dy      = zeros(Nx+1,1);
Ey      = zeros(Nx+1,1);
Hz      = zeros(Nx+1,1);
Py      = zeros(Nx+1,1);
Pyold   = zeros(Nx+1,1);
Pynew   = zeros(Nx+1,1);
Ptime   = zeros(Tot_steps,1);
time    = linspace(0,dt*Tot_steps,Tot_steps);

E0 = 1; 

%eps_s = 5;
eps_infty = 1;
epsEy   = eps_infty*eps0*ones(Nx+1,1);
Z2  = mu0/epsEy;


% polarisation parameters 
delta = tau^2/eps0 ;
xi    = Gamma*tau^2/eps0;

% Lorentz
den = 2*delta + xi*dt;
pc1 =   (4*delta-2*alpha*dt^2)/den;  
pc2 =   (xi*dt -2*delta)/den;
pc3 =   2*dt^2 / den; 
pc4 =   -2*beta*dt^2 / den;
    
px1 = Lx;
px2 = Lx + PML_size*dx;
px3 = Ux - PML_size*dx;
px4 = Ux;
ppx = [px1,px2,px3,px4];

sgx     = sg_max*find_sigma(xx,       ppx);
sgcx    = sg_max*find_sigma(xx+0.5*dx,ppx);
cHz1    = (2*mu0   - dt*sgcx*Z2)./(2*mu0    + dt*sgcx*Z2);
cHz2    =                   2*dt./(2*mu0    + dt*sgcx*Z2);
cEy1    = (2*epsEy - dt*sgx   )./(2*epsEy  + dt*sgx);
cEy2    =                  2*dt./(2*epsEy  + dt*sgx).*epsEy;
    
disp('program is running ... ')

for step=0:Tot_steps
            
    tt = (step)*dt;
    Ptime(step+1)=Py(round(Nx/2));
    if(mod(step,nsamp)==0)
        figure(1);
        plot(xx*1e6,Ey,'b-');
%        plot(xx(jl:jr)*1e6,Ey(jl:jr),'b-');
        xlabel('x(\mum)')
        ylabel('Ey (V/m)')
        set(gca,'ylim',[-3 3]*1.1)
        getframe;
        figure(2);
        plot(xx(jl:jr)*1e6,Py(jl:jr),'b-');
        xlabel('x(\mum)')
        ylabel('Py (V/m)')
        set(gca,'ylim',[-1e-10 1e-10]*1.1)
        getframe;
    end
%for graphics
%    if(tt > 5.e-13)
%        stop
%    end
    
% incident source
    if(tt < Start_time)
        Ey(PML_size+2,:) = Ey(PML_size+2,:)+E0*sin(omega*tt)*exp(-((tt-pcen)/(pwid))^2);
        Dy(PML_size+2,:) = eps0*Ey(PML_size+2,:);%
%%        Ey(Nx/2,:) = Ey(Nx/2,:)+E0*sin(omega*tt)*exp(-((tt-pcen)/(pwid))^2);
%%        Dy(Nx/2,:) = eps0*Ey(Nx/2,:);%
    end
    
% Update Hz: d(Hz)/dt = 1/mu * (d(Ex)/dy - d(Ey)/dx)
    Hz(1:Nx,1) =  cHz1(1:Nx,1).*Hz(1:Nx,1) - ...
                  cHz2(1:Nx,1).*(Ey(2:Nx+1,1) - Ey(1:Nx,1))/dx;

% Lorentz, Update Py
    Pynew(jl:jr) = pc1*Py(jl:jr) + pc2*Pyold(jl:jr) + pc3*Ey(jl:jr)+ pc4*Py(jl:jr).^3;
    Pyold(jl:jr) = Py(jl:jr);
    Py(jl:jr) = Pynew(jl:jr);

%% for debugging
    fprintf('%5.2e  %5.2e  %5.2e \n',tt,Ey(jl),Pynew(jl));
% write Ea Eb
    fprintf(fab,'%5.2e  %15.10g %15.10g \n',tt,Ey(ja),Ey(jb));
    
    % Update Dy: d(Dy)/dt = - d(Hz)/dx
    Dy(2:Nx,1) = cEy1(2:Nx,1).* Dy(2:Nx,1) - ...
                 cEy2(2:Nx,1).*((Hz(2:Nx,1) - Hz(1:Nx-1,1))/dx);
    
    % Update Ey   
    Ey = (Dy - Py)/eps_infty/eps0;
      
end

% end of main function


%% Subroutines

function sgs = find_sigma(crds, ran)
    sgs=zeros(Nx+1,1);
    sgs(:,1) = max(0,max(   (ran(2)-crds(:,1))/(ran(2)-ran(1)),...
                            (crds(:,1)-ran(3))/(ran(4)-ran(3)) ) ).^3;
end


end
\end{verbatim}

\end{document}


\section{Appendix First approach : 1D model}

Using (\ref{At}) and (\ref{J1}), we assume $\mathbf{A}$ to have
a single azimuthal component $A$
and write
\be\label{At1d}
{A}_{tt} - c^2 {A}_{xx}= I_\delta {A \over \tau^2} ,\ee
where $\tau = \lambda / c $ is the relaxation time of the superconductor and $ I_\delta$
is the indicator function of the spatial interval $[0 ;L]$.
To solve (\ref{At1d}), we separate time and space by writing
$$A(x,t) = e^{i \omega t} a(x)  .  $$
We obtain
\be\label{a1d}
a_{xx} + [k^2 - {1 \over \lambda^2} I_\delta ]a=0 . \ee


To model an incoming pulse of magnetic field, we set up a scattering 
experiment.
\begin{eqnarray}
\label{sca1} a = e^{-ikx} + R e^{ikx}, ~~~ x \le  0 ,\\
\label{sca2} a = A \cos k_0 x + B \sin k_0 x, ~~~ 0 \le x \le L,  \\
\label{sca3} a = T e^{-ikx} , ~~~  L \le x  ,
\end{eqnarray}
where $k_0= \sqrt{k^2 - {1 \over \lambda^2}}$.
From equation (\ref{a1d}) we have at the interfaces $x=0,L$
continuity of $a$ and $a_x$.
As in \cite{bcg19} this yields a linear system whose solution is
\begin{flalign}\label{ref_trans}
& R= {(k^2-k_0^2) \sin k_0 L \over D}, ~~~~
T= {2i k k_0 \exp i k L \over D}, && \\
& A= {2 k (-i k_0 \cos k_0 L + k \sin k_0 L) \over D},~~~~
B= {2 k ( k \cos k_0 L + i k_0 \sin k_0 L) \over D}, && \nonumber \\
& D = -2i k k_0 \cos k_0 L + (k^2+k_0^2) \sin k_0 L  . &
\end{flalign}


{\bf Continue this ...}


