# Gating ferromagnetic resonance by superconductors via modulated reflection
of magnetization-radiated electric fields

Xi-Han Zhou  School of Physics, Huazhong University of Science and Technology,
Wuhan 430074, China  Tao Yu<EMAIL_ADDRESS>School of Physics, Huazhong
University of Science and Technology, Wuhan 430074, China

###### Abstract

We predict that ferromagnetic resonance in insulating magnetic film with
inplane magnetization radiates electric fields polarized along the
magnetization with opposite amplitudes at two sides of the magnetic insulator,
which can be modulated strongly by adjacent superconductors. With a single
superconductor adjacent to the magnetic insulator this radiated electric field
is totally reflected with a $\pi$-phase shift, which thereby vanishes at the
superconductor side and causes no influence on the ferromagnetic resonance.
When the magnetic insulator is sandwiched by two superconductors, this
reflection becomes back and forth, so the electric field exists at both
superconductors that drives the Meissner supercurrent, which in turn shifts
efficiently the ferromagnetic resonance. We predict an ultrastrong coupling
between magnons in the yttrium iron garnet and Cooper pairs in NbN with the
frequency shift achieving tens of percent of the bare ferromagnetic resonance.

##  I Introduction

“Magnonics” exploits magnetic excitations, i.e., spin waves or their quanta,
magnons, as potential information carriers for spin transport in insulators
with low-energy consumption [1, 2, 3, 4, 5, 6, 7]. Interaction between magnons
and Cooper pairs in heterostructures composed of magnets and superconductors
may modulate the transport of spin information [8, 11, 17, 13, 15, 16, 9, 10,
14, 12], strongly enhance the magnon-photon interaction [18, 22, 19, 20, 21,
23], and lead to the emergence of triplet Cooper pairing [26, 24, 25, 28, 27],
which is potential to bring unprecedented functionalities in spintronics [26,
24, 25], quantum information [29, 30, 31, 32, 33], and topological quantum
computation [34]. In this heterostructure, the hybridized quantum states and
distribution of macroscopic electromagnetic fields govern its properties. For
example, the “ultrastrong coupling” [35] with the coupling strength close to
the ferromagnetic resonance (FMR) frequency unveils the importance of the
dipolar interaction in the superconductor(S)$|$metallic
ferromagnet(F)$|$superconductor(S) heterostructure [19, 20], where the photon
mode with a large mode density is localized in the nano-scale between two
superconductors [36].

The importance of the dipolar interaction also manifests in the superconductor
gating effect on magnons [37, 38, 39, 40, 12, 11, 17], in which the frequency
of magnons with finite wave number [41, 42, 43, 44] can be shifted up to tens
of GHz, as recently predicted [11, 12] and observed [17] in the
superconductor(S)$|$ferromagnet(F) insulator heterostructure. The stray
electric field of magnons drives the supercurrent in the adjacent
superconductor which in turn generates the Oersted magnetic field that affects
the low-frequency magnetization dynamics. This gating effect favors the spin
diode [7, 45] and magnon trap [46, 47, 48] in proper gating configurations.
The FMR frequency in this S$|$F bilayer is not affected, however.

Figure 1: Snapshots of magnetization-radiated electric fields in different
heterostructure configurations. The electric field changes linearly when
across the thickness of the ferromagnetic insulating film. (a) The electric-
field amplitude is opposite at two sides of the thin magnetic insulator. (b)
When fabricating a superconductor thin film on a ferromagnetic insulator, the
electric field is suppressed to vanish at the superconductor side but enhanced
at the other side of the magnet. (c) and (d) When the magnet is sandwiched by
two superconductors, the electric field exists but differs at both sides in
both symmetric and asymmetric configurations.

On the other hand, the FMR of the metallic ferromagnet sandwiched by two
superconductors was shifted up to 50 mT in the resonant field when the
thickness of two superconductor layers is larger than the London’s penetration
depth, as observed in several recent experiments [49, 50, 51]. Above the
superconducting transition temperature, the FMR frequency recovers to the
Kittel mode [52], which may be exploited to realize the magnetic logic gate
through a phase transition in the superconductor. This phenomenon may be
related to the frequency splitting induced by spin-triplet superconducting
state [49], Meissner screening [51], and giant demagnetization effects [53,
13]. It appears that this modulation could be absent for the FMR in the
ferromagnetic insulators, however, which has not been reported yet [49, 51,
53, 13]. The experiment [49] showed that inserting a thin insulator layer in
the heterostructures composed of a metallic ferromagnet sandwiched by two
superconductors completely suppresses the shift of FMR. This raises the issue
of whether the FMR can be gated or not in magnetic insulators by adjacent
superconductors in proper configurations.

In this work, we study this issue by going beyond the quasi-static
approximation for magnetostatic modes [54] and demonstrate that although the
stray magnetic field of Kittel magnon with uniform magnetization precession is
vanishingly small outside of the in-plane magnetized ferromagnetic insulating
film, the radiated electric field is significant with opposite amplitudes at
two sides of the magnetic film and polarization parallel to the magnetization
direction. This distribution of the radiated electric field is sensitive to
the adjacent superconductors due to the total reflection, as illustrated in
Fig. 1 for snapshots of the distribution of electric fields in different
heterostructure configurations. The electric field is opposite at two sides of
a single thin ferromagnetic insulator [Fig. 1(a)]; contra-intuitively, in the
S$|$F bilayer this electric field is suppressed to vanish at the
superconductor side [Fig. 1(b)]; nevertheless, when sandwiched by two
superconductors, the electric field is neither shifted to vanish nor screened
completely, as plotted in Figs. 1(c) and (d) for symmetric and asymmetric
configurations. These features are well understood by our mechanism of
modulated reflection of magnetization-induced electric fields by
superconductors, which predicts the absence of FMR shift in ferromagnetic
insulator$|$superconductor heterostructure and the ultrastrong modulation of
FMR, shifted up to tens of percent of the bare frequency when the
ferromagnetic insulator is sandwiched by two thin superconductors.

This paper is organized as follows. We address the model and general formalism
in Sec. II. In Sec. III, IV, and V, we analyze the distribution of the
electric fields from FMR of a single ferromagnetic insulator, S$|$F bilayer,
and S$|$F$|$S heterostructure, respectively, and address the ultrastrong
interaction between the FMR and supercurrent. We conclude and discuss in Sec.
VI.

##  II Model and general formalism

We consider a heterostructure composed of a ferromagnetic insulating film of
thickness $2d_{F}\sim O(50~{}{\rm nm})$ with inplane magnetization sandwiched
by two thin superconductor layers with thickness $d_{S}\lesssim\lambda$ and
$d_{S}^{\prime}\lesssim\lambda$, respectively, as illustrated in Fig. 2. Here
$\lambda\sim O(100~{}{\rm nm})$ is London’s penetration depth of conventional
superconductors. In the ferromagnetic insulators, the dynamics of
magnetization ${\bf M}=M_{x}\hat{\bf x}+M_{y}\hat{\bf y}+M_{0}\hat{\bf z}$,
where $M_{0}$ is the saturated magnetization, is phenomenologically governed
by the Landau-Lifshitz equation [55]

$\displaystyle\partial{\bf M}/\partial t=-\mu_{0}\gamma{\bf M}\times{\bf
H}_{\text{eff}},$ (1)

where $\mu_{0}$ is the vacuum permeability and $-\gamma$ is the electron
gyromagnetic ratio. The magnetization precesses around the effective magnetic
field ${\bf H}_{\text{eff}}={\bf H}_{\text{app}}+{\bf H}_{d}+{\bf H}_{s}$ that
contains the external static field ${\bf H}_{\text{app}}=H_{0}\hat{\bf z}$,
the dipolar magnetic field ${\bf H}_{d}$ generated by the magnetic charge
$-\nabla\cdot{\bf M}$, and the Oersted magnetic field ${\bf H}_{s}$ from the
superconductor that needs a self-consistent treatment with Maxwell’s equation
[11]. At low frequencies, it is sufficient to express the stray magnetic field
as [55, 54]

$\displaystyle{\bf H}_{d,\beta}({\bf
M})=\dfrac{1}{4\pi}\partial_{\beta}\partial_{\alpha}\int d{\bf
r^{\prime}}\dfrac{M_{\alpha}(\bf r^{\prime})}{|\bf r-r^{\prime}|}.$ (2)

The exchange interaction plays no role in the FMR since the gradient of ${\bf
M}$ vanishes for the uniform precession.

Figure 2: S(1)$|$F$|$S(2) heterostructure. The thickness of superconductors
above and beneath the thin ferromagnetic insulator of thickness $2d_{F}$ is
$d_{S}$ and $d_{S}^{\prime}$, respectively. The driven supercurrents ${\bf
J}_{s}$ and ${\bf J}_{s}^{\prime}$ by FMR flow oppositely along the
magnetization direction.

On the other hand, the oscillating magnetic induction ${\bf B}=\mu_{0}({\bf
M}+{\bf H})$ of frequency $\omega$ governs the radiation of electric fields
inside and outside the ferromagnetic insulator according to [56]

$\displaystyle\nabla\times{\bf E}=i\omega{\bf B},$
$\displaystyle\nabla\times{\bf H}={\bf J}_{s}-i\omega\varepsilon_{0}{\bf E},$
(3)

where $\varepsilon_{0}$ is the vacuum permittivity. When coupled with
superconductors, this electric field drives the supercurrent ${\bf J}_{s}$ via
London’s equation [57]

$\displaystyle\dfrac{\partial{\bf J}_{s}}{\partial
t}=\dfrac{1}{\mu_{0}\lambda^{2}}{\bf E},$ $\displaystyle\nabla\times{\bf
J}_{s}=-\dfrac{1}{\mu_{0}\lambda^{2}}{\bf B}.$ (4)

Here London’s penetration depth at different temperatures $T<T_{c}$ follows
the relation [57]

$\displaystyle\lambda(T)=\sqrt{\frac{m_{e}}{\mu_{0}n_{e}e^{2}}}\left(1-\left(\frac{T}{T_{c}}\right)^{4}\right)^{-1/2},$
(5)

where $m_{e}$ is the electron mass and $n_{e}$ is the electron density.

At low frequencies, we may apply the quasi-static approximation
$\nabla\times{\bf B}=\mu_{0}{\bf J}_{s}$ in superconductors. Taking the curl
of Eq. (3) and substituting Eq. (4) into it, the electric field inside the
superconductor obeys

$\displaystyle\nabla^{2}{\bf E}-{\bf E}/\lambda^{2}=0.$ (6)

On the other hand, taking the curl of $\nabla\times{\bf B}=\mu_{0}{\bf J}_{s}$
and combining with Eq. (4), the magnetic induction inside the superconductor
obeys $\nabla^{2}{\bf B}-{\bf B}/\lambda^{2}=0$.

The driven supercurrent then affects the magnetization dynamics. From Eq. (4),
the electric field drives supercurrent inside the superconductor, which then
generates the vector potential. With the uniform magnetization precession, the
system is translational invariant in the $y$-$z$ plane, so the supercurrent
only depends on $x$ and as it for the vector potential [56]

$\displaystyle{\bf A}(x)=\dfrac{\mu_{0}}{4\pi}\int d{\bf
r^{\prime}}\dfrac{{{\bf J}_{s}}(x^{\prime})}{|{\bf r}-{\bf r}^{\prime}|}.$ (7)

Accordingly, the Oersted magnetic field ${\bf
H}_{s}=(1/\mu_{0})\nabla\times{\bf A}$ only contains the $y$-component
${H}_{y}=-\partial_{x}A_{z}(x)/\mu_{0}$, which drives the magnetization.

The boundary condition describes the fields at the interfaces [56]. For the
magnetic induction and field, ${\bf B}_{\perp}$ and ${\bf H}_{\parallel}$ are
continuous at the boundaries. Since there is no surface current or charge
accumulation, the electric field $\bf E$ is continuous at interfaces.

##  III Single thin ferromagnetic insulator

We start with a single insulating ferromagnetic film to address the
significant radiated electric fields from the uniform magnetization
precession. For a single ferromagnetic insulator of thickness $2d_{F}$ biased
by a static magnetic field ${\bf H}_{\rm app}=H_{0}\hat{\bf z}$, the
magnetization $\bf M$ for the FMR is uniform inside the ferromagnetic layer by
the constant demagnetization factor $N_{xx}=-1$. Since the magnetic film is
sufficiently thin, we stick to the uniform precession throughout this work.
The opposite magnetic charges at the two surfaces of the film generate
opposite magnetic field outside, which results in vanished stray magnetic
field ${\bf H}_{d}=0$ outside the ferromagnetic layer, as also calculated from
Eq. (2); inside the ferromagnet, ${\bf H}_{d}=\\{-M_{x},0,0\\}$ and ${\bf
B}=\\{0,\mu_{0}M_{y},\mu_{0}(H_{0}+M_{0})\\}$, in which only the $y$-component
of $\bf B$ oscillates with frequency $\omega$ that can radiate the electric
field.

###  III.1 Full solution

Here we go beyond the quasi-static approximation and solve the radiated
electric field. According to Eq. (3), the oscillating electromagnetic field is
the source for radiating microwaves in space. Taking the curl of the first
equation in Eq. (3), the electric field obeys

$\displaystyle\nabla^{2}{\bf E}+\varepsilon_{0}\mu_{0}\omega^{2}{\bf
E}=-i\omega\mu_{0}\nabla\times{\bf M},$ (8)

which has the solution

$\displaystyle{\bf E}({\bf
r})=\dfrac{i\mu_{0}\omega}{4\pi}\int\dfrac{[\nabla^{\prime}\times{\bf M}({\bf
r^{\prime}})]e^{ik|{\bf r-r^{\prime}}|}}{|{\bf r-r^{\prime}}|}d{\bf
r^{\prime}},$ (9)

where $k=\omega\sqrt{\mu_{0}\varepsilon_{0}}$ is the wave number of
microwaves. Since only the $x$ and $y$ components of ${\bf M}$ oscillate with
frequency $\omega$ and $\bf M$ is uniform inside the ferromagnetic layer,
$(\nabla\times{\bf M})_{x,y}=0$ in all space, leading to $E_{x}=E_{y}=0$ and

$\displaystyle
E_{z}(x)=\dfrac{i\mu_{0}\omega}{4\pi}\int\dfrac{[\partial_{x^{\prime}}M_{y}({\bf
r^{\prime}})]e^{ik|{\bf r-r^{\prime}}|}}{|{\bf r-r^{\prime}}|}d{\bf
r^{\prime}}.$ (10)

Using Weyl identity [7]

$\displaystyle\frac{e^{ik|{\bf r}-{\bf r^{\prime}}|}}{|{\bf r}-{\bf
r}^{\prime}|}=\int
dk_{z}^{\prime}dk_{y}^{\prime}\frac{ie^{ik_{z}^{\prime}(z-z^{\prime})+ik_{y}^{\prime}(y-y^{\prime})}e^{i\sqrt{k^{2}-k_{z}^{\prime
2}-k_{y}^{\prime 2}}|x-x^{\prime}|}}{2\pi\sqrt{k^{2}-k_{z}^{\prime
2}-k_{y}^{\prime 2}}},$ (11)

we obtain the electric field

$\displaystyle E_{z}=\dfrac{-\mu_{0}\omega
M_{y}}{2k}\begin{cases}e^{ik(x+d_{F})}-e^{-ik(x-d_{F})},&-d_{F}<x<d_{F}\\\
e^{ik(x+d_{F})}-e^{ik(x-d_{F})},&x>d_{F}\\\
e^{-ik(x+d_{F})}-e^{-ik(x-d_{F})},&x<-d_{F}\end{cases}.$ (12)

From Eq. (3), we find the magnetic induction $B_{x}=0$,
$B_{z}=\mu_{0}(H_{0}+M_{0})$ is static, and
$B_{y}=-\partial_{x}E_{z}/(i\omega)$ follows

$\displaystyle
B_{y}=\dfrac{\mu_{0}M_{y}}{2}\begin{cases}e^{ik(x+d_{F})}+e^{-ik(x-d_{F})},&-d_{F}<x<d_{F}\\\
e^{ik(x+d_{F})}-e^{ik(x-d_{F})},&x>d_{F}\\\
-e^{-ik(x+d_{F})}+e^{-ik(x-d_{F})},&x<-d_{F}\end{cases}.$ (13)

We are interested in the field near the ferromagnet with a distance
$\sim\lambda$. In ferromagnetic insulators, $\omega\sim 100$ GHz [11], and
$\lambda\sim 100$ nm for conventional superconductors, so $k\lambda\sim
3\times 10^{-5}\ll 1$. When $kx\rightarrow 0$, we have

$\displaystyle E_{z}(x)=\begin{cases}-i\mu_{0}\omega M_{y}x,&-d_{F}<x<d_{F}\\\
-i\mu_{0}\omega M_{y}d_{F},&x>d_{F}\\\ i\mu_{0}\omega
M_{y}d_{F},&x<-d_{F}\end{cases},$ (14)

as plotted in Fig. 1(a) for a snapshot. The magnetic induction

$\displaystyle B_{y}(x)=\begin{cases}\mu_{0}M_{y},&-d_{F}<x<d_{F}\\\
0,&x>d_{F}\\\ 0,&x<-d_{F}\end{cases},$ (15)

recovers to the results from quasi-static approximation [54] with vanished
magnetic field $H_{y}$ outside of the ferromagnet.

###  III.2 Quasi-static approximation

The above analysis implies that when focusing on the near-field limit, we may
apply the quasi-static approximation that sets $\nabla\times{\bf H}=0$ in Eq.
(3). When focusing on the FMR case, $\bf E$ is translation invariant in the
$y$-$z$ plane. i.e., $\partial_{z}E_{x}=0$. Taking the $y$-component of Eq.
(1), the oscillation of $B_{y}$ only generate $E_{z}$ parallel to the
magnetization:

$\displaystyle-\partial_{x}E_{z}=i\omega\mu_{0}M_{y}.$ (16)

Integrating along $x$ across the ferromagnet yields

$\displaystyle E_{z}(x)=-i\omega\mu_{0}M_{y}(x+d_{F})+E_{z}(x=-d_{F}).$ (17)

Thereby, $E_{z}$ depends linearly on $x$ inside the ferromagnet. Outside the
ferromagnet,

$\displaystyle E_{z}(x)=-2i\omega\mu_{0}M_{y}d_{F}+E_{z}(x=-d_{F})$ (18)

is uniform, which is consistent with the vanished magnetic field
$H_{y|\text{outside}}=0$ in the quasi-static approximation. According to the
symmetry, $E_{z}(x=0)=0$, so the electric field is exactly the same as Eq.
(14).

##  IV S$|$F heterostructure

We consider the S$|$F heterostructure composed of a ferromagnetic film of
thickness $2d_{F}$ and a superconductor of thickness $d_{S}$, as shown in Fig.
3. We demonstrate the adjacent superconductors modulate strongly the radiated
electric field which explains the absence of the FMR shift in this
configuration [17, 49].

Figure 3: Radiated electric field of the F$|$S heterostructure.

###  IV.1 Full solution

Inside the ferromagnet, since $\nabla\times{\bf M}=0$ for uniform ${\bf M}$,
Eq. (8) has the solution $E_{z}(x)=E_{1}e^{ikx}+E_{1}^{\prime}e^{-ikx}$.
Inside the superconductor, according to Eqs. (1) and (4), the electric field
obeys

$\displaystyle\partial_{x}^{2}E_{z}+(\varepsilon_{0}\mu_{0}\omega^{2}-1/\lambda^{2})E_{z}=0,$
(19)

which has the solution
$E_{z}(x)=E_{2}e^{ik^{\prime}x}+E_{2}^{\prime}e^{-ik^{\prime}x}$, where
$k^{\prime}=\sqrt{\varepsilon_{0}\mu_{0}\omega^{2}-1/\lambda^{2}}$. Out of the
heterostructure, the electric fields $E_{3}e^{ikx}$ and $E_{4}e^{-ikx}$ are
radiated. These radiated electric fields are illustrated in Fig. 3.

The amplitudes $\\{E_{1},E_{1}^{\prime},E_{2},E_{2}^{\prime},E_{3},E_{4}\\}$
are governed by the boundary conditions, i.e., $E_{z}$ and $H_{y}$ are
continuous at interfaces. The continuous $E_{z}$ at interfaces requests

$\displaystyle
E_{1}e^{ikd_{F}}+E_{1}^{\prime}e^{-ikd_{F}}=E_{2}e^{ik^{\prime}d_{F}}+E_{2}^{\prime}e^{-ik^{\prime}d_{F}},$
$\displaystyle
E_{2}e^{ik^{\prime}(d_{F}+d_{S})}+E_{2}^{\prime}e^{-ik^{\prime}(d_{F}+d_{S})}=E_{3}e^{ik(d_{F}+d_{S})},$
$\displaystyle E_{1}e^{-ikd_{F}}+E_{1}^{\prime}e^{ikd_{F}}=E_{4}e^{ikd_{F}}.$
(20)

In the superconductors, $H_{y}=-1/(i\omega\mu_{0})\partial_{x}E_{z}$, while in
the ferromagnet, $H_{y}=-1/(i\omega\mu_{0})\partial_{x}E_{z}-M_{y}$, so the
continuous $H_{y}$ at interfaces leads to

$\displaystyle{k}(E_{1}e^{ikd_{F}}-E_{1}^{\prime}e^{-ikd_{F}})+\omega\mu_{0}M_{y}$
$\displaystyle=k^{\prime}(E_{2}e^{ik^{\prime}d_{F}}-E_{2}^{\prime}e^{-ik^{\prime}d_{F}}),$
$\displaystyle{k^{\prime}}(E_{2}e^{ik^{\prime}(d_{F}+d_{S})}-E_{2}^{\prime}e^{-ik^{\prime}(d_{F}+d_{S})})={k}E_{3}e^{ik(d_{F}+d_{S})},$
$\displaystyle
k(E_{1}e^{-ikd_{F}}-E_{1}^{\prime}e^{ikd_{F}})+\omega\mu_{0}M_{y}=-kE_{4}e^{ikd_{F}}.$
(21)

Combining Eqs. (20) and (21), we obtain all the amplitudes. In the
ferromagnetic insulator,

$\displaystyle E_{z}(-d_{F}<x<d_{F})$ $\displaystyle={\cal
R}E_{0}e^{-ik(x-d_{F})}-\frac{\omega\mu_{0}M_{y}}{2k}\left(e^{ik(x+d_{F})}-e^{-ik(x-d_{F})}\right),$
(22)

where the amplitude

$\displaystyle
E_{0}=-\frac{\omega\mu_{0}M_{y}}{2k}\left(e^{2ik(x+d_{F})}-1\right),$ (23)

and

$\displaystyle{\cal R}=\frac{e^{ik^{\prime}d_{S}}(k^{2}-k^{\prime
2})+e^{-ik^{\prime}d_{S}}(k^{\prime
2}-k^{2})}{e^{ik^{\prime}d_{S}}(k-k^{\prime})^{2}-e^{-ik^{\prime}d_{S}}(k+k^{\prime})^{2}}$
(24)

is the reflection coefficient of the electric field at the F$|$S interface.

When $d_{S}=0$, ${\cal R}=0$, the solution (22) recovers the single layer case
(12). On the other hand, even with a small $d_{S}\ll\lambda$, since
$|k|\ll|k^{\prime}|$ when $\omega\sim 100$ GHz, ${\cal R}\rightarrow-1$
implies the total reflection of the electric fields at the F$|$S interface
even with ultrathin conventional superconductor layer. As shown below, this
implies the absence of FMR shift in all the available experiments [17, 49].

Inside the superconductor,

$\displaystyle
E_{z}(d_{F}<x<d_{F}+d_{S})=\frac{2kE_{0}}{e^{ik^{\prime}d_{S}}(k-k^{\prime})^{2}-e^{-ik^{\prime}d_{S}}(k+k^{\prime})^{2}}$
$\displaystyle\times\left((k-k^{\prime})e^{-ik^{\prime}(x-d_{F}+d_{S})}-(k+k^{\prime})e^{ik^{\prime}(x-d_{F}-d_{S})}\right),$
(25)

which is indeed very weak since $|k|\ll|k^{\prime}|$. Out of the
heterostructure,

$\displaystyle
E_{z}(x>d_{F}+d_{S})=\frac{-4kk^{\prime}E_{0}e^{ik(x-d_{F}-d_{S})}}{e^{ik^{\prime}d_{S}}(k-k^{\prime})^{2}-e^{-ik^{\prime}d_{S}}(k+k^{\prime})^{2}},$
$\displaystyle E_{z}(x<-d_{F})={\cal R}E_{0}e^{-ik(x-d_{F})}$
$\displaystyle-\frac{\omega\mu_{0}M_{y}}{2k}(e^{-ik(x+d_{F})}-e^{-ik(x-d_{F})}).$
(26)

With low frequencies and near the heterostructure, $kx\rightarrow 0$,
$kd_{F}\rightarrow 0$, and $kd_{S}\rightarrow 0$, so the electric fields

$\displaystyle E_{z}(x)=\begin{cases}0,&x>d_{F}\\\
-i\omega\mu_{0}M_{y}(x-d_{F}),&-d_{F}<x<d_{F}\\\
2i\omega\mu_{0}M_{y}d_{F},&x<-d_{F}\end{cases},$ (27)

which is illustrated in Fig. 1(b) for a snapshot. The electric field vanishes
in the superconductor due to the total reflection with a $\pi$-phase shift
${\cal R}=-1$ that generates no supercurrent and thereby leads to no
modulation on the FMR.

###  IV.2 Quasi-static approximation

The full solution clearly shows the absence of electric fields at the
superconductor side of the S$|$F heterostructure, which can be well understood
within the quasi-static approximation $\nabla\times{\bf H}=0$ or ${\bf
J}_{s}$. Assuming $E_{z}(x=d_{F})=\tilde{E}_{0}$ at the F$|$S interface,
according to Eq. (6) the electric field in the adjacent superconductor

$\displaystyle
E_{z}(x)=\tilde{E}_{0}\dfrac{\cosh{((x-d_{S}-d_{F})/\lambda)}}{\cosh{(d_{S}/\lambda)}}$
(28)

drives the supercurrents. For a thin superconducting film of thickness
$O(\lambda)$, we are allowed to take an average of the supercurrent
$J_{s,z}=[J_{s,z}(x=d_{F})+J_{s,z}(x=d_{F}+d_{S})]/2$, and from the first
equation of Eq. (4)

$\displaystyle
J_{s,z}=\dfrac{i}{\mu_{0}\omega\lambda^{2}}\tilde{E}_{0}\dfrac{1+\cosh(d_{S}/\lambda)}{2\cosh(d_{S}/\lambda)}.$
(29)

The supercurrents generate the vector potential (7) and the Oersted magnetic
field according to ${H}_{y}=-\partial_{x}A_{z}/\mu_{0}$. Taking $k=0$ at low
frequencies in the Weyl identity (11), i.e., [7]

$\displaystyle\dfrac{1}{|{\bf r}-{\bf r}^{\prime}|}=\int
dk_{y}^{\prime}dk_{z}^{\prime}\dfrac{e^{ik_{y}^{\prime}(y-y^{\prime})+ik_{z}^{\prime}(z-z^{\prime})}e^{-\sqrt{k_{y}^{\prime
2}+k_{z}^{\prime 2}}|x-x^{\prime}|}}{2\pi\sqrt{k_{y}^{\prime 2}+k_{z}^{\prime
2}}},$ (30)

we obtain the Oersted magnetic field generated by the supercurrents

$\displaystyle
H_{s,y}(x)=\left\\{\begin{array}[]{cc}d_{S}{{J}}_{s,z}/2,&x>d_{F}+d_{S}\\\
-d_{S}{{J}}_{s,z}/2,&~{}x<d_{F}\end{array}\right..$ (33)

However, constant $H_{s,y}$ independent of $x$ should vanish out of the
heterostructure within the quasi-static approximation since a constant
magnetic field renders the radiated electric field divergent, which requests
$J_{s,z}=0$ when $d_{S}\neq 0$ and $E_{z}(x>d_{F})=0$. Since the electric
field is continuous at interfaces, $E_{z}(x=d_{F})=\tilde{E}_{0}=0$ and
according to Eq. (16) $E_{z}(x=-d_{F})=2id_{F}\omega\mu_{0}M_{y}$. These
simple calculations thereby capture precisely the key physics of the full
solution (27).

##  V S$|$F$|$S heterostructure

Further, we consider the S$|$F$|$S heterostructure as illustrated in Fig. 2
composed of the ferromagnetic insulator of thickness $2d_{F}$ and two adjacent
superconductor films of thickness $d_{S}$ and $d_{S}^{\prime}$, respectively.
In comparison to that of the S$|$F bilayer, the distribution of the electric
field in S$|$F$|$S heterostructure changes much due to its back-and-forth
reflection by the superconductors, as addressed in this section.

###  V.1 Full solution

Similar to the S$|$F heterostructure, inside the ferromagnet,
$E_{z}(x)=E_{1}e^{ikx}+E_{1}^{\prime}e^{-ikx}$; in the superconductor “1”,
$E_{z}(x)=E_{2}e^{ik^{\prime}x}+E_{2}^{\prime}e^{-ik^{\prime}x}$; and in the
superconductor “2”,
$E_{z}(x)=E_{3}e^{ik^{\prime}x}+E_{3}^{\prime}e^{-ik^{\prime}x}$. Out of the
heterostructure, the electric fields $E_{4}e^{ikx}$ and $E_{5}e^{-ikx}$ are
radiated. These electric fields are illustrated in Fig. 4.

Figure 4: Radiated electric field of the S$|$F$|$S heterostructure.

The amplitudes
$\\{E_{1},E_{1}^{\prime},E_{2},E_{2}^{\prime},E_{3},E_{3}^{\prime},E_{4},E_{5}\\}$
are governed by the boundary conditions. The continuous $E_{z}$ at interfaces
requests

$\displaystyle
E_{1}e^{ikd_{F}}+E_{1}^{\prime}e^{-ikd_{F}}=E_{2}e^{ik^{\prime}d_{F}}+E_{2}^{\prime}e^{-ik^{\prime}d_{F}},$
$\displaystyle
E_{1}e^{-ikd_{F}}+E_{1}^{\prime}e^{ikd_{F}}=E_{3}e^{-ik^{\prime}d_{F}}+E_{3}^{\prime}e^{ik^{\prime}d_{F}},$
$\displaystyle
E_{2}e^{ik^{\prime}(d_{F}+d_{S})}+E_{2}^{\prime}e^{-ik^{\prime}(d_{F}+d_{S})}=E_{4}e^{ik(d_{F}+d_{S})},$
$\displaystyle
E_{3}e^{-ik^{\prime}(d_{F}+d_{S}^{\prime})}+E_{3}^{\prime}e^{ik^{\prime}(d_{F}+d_{S}^{\prime})}=E_{5}e^{ik(d_{F}+d_{S}^{\prime})},$
(34)

and the continuous $H_{y}$ at interfaces leads to

$\displaystyle
k(E_{1}e^{ikd_{F}}-E_{1}^{\prime}e^{-ikd_{F}})+\omega\mu_{0}M_{y}=k^{\prime}(E_{2}e^{ik^{\prime}d_{F}}-E_{2}^{\prime}e^{-ik^{\prime}d_{F}}),$
$\displaystyle
k(E_{1}e^{-ikd_{F}}-E_{1}^{\prime}e^{ikd_{F}})+\omega\mu_{0}M_{y}=k^{\prime}(E_{3}e^{-ik^{\prime}d_{F}}-E_{3}^{\prime}e^{ik^{\prime}d_{F}}),$
$\displaystyle
k^{\prime}(E_{2}e^{ik^{\prime}(d_{F}+d_{S})}-E_{2}^{\prime}e^{-ik^{\prime}(d_{F}+d_{S})})=kE_{4}e^{ik(d_{F}+d_{S})},$
$\displaystyle
k^{\prime}(E_{3}e^{-ik^{\prime}(d_{F}+d_{S}^{\prime})}-E_{3}^{\prime}e^{ik^{\prime}(d_{F}+d_{S}^{\prime})})=-kE_{5}e^{ik(d_{F}+d_{S}^{\prime})}.$
(35)

Combining Eqs. (34) and (35), we obtain the electric-field distribution,
referring to Appendix A for the general solution. In particular, when
$d_{S}=d_{S}^{\prime}$, in the ferromagnetic film,

$\displaystyle E_{z}(-d_{F}<x<d_{F})$
$\displaystyle=\frac{-\omega\mu_{0}M_{y}\sinh{(ikx)}}{k\cosh{(ikd_{F})}-k^{\prime}f(u)\sinh{(ikd_{F})}},$
(36)

where $u=-[(k+k^{\prime})/(k-k^{\prime})]\exp(-2ik^{\prime}d_{S})$ and

$\displaystyle
f(u)=\frac{u-1}{u+1}=\frac{k^{\prime}\sinh{(ik^{\prime}d_{S})}-k\cosh{(ik^{\prime}d_{S})}}{k\sinh{(ik^{\prime}d_{S})}-k^{\prime}\cosh{(ik^{\prime}d_{S})}}.$
(37)

In the superconductor “1”,

$\displaystyle E_{z}(d_{F}<x<d_{F}+d_{S})$
$\displaystyle=\frac{-\omega\mu_{0}M_{y}(ue^{ik^{\prime}(x-d_{F})}+e^{-ik^{\prime}(x-d_{F})})}{k(1+u)\coth(ikd_{F})-k^{\prime}(u-1)},$
(38)

and in the superconductor “2”,

$\displaystyle E_{z}(-d_{F}-d_{S}<x<-d_{F})$
$\displaystyle=\frac{\omega\mu_{0}M_{y}(ue^{-ik^{\prime}(x+d_{F})}+e^{ik^{\prime}(x+d_{F})})}{k(1+u)\coth(ikd_{F})-k^{\prime}(u-1)}.$
(39)

They both exist, and $E_{z}(x=-d_{F})$ and $E_{z}(x=d_{F})$ are opposite. Out
of the heterostructure,

$\displaystyle E_{z}(x>d_{F}+d_{S})$
$\displaystyle=\frac{-\omega\mu_{0}M_{y}(ue^{ik^{\prime}d_{S}}+e^{-ik^{\prime}d_{S}})}{k(1+u)\coth(ikd_{F})-k^{\prime}(u-1)}e^{ikx},$
$\displaystyle E_{z}(x<-d_{F}-d_{S})$
$\displaystyle=\frac{\omega\mu_{0}M_{y}(ue^{ik^{\prime}d_{S}}+e^{-ik^{\prime}d_{S}})}{k(1+u)\coth(ikd_{F})-k^{\prime}(u-1)}e^{-ikx}.$
(40)

We illustrate in Fig. 5 the distribution of the electric fields ${\rm
Re}(E_{z}/(i\omega\mu_{0}M_{y}d_{F}))$ in the symmetric
$d_{S}^{\prime}=d_{S}=30$ nm and asymmetric $d_{S}^{\prime}=2d_{S}=60$ nm
S$|$F$|$S heterostructure, respectively, in the near-field limit. For NbN with
electron density $n_{e}=1.65\times 10^{28}$/m3 [58] and $T_{c}=6.5$ K, the
London penetration depth $\lambda(T=0)=42.0$ nm and $\lambda(T=0.5T_{c})=43.4$
nm. The fields are opposite at the two superconductors in the symmetric
heterostructure but are skewed when $d_{S}\neq d_{S}^{\prime}$. These fields
carrying energy are radiated out in the far zone [56]. When the
superconductors are sufficiently thick $\\{d_{S},d_{S}^{\prime}\\}\gg\lambda$,
these electric fields are confined between them, which corresponds to an
excellent waveguide with small size [36].

Figure 5: Distribution of electric fields in symmetric
$d_{S}=d_{S}^{\prime}=30$ nm [(a)] and asymmetric $d_{S}^{\prime}=2d_{S}=60$
nm [(b)] S$|$F$|$S heterostructure. The thickness of the ferromagnetic film
2$d_{F}=60$ nm and London’s penetration depth $\lambda=43.4$ nm.

###  V.2 Quasi-static approximation

As justified, the quasi-static approximation $\nabla\times{\bf H}=0$ or ${\bf
J}_{s}$ is allowed when solving the electric fields near the heterostructure
[56]. In the FMR case, the radiated electric field is uniform in the $y$-$z$
plane, so from $\nabla\times{\bf E}=i\omega{\bf B}$, the $x$-component
$B_{x}=H_{d,x}+M_{x}=0$ generates no electric field outside the magnet. On the
other hand, in the linear response regime $B_{z}=\mu_{0}(H_{0}+M_{z})$ is
static, so only $B_{y}=\mu_{0}M_{y}$ in the magnet radiates the time-dependent
electric field according to
$-\partial_{x}E_{z}=i\omega\mu_{0}(M_{y}+H_{s,y})$. Integrating along $x$
across the ferromagnet yields the net electric field at the interfaces obeying

$\displaystyle
E_{z}(x=d_{F})-E_{z}(x=-d_{F})=-2d_{F}i\omega\mu_{0}(M_{y}+H_{s,y}).$ (41)

Out of the heterostructure, from the $z$-component of $\nabla\times{\bf H}=0$,
$H_{y}|_{\rm outside}$ is a constant, which can be proved to vanish as in Sec.
IV.2.

In the quasi-static approximation, the electric field in the superconductors
$``1"$ and $``2"$ obeys Eq. (6). From the boundary conditions with continuous
$E_{z}$ and $H_{y}$ at interfaces and $H_{y}|_{\text{outside}}=0$, the
electric field in the superconductors reads

$\displaystyle E_{z}(d_{F}<x<d_{F}+d_{S})$
$\displaystyle=E_{z}(x=d_{F})\dfrac{\cosh((x-d_{S}-d_{F})/\lambda)}{\cosh(d_{S}/\lambda)},$
$\displaystyle E_{z}(-d_{F}-d_{S}<x<-d_{F})$
$\displaystyle=E_{z}(x=-d_{F})\dfrac{\cosh((x+d_{S}^{\prime}+d_{F})/\lambda)}{\cosh(d_{S}^{\prime}/\lambda)},$
(42)

which drive the supercurrents in the adjacent superconductors. For thin
superconducting films of thickness $O({\lambda})$, we are allowed to take an
average of the supercurrents ${\bf J}^{(1)}_{s}=\left[{\bf
J}_{s}(x=d_{F})+{\bf J}_{s}(x=d_{F}+d_{S})\right]/2$ and ${\bf
J}^{(2)}_{s}=\left[{\bf J}_{s}(x=-d_{F})+{\bf
J}_{s}(x=-d_{F}-d_{S}^{\prime})\right]/2$, i.e.,

$\displaystyle{{J}}_{s,z}^{(1)}=\dfrac{i}{\omega\mu_{0}\lambda^{2}}E_{z}(x=d_{F})\dfrac{1+\cosh(d_{S}/\lambda)}{2\cosh(d_{S}/\lambda)},$
$\displaystyle{{J}}^{(2)}_{s,z}=\dfrac{i}{\omega\mu_{0}\lambda^{2}}E_{z}(x=-d_{F})\dfrac{1+\cosh(d_{S}^{\prime}/\lambda)}{2\cosh(d_{S}^{\prime}/\lambda)}.$
(43)

The supercurrents generate the vector potential (7) and the Oersted magnetic
field according to ${H}_{s,y}=-\partial_{x}A_{z}/\mu_{0}$. Using the Weyl
identity (30) we obtain

$\displaystyle{H_{s,y}}(x)=\left\\{\begin{array}[]{cc}\left(d_{S}{{J}}^{(1)}_{s,z}+d_{S}^{\prime}{{J}}^{(2)}_{s,z}\right)/2,&x>d_{F}+d_{S}\\\
\left(-d_{S}{{J}}^{(1)}_{s,z}+d_{S}^{\prime}{{J}}^{(2)}_{s,z}\right)/2,&~{}-d_{F}<x<d_{F}\\\
\left(-d_{S}{{J}}^{(1)}_{s,z}-d_{S}^{\prime}{{J}}^{(2)}_{s,z}\right)/2,&~{}~{}x<-d_{F}-d_{S}^{\prime}\end{array}\right..$
(47)

$H_{s,y}|_{\text{outside}}=0$ requests

$\displaystyle d_{S}{{J}}^{(1)}_{s,z}+d_{S}^{\prime}{{J}}^{(2)}_{s,z}=0,$ (48)

so the Oersted magnetic field inside the ferromagnetic slab is reduced to

$\displaystyle{H}_{s,y}(-d_{F}<x<d_{F})=d_{S}^{\prime}{J}_{s,z}^{(2)}=-d_{S}J_{s,z}^{(1)}.$
(49)

Thereby, when $d_{S}=d_{S}^{\prime}$, the supercurrents are opposite in the
two superconductors. When $d_{S}^{\prime}\rightarrow 0$, $H_{s,y}$ vanishes in
the magnet.

Substituting Eqs. (43) and (41) into (48), we obtain the electric field at the
surface of the ferromagnetic film:

$\displaystyle E_{z}(x=-d_{F})=i\mu_{0}\omega
d_{S}d_{F}(M_{y}+H_{s,y})\dfrac{\cosh(d_{S}/\lambda)+1}{\cosh(d_{S}/\lambda)}$
$\displaystyle\times\left(\dfrac{d_{S}(\cosh(d_{S}/\lambda)+1)}{2\cosh(d_{S}/\lambda)}+\dfrac{d_{S}^{\prime}(\cosh(d_{S}^{\prime}/\lambda)+1)}{2\cosh(d_{S}^{\prime}/\lambda)}\right)^{-1}.$
(50)

Substituting it into Eq. (49), the Oersted magnetic field in the ferromagnetic
film

$\displaystyle{H}_{s,y}(-d_{F}<x<d_{F})$
$\displaystyle=-M_{y}\dfrac{d_{F}d_{S}^{\prime}d_{S}G(d_{S},d_{S}^{\prime},\lambda)}{\lambda^{2}+d_{F}d_{S}^{\prime}d_{S}G(d_{S},d_{S}^{\prime},\lambda)},$
(51)

where

$\displaystyle
G(d_{S},d_{S}^{\prime},\lambda)=\dfrac{(\cosh(d_{S}/\lambda)+1)}{\cosh(d_{S}/\lambda)}\dfrac{(\cosh(d_{S}^{\prime}/\lambda)+1)}{\cosh(d_{S}^{\prime}/\lambda)}$
$\displaystyle\times\left(\dfrac{d_{S}(\cosh(d_{S}/\lambda)+1)}{\cosh(d_{S}/\lambda)}+\dfrac{d_{S}^{\prime}(\cosh(d_{S}^{\prime}/\lambda)+1)}{\cosh(d_{S}^{\prime}/\lambda)}\right)^{-1}.$
(52)

These results capture precisely the key physics of the full solution and are
convenient for the calculation of the interaction between Kittel magnon and
Cooper pairs.

###  V.3 Ultrastrong interaction between Kittel magnon and Cooper pairs

Above we address that the dynamics of magnetization ${\bf M}$ generates
$H_{s,y}$ via the backaction of superconductors, which, in turn, drives ${\bf
M}$ in the ferromagnet, imposing a self-consistent problem that is solved by
combining the Landau-Lifshitz and Maxwell’s equations.

In the linear regime, the Landau-Lifshitz equation

$\displaystyle-i\omega M_{x}+\mu_{0}\gamma M_{y}H_{0}=\mu_{0}\gamma
M_{0}H_{s,y},$ $\displaystyle i\omega M_{y}+\mu_{0}\gamma
M_{x}H_{0}=\mu_{0}\gamma M_{0}H_{d,x}.$ (53)

Substituting $B_{x}=M_{x}+H_{d,x}=0$ into Eq. (53), $M_{y}$ relates to
$H_{s,y}$ via

$\displaystyle
M_{y}=\dfrac{\mu_{0}^{2}\gamma^{2}M_{0}(H_{0}+M_{0})}{\mu_{0}^{2}\gamma^{2}H_{0}(H_{0}+M_{0})-\omega^{2}}H_{s,y}.$
(54)

When $d_{S}^{\prime}\rightarrow 0$, $H_{s,y}=0$ according to Eq. (51), and the
FMR frequency recovers the Kittel formula $\tilde{\omega}_{\rm
K}=\mu_{0}\gamma\sqrt{H_{0}(H_{0}+M_{0})}$ [52]. With finite $d_{S}$ and
$d_{S}^{\prime}$, the FMR frequency is self-consistently solved via combining
Eqs. (51) and (54), leading to the modified FMR frequency

$\displaystyle\omega_{\rm K}=\mu_{0}\gamma$
$\displaystyle\times\sqrt{\frac{\lambda^{2}H_{0}(H_{0}+M_{0})+d_{S}d_{S}^{\prime}d_{F}G(d_{S},d_{S}^{\prime},\lambda)(H_{0}+M_{0})^{2}}{d_{S}d_{S}^{\prime}d_{F}G(d_{S},d_{S}^{\prime},\lambda)+\lambda^{2}}}.$
(55)

In particular, when $d_{S}=d_{S}^{\prime}$,

$\displaystyle\omega_{\rm K}$
$\displaystyle={\mu_{0}\gamma}\left(\dfrac{2\lambda^{2}\cosh{(d_{S}/\lambda)}H_{0}(H_{0}+M_{0})}{d_{S}d_{F}\left(\cosh{(d_{S}/\lambda)}+1\right)+2\lambda^{2}\cosh{(d_{S}/\lambda)}}\right.$
$\displaystyle\left.+\dfrac{d_{S}d_{F}(\cosh{(d_{S}/\lambda)}+1)(H_{0}+M_{0})^{2}}{d_{S}d_{F}\left(\cosh{(d_{S}/\lambda)}+1\right)+2\lambda^{2}\cosh{(d_{S}/\lambda)}}\right)^{1/2}.$
(56)

Approaching $T_{c}$, $\lambda\rightarrow\infty$,
$\cosh(d_{S}/\lambda)\rightarrow 1$, so the FMR frequency (56) recovers the
Kittel formula $\omega_{\rm K}\rightarrow\tilde{\omega}_{\rm K}$; otherwise
$T<T_{c}$, it is shifted.

To show the FMR shift, we assume an oscillating magnetic field
$\tilde{H}e^{-i\omega_{0}t}\hat{\bf y}$ of frequency $\omega_{0}$ applied
along the $\bf\hat{y}$-direction (the associated microwave electric field is
along the normal $\hat{\bf x}$-direction). The wavelength of this microwave is
much larger than the thickness of the heterostructure, so it can be treated as
uniform across the heterostructure thickness. It can penetrate the
superconductor easily when $\\{d_{S},d_{S}^{\prime}\\}\sim\lambda$. With the
wave vector$\parallel\hat{\bf z}$ parallel to the film, it only excites $\bf
M$ in the ferromagnet but does not drive the superconductor.

Including the external pump field $\tilde{H}e^{-i\omega_{0}t}\hat{\bf y}$ and
$H_{s,y}\hat{\bf y}$ from the superconductor (51) into ${\bf H}_{\rm eff}$ and
incorporating the Gilbert damping $\alpha_{G}$, the linearized Landau-
Lifshitz-Gilbert equation reads

$\displaystyle-i\omega_{0}M_{x}+\mu_{0}\gamma M_{y}H_{0}$
$\displaystyle=\mu_{0}\gamma
M_{0}H_{\text{eff},y}+i\alpha_{G}\omega_{0}M_{y},$ $\displaystyle\mu_{0}\gamma
H_{0}M_{x}+i\omega_{0}M_{y}$ $\displaystyle=\mu_{0}\gamma
M_{0}H_{\text{eff},x}+i\alpha_{G}\omega_{0}M_{x},$ (57)

from which we solve with $\alpha_{G}\ll 1$

$\displaystyle M_{y}$
$\displaystyle=\dfrac{\mu_{0}^{2}\gamma^{2}M_{0}(H_{0}+M_{0})}{\omega_{\rm
K}^{2}-\omega_{0}^{2}-i\Gamma}\tilde{H},$ $\displaystyle M_{x}$
$\displaystyle=-iM_{y}\left[\dfrac{\omega_{0}}{\mu_{0}\gamma(H_{0}+M_{0})}+\dfrac{i\alpha_{G}\omega_{0}^{2}}{(\mu_{0}\gamma(H_{0}+M_{0}))^{2}}\right],$
(58)

where

$\displaystyle\Gamma$
$\displaystyle=\dfrac{\alpha_{G}\omega_{0}(\mu_{0}^{2}\gamma^{2}(H_{0}+M_{0})^{2}+\omega_{0}^{2})}{\mu_{0}\gamma(H_{0}+M_{0})}.$
(59)

The Oersted field in the thin magnetic film induced by the supercurrent

$\displaystyle H_{s,y}=$
$\displaystyle-\tilde{H}\dfrac{\mu_{0}^{2}\gamma^{2}M_{0}(H_{0}+M_{0})}{\omega_{\rm
K}^{2}-\omega_{0}^{2}-i\Gamma}\dfrac{d_{F}d_{S}d_{S}^{\prime}G(d_{S},d_{S}^{\prime},\lambda)}{\lambda^{2}+d_{S}d_{F}d_{S}^{\prime}G(d_{S},d_{S}^{\prime},\lambda)}.$
(60)

The average supercurrent density in (one of) the thin superconductor

$\displaystyle
J_{s}^{(1)}=\tilde{H}\dfrac{\mu_{0}^{2}\gamma^{2}M_{0}(H_{0}+M_{0})}{\omega_{\rm
K}^{2}-\omega_{0}^{2}-i\Gamma}\dfrac{d_{F}d_{S}^{\prime}G(d_{S},d_{S}^{\prime},\lambda)}{\lambda^{2}+d_{S}d_{F}d_{S}^{\prime}G(d_{S},d_{S}^{\prime},\lambda)}.$
(61)

The average electric field $E_{z}$ in (one of) the superconductors

$\displaystyle
E^{(1)}_{z}=-i\tilde{H}\dfrac{\mu_{0}^{3}\gamma^{2}M_{0}(H_{0}+M_{0})}{\omega_{\rm
K}^{2}-\omega_{0}^{2}-i\Gamma}\dfrac{\omega_{0}\lambda^{2}d_{F}d_{S}^{\prime}G(d_{S},d_{S}^{\prime},\lambda)}{\lambda^{2}+d_{S}d_{F}d_{S}^{\prime}G(d_{S},d_{S}^{\prime},\lambda)}.$
(62)

Here we illustrate the numerical results considering a yttrium iron garnet
(YIG) film of thickness $2d_{F}=60$ nm sandwiched by two NbN superconductors
of thickness $d_{S}=d_{S}^{\prime}=30$ nm. Insulating EuS thin magnetic film
[59, 60] is also a possible candidate to test our prediction. For YIG,
$\mu_{0}M_{0}=0.2$ T and $\alpha_{G}=5\times 10^{-4}$ [61, 62]. We use
$\lambda(T=0.5T_{c})=43.4$ nm for NbN [58] . We take the bias field
$\mu_{0}H_{0}=0.05$ T and the excitation field $\mu_{0}\tilde{H}=0.01$ mT.
Figure 6 shows the radiated electric field in (one of) the superconductors and
the excited amplitudes of $\bf M$ as a function of the excitation frequency
$\omega_{0}$. The frequency shift is $2\pi\times 1.6$ GHz, comparable to half
of the bare FMR frequency $\tilde{\omega}_{\rm K}=2\pi\times 3.2$ GHz,
corresponding to the decrease of the resonant magnetic field as large as 350
mT. This demonstrates the potential to achieve ultrastrong interaction between
magnons and Cooper pairs even with magnetic insulators.

Figure 6: FMR spectra with the excitation field $\mu_{0}\tilde{H}=0.01$ mT.
(a) plots the excited electric field amplitude in (one of) the superconductors
in the symmetric S$|$F$|$S heterostructure. The amplitude of the resonance
electric field $E_{z}\sim 9$ V/m. (b) is the excited amplitudes of the
magnetization $M_{x,y}$ with and without two adjacent superconductors. The
frequency shift is as large as $2\pi\times 1.6$ GHz$\sim\tilde{\omega}_{\rm
K}/2$.

##  VI Conclusion and discussion

Magnetic insulators are ideal candidates for long-range spin transport [1, 2,
3, 4, 5, 6, 7], strong coupling between magnons and microwaves [32], and
quantum information processing [29, 30, 33], gating which by superconductors
may bring new control dimensions. In comparison to metallic magnets, the
mutual proximity effect may differ between magnetic insulators and
superconductors, which may be helpful to distinguish different competitive
mechanisms [26] in the future studies. Our model system differs from the
metallic ferromagnets since there are no electric currents flowing in the
insulators that, if large, may affect the field distribution via radiation.
Our theory may apply to the antiferromagnet as well, in which case we need to
replace the characterized frequency with terahertz ones.

In conclusion, we analyze the interaction between the Kittel magnons in
insulating magnetic film and Cooper pairs in superconductors mediated by the
radiated electric fields from the magnetization dynamics. Via highlighting the
role of the total reflection of the electric fields at the ferromagnet-
superconductor interface that are solved beyond the quasi-static
approximation, we provide a comprehensive understanding of the absence of the
FMR shift in the F$|$S heterostructure and predict its existence in the
S$|$F$|$S heterostructure with the Meissner screening. The coupling between
magnons and Cooper pairs is ultrastrong with the frequency shift achieving
tens of percent of the bare FMR frequency, which may bring superior advantage
in information processing in on-chip magnonics and quantum magnonics.

###### Acknowledgements.

We gratefully acknowledge Prof. Guang Yang and Prof. Lihui Bai for many
inspiring discussions. This work is financially supported by the National
Natural Science Foundation of China, and the startup grant of Huazhong
University of Science and Technology (Grants No. 3004012185 and 3004012198).

##  Appendix A General solution of $E_{z}$ in S$|$F$|$S heterostructure

Here we list the general solution of $E_{z}(x)$ in the S$|$F$|$S
heterostructure when $d_{S}\neq d_{S}^{\prime}$ in Fig. 4. Inside the
ferromagnet,

$\displaystyle E_{z}(-d_{F}<x<d_{F})$
$\displaystyle=\frac{-\omega\mu_{0}M_{y}(Ge^{ikx}+e^{-ikx})}{k(Ge^{ikd_{F}}-e^{-ikd_{F}})-k^{\prime}f(u)(Ge^{ikd_{F}}+e^{-ikd_{F}})},$
(63)

where

$\displaystyle
G=-\frac{-2k\sinh(ikd_{F})+k^{\prime}(f(u)e^{-ikd_{F}}+f(u^{\prime})e^{ikd_{F}})}{-2k\sinh(ikd_{F})+k^{\prime}(f(u)e^{ikd_{F}}+f(u^{\prime})e^{-ikd_{F}})},$
(64)

and
$u^{\prime}=-[(k+k^{\prime})/(k-k^{\prime})]\exp(-2ik^{\prime}d_{S}^{\prime})$.
In the superconductor “1”,

$\displaystyle
E_{z}(d_{F}<x<d_{F}+d_{S})=\frac{ue^{ik^{\prime}(x-d_{F})}+e^{-ik^{\prime}(x-d_{F})}}{1+u}$
$\displaystyle\times\frac{-\omega\mu_{0}M_{y}(Ge^{ikd_{F}}+e^{-ikd_{F}})}{k(Ge^{ikd_{F}}-e^{-ikd_{F}})-k^{\prime}f(u)(Ge^{ikd_{F}}+e^{-ikd_{F}})}.$
(65)

In the superconductor “2”,

$\displaystyle
E_{z}(-d_{F}-d_{S}^{\prime}<x<-d_{F})=\frac{e^{ik^{\prime}(x+d_{F})}+u^{\prime}e^{-ik^{\prime}(x+d_{F})}}{1+u^{\prime}}$
$\displaystyle\times\frac{-\omega\mu_{0}M_{y}(Ge^{-ikd_{F}}+e^{ikd_{F}})}{k(Ge^{ikd_{F}}-e^{-ikd_{F}})-k^{\prime}f(u)(Ge^{ikd_{F}}+e^{-ikd_{F}})}.$
(66)

Out of the heterostructure,

$\displaystyle
E_{z}(x>d_{F}+d_{S})=\frac{ue^{ik^{\prime}d_{S}}+e^{-ik^{\prime}d_{S}}}{1+u}$
$\displaystyle\times\frac{-\omega\mu_{0}M_{y}(Ge^{ikd_{F}}+e^{-ikd_{F}})e^{ik(x-d_{F}-d_{S})}}{k(Ge^{ikd_{F}}-e^{-ikd_{F}})-k^{\prime}f(u)(Ge^{ikd_{F}}+e^{-ikd_{F}})},$
$\displaystyle
E_{z}(x<-d_{F}-d_{S}^{\prime})=\frac{e^{-ik^{\prime}d_{S}^{\prime}}+u^{\prime}e^{ik^{\prime}d_{S}^{\prime}}}{1+u^{\prime}}$
$\displaystyle\times\frac{-\omega\mu_{0}M_{y}(Ge^{-ikd_{F}}+e^{ikd_{F}})e^{-ik(x+d_{F}+d_{S})}}{k(Ge^{ikd_{F}}-e^{-ikd_{F}})-k^{\prime}f(u)(Ge^{ikd_{F}}+e^{-ikd_{F}})}.$
(67)

## References

  * [1] B. Lenk, H. Ulrichs, F. Garbs, and M. Münzenberg, The building blocks of magnonics, Phys. Rep. 507, 107 (2011).
  * [2] A. V. Chumak, V. I. Vasyuchka, A. A. Serga, and B. Hillebrands, Magnon spintronics, Nat. Phys. 11, 453 (2015).
  * [3] D. Grundler, Nanomagnonics around the corner, Nat. Nanotechnol. 11, 407 (2016).
  * [4] V. E. Demidov, S. Urazhdin, G. de Loubens, O. Klein, V. Cros, A. Anane, and S. O. Demokritov, Magnetization oscillations and waves driven by pure spin currents, Phys. Rep. 673, 1 (2017).
  * [5] A. Brataas, B. van Wees, O. Klein, G. de Loubens, and M. Viret, Spin Insulatronics, Phys. Rep. 885, 1 (2020).
  * [6] Barman et al., The 2021 Magnonics Roadmap, J. Phys. Condens. Matter 33, 413001 (2021).
  * [7] T. Yu, Z. C. Luo, and G. E. W. Bauer, Chirality as generalized spin–orbit interaction in spintronics, Phys. Rep. 1009, 1 (2023).
  * [8] O. V. Dobrovolskity, R. Sachser, T. Brächer, T. Böttcher, V. V. Kruglyak, R. V. Vovk, V. A. Shklovskij, M. Huth, B. Hillebrands, and A. V. Chumak, Magnon-fluxon interaction in a ferromagnet/superconductor heterostructure, Nat. Phys. 15, 477 (2019).
  * [9] Y. Yao, R. Cai, T. Yu, Y. Ma, W. Xing, Y. Ji, X.- C. Xie, S.-H. Yang, and W. Han, Giant oscillatory Gilbert damping in superconductor/ferromagnet/superconductor junctions, Sci. Adv. 7, eabh3686 (2021).
  * [10] L. G. Johnsen, H. T. Simensen, A. Brataas, and J. Linder, Magnon Spin Current Induced by Triplet Cooper Pair Supercurrents, Phys. Rev. Lett. 127, 207001 (2021).
  * [11] T. Yu and G. E. W. Bauer, Efficient Gating of Magnons by Proximity Superconductors, Phys. Rev. Lett. 129, 117201 (2022).
  * [12] M. A. Kuznetsov and A. A. Fraerman, Temperature-sensitive spin-wave nonreciprocity induced by interlayer dipolar coupling in ferromagnet/paramagnet and ferromagnet/superconductor hybrid systems, Phys. Rev. B 105, 214401 (2022).
  * [13] M. Silaev, Anderson-Higgs Mass of Magnons in Superconductor-Ferromagnet-Superconductor Systems, Phys. Rev. Appl. 18, L061004 (2022).
  * [14] I. V. Bobkova, A. M. Bobkov, A. Kamra, and W. Belzig, Magnon-cooparons in magnet-superconductor hybrids, Commun. Mater. 3, 95 (2022).
  * [15] A. M. Bobkov, S. A. Sorokin, and I. V. Bobkova, Phys. Rev. B 107, 174521 (2023).
  * [16] A. S. Ianovskaia, A. M. Bobkov, and I. V. Bobkova, Magnon influence on the superconducting DOS in FI/S bilayers, arXiv:2307. 03954.
  * [17] M. Borst, P. H. Vree, A. Lowther, A. Teepe, S. Kurdi, I. Bertelli, B. G. Simon, Y. M. Blanter, and T. van der Sar, Observation and control of hybrid spin-wave–Meissner-current transport modes, arXiv:2307.07581.
  * [18] A. T. G. Janssønn, H. T. Simensen, A. Kamra, A. Brataas, and S. H. Jacobsen, Macroscale nonlocal transfer of superconducting signatures to a ferromagnet in a cavity, Phys. Rev. B 102, 180506(R) (2020).
  * [19] I. A. Golovchanskiy, N. N. Abramov, V. S. Stolyarov, M. Weides, V. V. Ryazanov, A. A. Golubov, A. V. Ustinov, and M. Y. Kupriyanov, Ultrastrong photon-to-magnon coupling in multilayered heterostructures involving superconducting coherence via ferromagnetic layers, Sci. Adv. 7, eabe8638 (2021).
  * [20] I. A. Golovchanskiy, N. N. Abramov, V. S. Stolyarov, A. A. Golubov, M. Y. Kupriyanov, V. V. Ryazanov, and A. V. Ustinov, Approaching Deep-Strong On-Chip Photon-To-Magnon Coupling, Phys. Rev. Appl. 16, 034029 (2021).
  * [21] M. Silaev, Ultrastrong magnon-photon coupling, squeezed vacuum, and entanglement in superconductor/ferromagnet nanostructures, Phys. Rev. B 107, L180503 (2023).
  * [22] A. T. G. Janssønn, H. G. Hugdal, A. Brataas, and S. H. Jacobsen, Cavity-mediated superconductor–ferromagnetic-insulator coupling, Phys. Rev. B 107, 035147 (2023).
  * [23] A. Ghirri, C. Bonizzoni, M. Maksutoglu, A. Mercurio, O. D. Stefano, S. Savasta, and M. Affronte, Ultrastrong magnon-photon coupling achieved by magnetic films in contact with superconducting resonators, Phys. Rev. Appl. (2023).
  * [24] M. Eschrig, Spin-polarized supercurrents for spintronics: a review of current progress, Rep. Prog. Phys. 78, 104501 (2015).
  * [25] J. Linder and J. W. A. Robinson, Superconducting spintronics, Nat. Phys. 11, 307 (2015).
  * [26] F. S. Bergeret, M. Silaev, P. Virtanen, and T. T. Heikkilä, Colloquium: Nonequilibrium effects in superconductors with a spin-splitting field, Rev. Mod. Phys. 90, 041001 (2018).
  * [27] J. Linder and A. V. Balatsky, Odd-frequency superconductivity, Rev. Mod. Phys. 91, 045005 (2019).
  * [28] M. Amundsen, J. Linder, J. W. A. Robinson, I. Žutić, and N. Banerjee, Colloquium: Spin-orbit effects in superconducting hybrid structures, arXiv: 2210.03549.
  * [29] Y. Tabuchi, S. Ishino, A. Noguchi, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura, Coherent coupling between a ferromagnetic magnon and a superconducting qubit, Science 349, 405 (2015).
  * [30] D. L. Quirion, S. P. Wolski, Y. Tabuchi, S. Kono, K. Usami, and Y. Nakamura, Entanglement-based single-shot detection of a single magnon with a superconducting qubit, Science  367, 425 (2020).
  * [31] T. Yu, M. Claassen, D. M. Kennes, and M. A. Sentef, Optical manipulation of domains in chiral topological superconductors, Phys. Rev. Research 3, 013253 (2021).
  * [32] B. Z. Rameshti, S. V. Kusminskiy, J. A. Haigh, K. Usami, D. Lachance-Quirion, Y. Nakamura, C. -M. Hu, H. X. Tang, G. E. W. Bauer, and Y. M. Blanter, Cavity magnonics, Phys. Rep. 979, 1 (2022).
  * [33] D. Xu, X. -K. Gu, H. -K. Li, Y. -C. Weng, Y. -P. Wang, J. Li, H. Wang, S. -Y. Zhu, and J. Q. You, Quantum Control of a Single Magnon in a Macroscopic Spin System, Phys. Rev. Lett. 130, 193603 (2023).
  * [34] K. M. D. Hals, M. Schecter, and M. S. Rudner, Composite Topological Excitations in Ferromagnet-Superconductor Heterostructures, Phys. Rev. Lett. 117, 017001 (2016).
  * [35] A. F. Kockum, A. Miranowicz, S. D. Liberato, S. Savasta, and F. Nori, Ultrastrong coupling between light and matter, Nat. Rev. Phys. 1, 19 (2019).
  * [36] J. C. Swihart, Field Solution for a Thin‐Film Superconducting Strip Transmission Line, J. Appl. Phys. 32, 461 (1961).
  * [37] I. A. Golovchanskiy, N. N. Abramov, V. S. Stolyarov, V. V. Bolginov, V. V. Ryazanov, A. A. Golubov, and A. V. Ustinov, Ferromagnet/Superconductor Hybridization for Magnonic Applications, Adv. Funct. Mater. 28, 1802375 (2018).
  * [38] I. A. Golovchanskiy, N. N. Abramov, V. S. Stolyarov, V. V. Ryazanov, A. A. Golubov, and A. V. Ustinov, Modified dispersion law for spin waves coupled to a superconductor, J. Appl. Phys. 124, 233903 (2018).
  * [39] I. A. Golovchanskiy, N. N. Abramov, V. S. Stolyarov, P. S. Dzhumaev, O. V. Emelyanova, A. A. Golubov, V. V. Ryazanov, and A. V. Ustinov, Ferromagnet/Superconductor Hybrid Magnonic Metamaterials, Adv. Sci. 6, 1900435 (2019).
  * [40] I. A. Golovchanskiy, N. N. Abramov, V. S. Stolyarov, A. A. Golubov, V. V. Ryazanov, and A. V. Ustinov, Nonlinear spin waves in ferromagnetic/superconductor hybrids, J. Appl. Phys. 127, 093903 (2020).
  * [41] Seshadri, Surface magnetostatic modes of a ferrite slab, Proc. IEEE 58, 506 (1970).
  * [42] C. Bayer, J. Jorzick, B. Hillebrands, S. O. Demokritov, R. Kouba, R. Bozinoski, A. N. Slavin, K. Y. Guslienko, D. V. Berkov, N. L. Gorn, and M. P. Kostylev, Spin-wave excitations in finite rectangular elements of ${\mathrm{Ni}}_{80}{\mathrm{Fe}}_{20}$, Phys. Rev. B 72, 064427 (2005).
  * [43] T. Yu, C. P. Liu, H. M. Yu, Y. M. Blanter, and G. E. W. Bauer, Chiral excitation of spin waves in ferromagnetic films by magnetic nanowire gratings, Phys. Rev. B 99, 134424 (2019).
  * [44] T. Yu, Y. M. Blanter, and G. E. W. Bauer, Chiral Pumping of Spin Waves, Phys. Rev. Lett. 123, 247202 (2019).
  * [45] T. Yu, J. Zou, B. Zeng, J. W. Rao, and K. Xia, Non-Hermitian Topological Magnonics, arXiv:2306.04348.
  * [46] A. V. Chumak, A. A. Serga, and B. Hillebrands, Nat. Commun. 5, 4700 (2014).
  * [47] T. Yu, H. Wang, M. A. Sentef, H. Yu, and G. E. W. Bauer, Magnon trap by chiral spin pumping, Phys. Rev. B 102, 054429 (2020).
  * [48] O. A. Santos and B. J. van Wees, Magnon confinement in an all-on-chip YIG cavity resonator using hybrid YIG/Py magnon barriers, arXiv:2306.14029.
  * [49] L. -L. Li, Y. -L. Zhao, X. -X. Zhang, and Y. Sun, Possible evidence for spin-transfer torque induced by spin-triplet supercurrents, Chin. Phys. Lett. 35, 077401 (2018).
  * [50] I. A. Golovchanskiy, N. N. Abramov, V. S. Stolyarov, V. I. Chichkov, M. Silaev, I. V. Shchetinin, A. A. Golubov, V. V. Ryazanov, A. V. Ustinov, and M. Y. Kupriyanov, Magnetization Dynamics in Proximity-Coupled Superconductor-Ferromagnet-Superconductor Multilayers, Phys. Rev. Appl. 14, 024086 (2020).
  * [51] K. -R. Jeon, C. Ciccarelli, H. Kurebayashi, L. F. Cohen, X. Montiel, M. Eschrig, T. Wagner, S. Komori, A. Srivastava, J. W. A. Robinson, and M. G. Blamire, Effect of Meissner Screening and Trapped Magnetic Flux on Magnetization Dynamics in Thick $\mathrm{Nb}/{\mathrm{Ni}}_{80}{\mathrm{Fe}}_{20}/\mathrm{Nb}$ Trilayers, Phys. Rev. Appl. 11, 014061 (2019).
  * [52] C. Kittel, On the Theory of Ferromagnetic Resonance Absorption, Phys. Rev. 73, 155 (1948).
  * [53] S. V. Mironov and A. I. Buzdin, Giant demagnetization effects induced by superconducting films, Appl. Phys. Lett. 119, 102601 (2021).
  * [54] S. M. Rezende, Fundamentals of Magnonics (Springer, Cham, 2020).
  * [55] L. D. Landau and E. M. Lifshitz, Electrodynamics of Continuous Media, 2nd ed. (Butterworth-Heinenann, Oxford, U.K., 1984).
  * [56] J. D. Jackson, Classical Electrodynamics (Wiley, New York, 1998).
  * [57] J. R. Schrieffer, Theory of Superconductivity, (W. A. Benjimin, New York, 1964).
  * [58] S. P. Chockalingam, M. Chand, J. Jesudasan, V. Tripathi, and P. Raychaudhuri, Superconducting properties and Hall effect of epitaxial NbN thin films, Phys. Rev. B 77, 214503 (2008).
  * [59] O. W. Dietrich, A. J. Henderson, Jr., and H. Meyer, Spin-wave analysis of specific heat and magnetization in Euo and Eus, Phys. Rev. B 12, 2844 (1975).
  * [60] Y. Hou, F. Nichele, H. Chi, A. Lodesani, Y. Wu, M. F. Ritter, D. Z. Haxell, M. Davydova, S. Ilić, O. Glezakou-Elbert, A. Varambally, F. S. Bergeret, A. Kamra, L. Fu, P. A. Lee, and J. S. Moodera , Ubiquitous Superconducting Diode Effect in Superconductor Thin Films, Phys. Rev. Lett. 131, 027001 (2023).
  * [61] X. Y. Wei, O. A. Santos, C. H. S. Lusero, G. E. W. Bauer, J. B. Youssef, and B. J. van Wees, Giant magnon spin conductivity in ultrathin yttrium iron garnet films, Nat. Mater. 21, 1352 (2022).
  * [62] S. Knauer, K. Davídková, D. Schmoll, R. O. Serha, A. Voronov, Q. Wang, R. Verba, O. V. Dobrovolskiy, M. Lindner, T. Reimann, C. Dubs, M. Urbánek, and A. V. Chumak, Propagating spin-wave spectroscopy in a liquid-phase epitaxial nanometer-thick YIG film at millikelvin temperatures, J. Appl. Phys. 133, 143905 (2023).


A closed parameterization of DNA–damage by charged particles, as a function of
energy — A geometrical approach

Frank Van den Heuvel PhD∗1,2,

1 CRUK/MRC Oxford Institute for Radiation Oncology, Department of Oncology,
University of Oxford, Oxford, UK
2 Laboratory for experimental radiotherapy, Department of Oncology, University
of Leuven, Leuven, Belgium
$\ast$ E-mail<EMAIL_ADDRESS>
## Abstract

Purpose: To present a closed formalism calculating charged particle radiation
damage induced in DNA. The formalism is valid for all types of charged
particles and due to its closed nature is suited to provide fast conversion of
dose to DNA-damage.
Methods: The induction of double strand breaks in DNA–strings residing in
irradiated cells is quantified using a single particle model. This leads to a
proposal to use the cumulative Cauchy distribution to express the mix of high
and low LET type damage probability generated by a single particle. A
microscopic phenomenological Monte Carlo code is used to fit the parameters of
the model as a function of kinetic energy related to the damage to a DNA
molecule embedded in a cell. The model is applied for four particles:
electrons, protons, alpha–particles, and carbon ions. A geometric
interpretation of this observation using the impact ionization mean free path
as a quantifier, allows extension of the model to very low energies.
Results: The mathematical expression describes the model adequately using a
chi–square test ($\chi^{2}/NDF<1$). This applies to all particle types with an
almost perfect fit for protons, while the other particles seem to result in
some discrepancies at very low energies. The implementation calculating a
strict version of the RBE based on complex damage alone is corroborated by
experimental data from the measured RBE. The geometric interpretation
generates a unique dimensionless parameter $k$ for each type of charged
particle. In addition, it predicts a distribution of DNA damage which is
different from the current models.

## Introduction

The biological effect of ionizing radiation on human cells is believed to be
related to the generation of damage in the DNA–molecule located in the cell’s
nucleus[1]. The physical mechanism is the ionization of the DNA macro
molecule, generating lesions in the molecular structure, either by direct
ionization or by the generation of radicals in the vicinity of the DNA which
then indirectly damage it. These events (direct or indirect) can create
several types of damage to the DNA by combining a number of lesions into a
cluster, which can only happen if they occur in close proximity (typically
within one turn of the DNA–helix). The most prevalent of these damage types
are base damage (2 lesions), followed by single strand breaks (SSB) (3
lesions) , double strand breaks (DSB) (4 lesions), and locally multiple damage
sites (LMDS). The latter are clusters of different types of damage occurring
close to each other. It is shown that base as well as SSB damage is not likely
to be a deciding factor in the destruction of cells, due to the efficient
repair mechanisms which exist in the cell[2]. The combination of double strand
breaks and LMDS’s is likely to be the root cause for cell kill[3].

To quantify the amount of ionizing interactions in a medium, the physical
notion of dose can be used. Dose is defined as the amount of energy deposited
in a medium per unit mass and is expressed in Joule(J) per kg or Gray (Gy). In
the case of dose deposition by charged particles the Bethe–formalism is used.
This describes ionization events in a medium in terms of energy loss of the
charged particles in inelastic collisions with the electrons of the medium,
through the notion of mass stopping power ($dE/\rho dx$). In his seminal work
already in 1930, Bethe showed that there is an intimate relationship between
stopping power on the one hand, and energy (i.e. speed), charge, and the
medium in which the interaction takes place on the other hand[4]. A further
extension taking into account the possibility of the charged particle picking
up electrons, thereby changing the stopping power was introduced by Barkas[5],
using the concept of an effective charge.

In radiation biology, linear energy transfer (LET) is used rather than
stopping power. LET is identical to stopping power with the energy delivered
to $\delta$–rays (i.e. highly energetic knock on electrons) subtracted. This
quantity is called restricted stopping power. As such, LET is a measure for
the density of ionization taking place along the track of a charged particle
through a medium. Due to its close relationship with stopping power, it
follows that there is a close relationship between LET and the kinetic energy
of the depositing particle. From observation a dearth of DSB’s and LMDS’s was
shown to be related to low LET irradiations, while an increased number of both
for the same dose is seen high in hight LET irradiations[1]. Brenner and
Ward[6] argued that DSB and LMDS damage was related to multiple interactions
by single particles, rather than the combination of single strand breaks
generated by single particles. In the field of microdosimetry, this is taken a
step further by defining the notion of lineal energy which introduces the
amount of energy deposited along lines confined in a convex geometric shape
with a given distribution of cord lengths estimating the energy deposited in
various shapes, which can be used for measurement (i.e. spheres, cylinders).

Extending this, it is natural to propose a model where distance between
ionizations along these lines plays a significant role in the generation of
DNA–damage. A full listing and treatment of these quantities can be found in
the ICRU reports 16, 19, and 36 [7, 8, 9].

To describe the damage impact of charged particles on the DNA–structure, the
science community has taken its recourse to using Monte Carlo simulations to
quantify the damage introduced[10, 11]. A more fundamental analytical approach
is currently lacking, due to the underlying complexity of the DNA molecule,
and the paucity of the available experimental data. The data which is
available is mainly provided in terms of relative biological effective dose
(RBE), a quantity combining physical, spectral, chemical, and biological
factors, all of which hamper ab–initio calculations.

Monte Carlo calculations are able to predict the induction of simple or
complex damage as well as induction of single and double strand breaks in
DNA–molecules. These findings are interpreted using the Bethe–Barkas formalism
in terms of LET and show that high LET particles indeed introduce more complex
damage.

In this paper we develop a parameterization using a simple geometrical model,
that describes the behavior as calculated by the Monte Carlo codes. We also
show that this formalism describes the current knowledge well.

## Methods and Materials

### Theory

We use the single charged particle model as proposed by Brenner and Ward,
distinguishing three types of interaction results: Low LET, high LET, and
intermediate LET mode. The specifics of each mode are explained below.

  1. 1.

Low LET: A single particle is generally not able to generate lesions close
enough together to induce double strand breaks at each interaction. It is
clear that DSB’s can be generated but in a limited fashion and that we use the
word lesion in liberal fashion to indicate an interactive event which has
damage as a consequence.

  2. 2.

High LET: The particle has the possibility to generate multiple lesions
irrespective of any geometrical considerations. We implicitly assume that the
double strand break damage is the result of multiple interactions by one
particle. How exactly this damage is introduced (direct or indirect) is
outside the scope of this article. An implicit assumption however is that
ionizing events need to be geometrically close to the DNA structure.

  3. 3.

Intermediate: In given geometric circumstances it is possible for the charged
particle to generate DSB–damage, in a high–LET manner, depending on the angle
under which the particle hits the sensitive volume (Fig. 1).

As a surrogate to categorize the charged particle in one of the types defined
above, we use the mean path length between ionizing interactions in a medium
consistent with the atomic make up of a DNA–molecule for the type of particle
under consideration. In the remainder, we denote this with $\lambda(E)$, where
$E$ is the kinetic energy of the particle. If $\lambda(E)$ is large relative
to the sensitive volume, then the lesions on average are too far apart and
only damage types related to a few lesions can occur (i.e. SSB and base
damage). Charged particles with such energies will be part of the first
category. If on the other hand $\lambda(E)$ is small then the probability of
lesions creating more complex clusters of damage close together will be
higher. Charged particles with this property will be in the high LET category.
Finally, charged particles with intermediate distances between ionizing events
have the capability of generating DSB and LMDS damage depending on other
factors than $\lambda(E)$ alone. In this model we use the geometric direction
of the path of the charged particle as a parameter. In Figure 1 a schematic
model of this approach is shown. This implies that only a limited amount of
directions are available to contribute to the amount of complex damage in the
manner as outlined for the high–LET type interactions. This happens when for a
particle of a given energy the quantity $\lambda(E)$ is slightly larger than
the maximal distance between two DNA–damage lesions to be considered as being
in the same cluster (usually about 10 base pairs (bp)). Due to the finite
thickness of the sensitive volume it is possible to behave in a high–LET
fashion depending on the angle with which the particle’s path crosses the
volume. This occurs when the projection of the path is smaller than the
previously determined maximum.

Figure 1:  A schematic model of a source of charged particles with a given
mean free path length (i.e. a given energy), which is comparable with the
diameter of the sensitive cell volume. As the angle ($\theta$) of the
particle’s path with respect to the normal to the axis of the structure
increases the chance that more than a single event will occur in the volume.
This implies that if the angle ($\theta$) is larger than the one for which the
projection of the average length between interactions equals the diameter,
more high LET events will be registered.

### Equivalence Principle

In the case of irradiation with charged particles all directions of the
particle’s paths are possible as are all rotational positions of the
DNA–structure. A particle that interacts (i.e. that creates a lesion) at the
surface of a given sensitive volume has limited possibilities to interact
again given that on average, a specific distance (which depends on the
particle energy) has to be travelled before it interacts again. The next
interaction’s position is then limited by the constraints outlined above if it
is to fall within the sensitive volume. This first interaction can happen
anywhere along the volume, but the constraints are relative to the position of
that point. This implies that we can invoke an equivalence principle and
reduce the problem to that of an isotropic point source positioned at the
surface of a sensitive volume.

### Mathematical expression of the equivalence principle

We need to calculate what fraction of the paths starting in the given point
can interact with the sensitive volume given the fact that there is a length
within which this is not likely, provided by $\lambda(E)$, and that there is a
maximal distance ($H$) that disqualifies the generated lesion to be registered
in the same cluster. We have reduced this problem to that of the distribution
of projections of a point source on a line–piece, the solution to which is
known as the Cauchy–distribution[12], and is described by the Lorenz function
$f(x)$ with $x\in\Re$ expressed as follows:

$f(x)~{}=~{}\frac{1}{1+(\frac{x}{r})^{2}}$ (1)

Figure 2:  The abstracted version of Figure 1 describes the distribution of
horizontal distances at which a line segment tilted at a random angle $\theta$
cuts the x–axis. Only particles with large angles contribute to double strand
break events by combining damage generated by a single particle. The red line
indicates the “forbidden” area as (on average) this distance between the
ionization events is observed.

From Figure 2, it follows that the contribution $P$ for a given energy of
charged particles to high LET events is proportional to:

$P~{}\sim~{}2\int_{H-\lambda(E)}^{H}\frac{1}{1+(\frac{x}{r})^{2}}dx$ (2)

Performing the calculation, we obtain:

$\displaystyle P~{}$
$\displaystyle\sim~{}\frac{2}{\pi}[\tan^{-1}(\frac{x}{r})\biggm{|}_{H-\lambda(E)}^{H}]+C$
(3)
$\displaystyle\sim~{}\frac{2}{\pi}[\tan^{-1}(\frac{H}{r})-\tan^{-1}(\frac{H-\lambda(E)}{r})]+C$
(4)
$\displaystyle\sim~{}\frac{2}{\pi}[\tan^{-1}(\frac{\lambda(E)-H}{r})+\tan^{-1}(\frac{H}{r})]+C$
(5)

This implies that the amount of DSB–damage for a given dose and given energy
of the charged particle is governed by the following expression.

$F_{cd}(E)~{}=~{}(a-b)\frac{2}{\pi}[\tan^{-1}(\frac{\lambda(E)-H}{r})]+b$ (6)

The change from a low to a high LET regimen occurs over a small energy
interval. In such a small interval the average distance dependence on the
energy of the particle can be approximated with a linear function. Therefore,
we expect the energy dependence of the contribution of complex damage to
follow the same form as in Equation 6 yielding the following expression, with
H=$\lambda(E_{0})$ and r=$\lambda(\Gamma/2)$, $E_{0}$ being the energy, where
the change in DSB is maximized and $\Gamma$ a measure for the width of the
slope (i.e. the full width at half maximum in differential energy space).

$F_{cd}(E)~{}=~{}(a-b)\frac{2}{\pi}[\tan^{-1}(\frac{E-E_{0}}{\Gamma/2})]+b$
(7)

With the parameters $a$ and $b$ related to the levels 1 and 2 as outlined
above. From boundary conditions we find that at very large energies (i.e.
$E>>E_{0}$) the expression is reduced to minimal number of double interactions
($D_{min}$) which is equal to $a$. The value of $b$ is related to the maximal
number of double interactions ($D_{max}$) as follows:

$D_{max}=D_{min}+(1-\tan^{-1}(-\frac{E_{0}}{\Gamma/2}))b.$ (8)

The formalism using energy alone allows us to forego specific assumptions
regarding the dimensions of the DNA–molecule. Furthermore, it also allows us
to apply this technique to particles where the values for $\lambda(E)$ are
less well known. In addition, it allows us to test this formalism using
experimentally available data which is available as a function of energy.

###  Validation using Monte Carlo Simulations

The use of microdosimetric calculations has provided important insight into
the mechanisms and effects of radiation deposition. In the past, Monte Carlo
simulations of charged particle deposition by various modalities were used to
quantify and typify the kinds of damage introduced by the different
modalities[10].

The Monte Carlo Damage Simulation code (MCDS) developed by Semenenko and
Stewart, generates spatial maps of the damaged nucleotides forming many types
of clustered DNA lesion, including single-strand breaks (SSB), double strand
breaks (DSB), and individual or clustered base damages[11]. This approach has
been shown to yield a linear relationship of the number of generated DSB’s up
to a high dosage. It follows that this parameterization also provides the
possibility to link dose to damage. In this paper, MCDS version 3.0 was used
with the parameters described below. The DNA length which was chosen to be
1Gbp (Giga base pairs) and a nucleus diameter of 5$\mathrm{\mu m}$. In the
MCDS software, the geometry of the DNA–molecule is not an explicit parameter.
Here four parameters are used: 1) the DNA–segment length $n_{seg}$, which is
an ad hoc parameter expressed as base pairs $Gy^{-1}cell^{-1}$, 2) the number
of strand breaks generated $\sigma_{Sb}$, 3) the number of base pair damages
generated $\sigma_{Bb}$ by defining $f=\sigma_{Bb}/\sigma_{Sb}$, and 4) a
parameter $N_{min}$ (bp) describing the minimal separation for damage to be
apart not to be counted as being in the same cluster. The values of these
parameters is determined on the basis of other simulations and measurements.
For a more in–depth treatment of these parameters we refer to the work by
Semenenko and Stewart[13]. Variable input parameters MCDS were; the modality
(i.e. energy depositing particle (electron, proton,…)), the energy (in MeV),
and the oxygen concentration in %. In the implementation described here we
chose to omit any oxygen enhancement as this could be a confounding factor and
is the subject of another study. In this study it was found that oxygen only
changed the amount of damage in the low LET regimen, leaving the formalism
unchanged (data not shown). Therefore, a concentration of 0% oxygen was used.
For every particle type at the relevant kinetic energies, all complex damage
was noted per Gy, per cell and per kinetic energy.

### Fitting procedure

The ultimate goal was to fit the complex damage function to the data obtained
by the Monte Carlo simulation. The parameters that need fitting are the energy
position ($E_{0}$) the width of the underlying Cauchy distribution ($\Gamma$)
and the parameters $a$ and $b$. If a regular fit (i.e. all parameters fit at
the same time) is performed we see strong co–variances between the parameters.
To come to meaningful results we opted to perform a two step procedure: First,
we eliminate the parameters $a$ and $b$ by fitting the differential, thereby
reducing expression 6 to the Lorenz function.

$\frac{dF_{cd}}{dE}=\frac{\Gamma^{2}/4}{(E-E_{0})^{2}+\Gamma^{2}/4}$ (9)

This is also mathematically equivalent to the fit of a Breit–Wigner resonance
in high energy physics[14]. In a second fit–procedure, the remaining variables
$a$ and $b$ are fit using the cumulative Cauchy function. The fitting
procedures were performed in the gnuplot111http://www.gnuplot.info–software
using a Levenberg–Marquardt minimization routine.

## Results

In Figure 3, the Lorenz expression as outlined in Equation 6 together with a
normalization factor, is used to fit the energy differential probability for
the generation of DSBs. The fit is performed to minimize the $\chi^{2}$–value.
In all cases, the resulting $\chi^{2}/NDF$ (NDF = Number of Degrees of
Freedom) are lower than 1. The values of the parameters are provided in Table
1.

(a) Electrons (b) Protons (c) $\alpha$–particles (d) Carbon Ions

Figure 3:  Fitting the Cauchy expression to the energy differential probability of generating DSB’s denoted $\frac{d\sigma}{dE}$ Particle | $\Gamma$ | $E_{0}$ | $a$ | $b$
---|---|---|---|---
$e^{-}$ | (2.854$\pm$0.051)$10^{-04}$ MeV | (1.05736$\pm$0.036)$10^{-04}$ MeV | 2.9061 | 21.460
$p^{+}$ | 0.5575$\pm$0.0094 MeV | 0.1642 $\pm$ 0.0037 MeV | 2.89068 | 21.4273
$\alpha^{++}$ | 8.20$\pm$0.17 MeV | 3.1850$\pm$0.056 MeV | 3.0856 | 20.7933
C6+ | 201.7$\pm$8.4 MeV | 95.4$\pm$2.5 MeV | 3.01459 | 21.8489
Table 1:  The different values for $\Gamma$ and $E_{0}$ as defined by Equation
9 and obtained from a fitting procedure together with the asymptotic standard
error of the fitted parameter. All fits exhibited minimal values of
$\chi^{2}/NDF$ (NDF = Number of Degrees of Freedom). The columns $a$ and $b$
are the parameters indicating the levels of DSB at low, resp. high LET. Note
that even in low LET the number of DSB’s is not zero as complex damage can
occur due to the combination of simple damage events.

All fits are completely satisfactory at energies higher than $E_{0}$. On the
lower energy side some discrepancies can be observed depending on the incoming
particles, particularly in the case of electrons and carbon ions. We refer the
reader to the discussion section. For protons we see a satisfactory fit over
the full energy range.

Figure 4 shows the final results with all parameters fit. Again, all fits have
$\chi^{2}$–values commensurate with a positive goodness of fit. The final
values and the standard errors for the fitted parameters are listed in Table
1. Note, that the noise in the differential curves increases as the particles
become heavier. The random-seeming errors in the estimates of the derivative
arise in part from the Monte Carlo estimates of the mean number of DSB per Gy
per Gbp and from numerical instabilities associated with the calculation of
the derivative using finite difference methods.

(a) Electrons (b) Protons (c) $\alpha$–particles (d) Carbon Ions

Figure 4:  The prediction of the number of double strand breaks or more
complex damage as a function of energy for 4 relevant charged particles. This
provides the number of Double Strand breaks (DSB) per Gy, Gbp and per cell.
The prediction for protons and alpha particles is almost perfect. For
electrons and carbon ions some discrepancies exist at lower energies.

### Geometric approach

Now is the time to investigate the geometric interpretation further. To
quantify the function $\lambda(E)$ we can use the inelastic mean free path as
measure (IMFP). Values for IMFP for electrons are well known in the
literature, not in the least as they are important in solid state physics and
electron microscopy. They can be found in freely available databases for a
variety of elements and compounds, even for organic molecules like DNA[15].
Proton values can be found in a publication by Zhen–Yu and colleagues[16]. For
heavier particles such as $\alpha$–particles and carbon–ions, the data is more
difficult to find. We therefore opt not to use the data for these particles
and restrict ourselves to electrons and protons in this further treatment.

In all current microdosimetric codes, the Bethe formalism is used which is
valid for higher energies (i.e. above 500eV for electrons). This implies that
changes in IMFP, denoted by $\lambda$, which impact the damage calculated
using these codes, also reflect the limitations of the Bethe formalism. From
the theory the following expression is used:

$\lambda(E)=\frac{E}{A\log(E/E_{0})+B}$ (10)

Particle | A | B
---|---|---
Electrons | 69.200 eV/nm | -153.94 eV/nm
Protons | 115.231 keV/nm | -301.45 keV/nm
Table 2:  Parameters obtained by fitting Eq. 10 to data obtained from NIST
(electrons) and Zhen Yu et al. (protons)

In this work the parameters $H$ and $r$ have thus far not been linked to any
physical property but were fit. An interesting proposition could be to link
these to dimensions of the target structure. Indeed, the choice of a cylinder
as a geometric representation is not an accident. It is natural to use the
diameter of a DNA–molecule as a measure of the cylinder’s diameter. The length
of the cylinder is then related to the maximal distance we allow to classify
two damage events, being part of the same cluster of complex damage. Both
values can readily be found in the literature and text books[17]. For the most
prevalent form of cellular DNA (B–DNA), the values are $3.4$nm (i.e. the
height of a spiral of 10 base pairs), and $2.37$nm as the diameter. We now
define a dimensionless quantity $k$ which is specific to the type of charged
particle used. It is clear that this parameter acts as a scaling parameter but
also depends on the ratio of both fixed parameters. Equation 6 now reads as
follows:

$F_{d}(E)~{}=~{}(a-b)\frac{2}{\pi}[\tan^{-1}(\frac{k\lambda(E)-3.4}{2.37})]+b$
(11)

This reduces the impact of the charged particle’s energy on the induction of
complex damage in a DNA–molecule to three parameters $a$, $b$, and $k$. Figure
5 illustrates the use of these parameters and shows that comparable results to
the energy–based formalism can be obtained. It follows that we can repeat the
fitting procedure keeping $a$ and $b$ from the expression based on energy (Eq
7). We find values of k=5.18 for electrons and k=4.82 for protons.

### Extending the model

In the work presented above as well as in the used Monte Carlo simulations,
the Bethe–Barkas formalism together with its flawed approach in the lower
energy regions has always been used. It is well established that the IMFP does
not follow the expression outlined in Equation 10, where $\lambda(E)$ keeps
diminishing as the energy diminishes. Indeed, when the energy is lower than
200eV an increase in IMFP is observed due to plasmonic effects[18]. Ziaja et
al[19] showed that it is possible to describe this behavior analytically by
extending Equation 10 with a second term as follows:

$\lambda(E)_{Z}~{}=~{}\frac{\sqrt{E}}{A_{1}(E-E_{th})^{B_{1}}}~{}+~{}\frac{E-E_{0}\exp(-B/A)}{A\log(E/E_{0})+B}$
(12)

In this equation the parameter $E_{th}$ serves as a threshold separating the
behavior as described by Bethe from the plasmonic interactions. Using the data
provided in the work from Zhen–Yu and colleagues[16] it is straightforward to
obtain parameters for the behavior of protons. These are presented in Table 3.

Particle | $A_{1}$ | $B_{1}$ | $E_{th}$ | $A$ | $B$ | $E_{0}$
---|---|---|---|---|---|---
Electrons | 0.6560 | 1.0100 | 24.2838 | 65.898 | -128.23 | 1.0
Protons | 0.681 | 1.249 | 42.38 | 117.01 | -318.7 | $1.0\times 10^{3}$
Table 3:  Parameters as in Table 2 with added lower energy factors. The
fitting was performed using Eq. 12

To extend our model to incorporate the behavior of very low energy particles
it is sufficient to replace the expression $\lambda(E)$ by $\lambda(E)_{Z}$ in
equation 11. In Figure 5, the modified curves show the difference with the
calculations based on the Bethe formalism only. This also shows that there is
an upper limit to the increase in DSB’s which depends on the type of particle.
It is conceivable that this approach also works for the heavier particles
which can be seen when using the IFMP’s in water for these (not shown).

Figure 5:  Using the quantities for $H$ and $r$, the dimensionless constant
$k$ for electrons (left) and protons (right) is determined. Using both the
limited expression for $\lambda(E)$ and the more accurate estimate
$\lambda_{Z}(E)$. The former provides a fit to the Monte Carlo data comparable
with the results obtained using the energy–based formalism. The second
approach provides a maximal complex damage yield which differs for electrons
and protons.

## Discussion

We developed an approach to predict damage in complicated situations where
fields of different charged particles and their respective energy spectra
impact on living cells. The approach, due to its analytical nature allows very
fast calculation of damage in otherwise long simulations. In the derivation of
this approach using energy alone there are no assumptions on the mechanics
with which DNA–damage is caused by the charged particles. The only assumption
is that there is a sensitive volume where, if ionizations take place, damage
is introduced in the DNA. How exactly this damage is caused is not specified.
In the remainder of the text a parameter is identified, the average distance
between ionizations for the given charged particle in the medium ($\lambda$).
We show that this approach adequately quantifies the results from Monte Carlo
simulations based on phenomenological data and reduces these to a closed
analytical expression whereby the type of charged particle is expressed by a
single parameter ($k$). On the other hand we should be aware that issues like
repair mechanisms and oxygen effects are not present in the model, making its
applicability limited. However, if all things are identical (i.e. the type of
cells, oxygenation, etc…) and the only thing different is the type of charged
particle and its energy, then the original damage introduced in the DNA
structure should correlate with the outcome. An underlying assumption here is
that the repair processes are somehow independent from the modality with which
the cell is irradiated.

The results of this approach can be applied to determine the biological impact
of radiation in mixed environments, as in the case of proton therapy, where
protons, electrons and heavier ions (due to neutrons), deposit energy. Other
approaches have been proposed to try to predict outcomes from mixed fields,
which are based on available clinical response data. Most notably, an approach
based on the local effect model (LEM), where macroscopic response data in the
form of dose–effect curves is used to quantify the relative effect of the dose
delivered[20]. The parameterization, however, of the latter approach is
extensive due to the fact that every effect curve has two parameters for a
given $\alpha/\beta$–value, making the model over–parameterized. As such, it
is possible to have this model reflect the current knowledge of dose and
modality response adequately, which forms an important, albeit controversial
tool[21, 22]. Its power to predict the behavior outside of the current
knowledge therefore seems to be limited.

Cucinotta et al. attempted to incorporate the volumetric properties of the
dose deposition[23] to account for differences in track structure. They
observed that: “LET is a poor descriptor of energy deposition in small volumes
because of the diffusion of secondary electrons out of the volume and
contribution of $\delta$–rays that pass outside of the volume”. To address
this problem a quantification of the energy distribution of generated
secondary particles, or $\delta$–rays was proposed.

Such a secondary charged particle indirectly changes the behavior with respect
to the DNA damage induced. Indeed, depending on the median energy of the
spectrum the DNA damage changes accordingly if the dose is kept constant. In
the paper presented here this behavior could be easily incorporated by
considering the DNA damage for all the particles (i.e. ions and $\delta$–rays)
separately using a methodology modeled on the use of the electronic
equilibrium concept in photon cavity theory. Currently this behavior is hidden
in the $k$ parameter and it would be interesting to see if such an approach
will lead to a convergence of all $k$–values for all particles.

To take these actions fully into account an approach to provide a more
detailed model of the biological effect directly in the Monte Carlo simulation
is proposed by Sato et al.[24] This would, in theory, allow a direct
calculation of the effect in terms of energy deposited. However, as outlined
by Cucinotta this is not without problems as the behavior of low energy
electrons needs to be adequately modelled. This work predicts that the current
knowledge using the Bethe formalism, might not be suitably extended.

The results from the geometric interpretation indicate that the overall
behavior of the DNA damage induction is identical for all types of charged
particles. The only difference is in the dimensionless parameter $k$. The
latter seems to change as the ion used is heavier. Preliminary calculations
using the IMFP in water indicate that the value of $k$ diminishes as the
charged particles used are heavier (or more charged, data not shown). A
possible reason for this is that the track structure can be quite different
for different charged particles. This fact could also be an explanation for
the discrepancy found at very low energies for carbon–ions. Indeed, allowing
the parameter $k$ to be covariant with the other parameters, does provide a
more adequate fit (data nor shown).

The results for the electrons also shows a discrepancy with regard to the
generation of complex damage at lower energies. For electrons, the data on
very low energy electrons are not available in terms of energy deposition.
Indeed, the model proposed here shows a much lower incidence of complex damage
due to plasmonic effects in that region.

In summary, the model proposed here allows extension to very low energies for
electrons and protons.The fact that there are indications that the induction
of DSB’s varies linearly with dose, provides an easy implementation to dose
planning systems, given the knowledge of dose deposition spectra in a
treatment beam. An example of such implementation is provided in the Appendix.

## Appendix: Implementation in dose deposition calculations

### Mono–energetic treatment

In dose calculations a dose matrix is obtained on a dose grid Let
$\mathbf{D}=D[i,j,k]$ be the dose matrix provided. Then we can write the
amount of complex damage incurred by particles with an energy ($E$) as a
damage matrix, denoted as ($M_{cd}$). as follows:

$\displaystyle\mathbf{M_{cd}}=M_{cd}[i,j,k]$ $\displaystyle=\mathbf{D}\times
F_{cd}(E)$ (13)

$F_{cd}(E)$ then denotes a response function converting dose to damage.

### Poly–energetic treatment

Dose deposition spectra rarely consist of a field of mono–energetic electrons.
For a photon source with a given photon spectrum, an energy depositing
electron fields exists, which is roughly constant throughout the target
volume. Using Monte Carlo simulations it is possible to calculate this field
and its spectrum $\Psi(E)$. It then becomes possible to include the spectrum
in the calculation of the damage matrices. This approach has been used already
by different authors [25, 26].

$\mathbf{M_{cd}}~{}=~{}\mathbf{D}\times\frac{\int_{0}^{E_{max}}\Psi(E)F_{cd}(E)dE}{\int_{0}^{E_{max}}\Psi(E)dE}$
(14)

In the case of charged particle treatment, the particles are moderated and the
energy spectrum changes depending on the position of the point where the dose
is being deposited. It is therefore necessary to apply Equation 14 to each
point separately with knowledge of the depositing energy spectrum in that
point. Due to the closed nature of the formalism developed in this paper, it
becomes feasible to use off the shelf computing equipment.

### Application: Proton treatment

Recently, the coupling of Monte Carlo simulations in dose deposition to micro-
dosimetric code has been proposed and applied by several groups[26, 25]. Here
a two step approach is followed; 1) a general purpose Monte Carlo code (MCNPX
2.7b)[27] is used to estimate the spectrum of all different dose contributing
particles, 2) a micro dosimetric code[13] is used to determine the biological
damage.

The framework for conversion of dose to biological effect is implemented on a
simulation of a pristine 200MeV proton beam, taking into account the changing
proton spectrum. The proton simulation is performed using MCNPX. Figure 6(a)
shows the variation of the number of complex damage events as a function of
energy of the proton. In addition, the spectrum of depositing protons is shown
at a position before the Bragg peak and at the Bragg peak. In Figure 6(b) the
effect on the dose deposition is shown together with the $RBE_{cd}$ calculated
as the complex damage yield generated by the protons at that particular
position, divided by the complex damage induced by a 6MV photon beam with the
same spatial characteristics. Note, that the $RBE_{cd}$ is of the order of 1.1
with larger value of 2 a few mm distal from the Bragg peak. This is
commensurate with the cell data reported by Paganetti et al.[28] and Chaudhary
et al.[29], who showed that the radiobiological effect at the distal end of a
spread out bragg peak increases, a fact predicted by Goitein[30]. Currently,
data of direct measurement of DNA–damage in–vitro along a proton beam are
scarce. The advent of $\gamma$–H2AX measurements, as a marker for DSB–damage
is promising in this regard and has been used to investigate anti–protons[31].

(a) Impact of the spectrum (b) Proton dose deposition

Figure 6: In Fig 6(a) the spectrum at the beginning of the Bragg peak (scaled
by 0.1) is completely within the low LET regimen. While at the end of the
Bragg peak a significant part of the dose depositing protons exhibits high LET
characteristics. Fig.6(b) shows the $\mathrm{RBE_{cd}}$ (red line) together
with the damage induced by a mono–energetic proton beam.

## Acknowledgements

I am indebted to Rob Stewart, not only for generously providing the MCDS–code
for everyone to use, but also for providing much needed input on the
biological aspects of the MCDS implementation. Both Ricardo Raabe and Mike
Partridge reviewed some of the physics and made it a more rigorous paper.
Sandra Nuyts and Dirk De Ruysscher’s clinical input was also highly
appreciated.

## References

  * 1. Hall EJ (1978) LET and RBE.  In: Radiobiology for the radiologist, Philadelphia: Harper & Row, Publishers.
  * 2. Caldecott KW (2008) Single-strand break repair and genetic disease.  Nat Rev Genet 9: 619–631.
  * 3. Ward JF (1985) Biochemistry of dna lesions.  Radiation Research Supplement 8: pp. S103-S111.
  * 4. Bethe H (1930) Zur theorie des durchgangs schneller korpuskularstrahlen durch materie.  Ann Phys 5.
  * 5. Barkas W, Dyer J, Heckman H (1963) Resolution of the $\sigma$–mass anomaly.  Phys Rev Lett 11: 26.
  * 6. Brenner DJ, Ward J (1992) Constraints on energy deposition and target size of multiply damaged sites associated with dna double-strand breaks.  International Journal of Radiation Biology 61: 737-748.
  * 7. ICRU, International Commission on Radiation Units and Measurements (1970) Linear energy transfer.  ICRU Report 16.
  * 8. ICRU, International Commission on Radiation Units and Measurements (1971) Radiation quantities and units.  ICRU Report 19.
  * 9. ICRU, International Commission on Radiation Units and Measurements (1983) Microdosimetry.  ICRU Report 36.
  * 10. Nikjoo H, O’Neill P, Goodhead DT, Terrissol M (1997) Computational modelling of low-energy electron-induced DNA damage by early physical and chemical events.  Int J Radiat Biol 71: 467-473.
  * 11. Stewart RD, Yu VK, Georgakilas AG, Koumenis C, Park JH, et al. (2011) Effects of radiation quality and oxygen on clustered DNA lesions and cell death.  Radiation Research 176: 587-602.
  * 12. Cauchy A (1853) Sur les résultats moyens d’observations de même nature, et sur les résultats les plus probables.  Comptes Rendus de l’Académie des Sciences : 198-206.
  * 13. Semenenko V, Stewart R (2006) Fast Monte Carlo simulation of DNA damage formed by electrons and light ions.  Phys Med Biol 51: 1693-1706.
  * 14. Breit G (1959) Handbuch der Physik XLI/1.  Berlin, Heidelberg: Springer.
  * 15. Powell C, Jablonski A (2010) Nist standard reference database 71.  In: NIST Electron Inelastic–Mean–Free–Path Database — Version1.2, NIST Gaithersburg, MD.
  * 16. Zhen-Yu T, Yue-Yuan X, Ming-Wen Z, Xiang-Dong L (2010) Proton inelastic mean free path in a group of organic materials in 0.05-10mev range.  Chinese Physics Letters 27: 113403.
  * 17. Sinden R (1994) DNA structure and function.  San Diego: Academic Press.
  * 18. Tanuma S, Powell CJ, Penn DR (2005) Calculations of electron inelastic mean free paths.  Surface and Interface Analysis 37: 1–14.
  * 19. Ziaja B, London RA, Hajdu J (2006) Ionization by impact electrons in solids: Electron mean free path fitted over a wide energy range.  Journal of Applied Physics 99: 033514.
  * 20. Kraemer M, Scholz M (2006) Rapid calculation of biological effects in ion radiotherapy.  Phys Med Biol 51: 1959-1970.
  * 21. Katz R (2003) The parameter-free track structure model of Scholz and Kraft for heavy-ion cross sections.  Radiation Research 160: 724–728.
  * 22. Paganetti H, Goitein M (2001) Biophysical modelling of proton radiation effects based on amorphous track models.  International Journal of Radiation Biology 77: 911-928.
  * 23. Cucinotta FA, Nikjoo H, Goodhead DT (2000) Model for radial dependence of frequency distributions for energy imparted in nanometer volumes from hze particles.  Radiation Research 153: 459–468.
  * 24. Sato T, Kase Y, Watanabe R, Niita K, Sihver L (2009) Biological dose estimation for charged-particle therapy using an improved phits code coupled with a microdosimetric kinetic model.  Radiation Research 171: 107–117.
  * 25. Van den Heuvel F, Locquet JP, Nuyts S (2010) Beam energy considerations for gold nano-particle enhanced radiation treatment.  Physics in Medicine and Biology 55: 4509.
  * 26. Hsiao Y, Stewart RD (2008) Monte carlo simulation of DNA damage induction by X-rays and selected radioisotopes.  Phys Med Biol 53: 233.
  * 27. Waters LS, McKinney GW, Durkee JW, Fensin ML, Hendricks JS, et al. (2007) The MCNPX Monte Carlo radiation transport code.  In: Albrow, M and Raja, R, editor, Hadronic Shower Simulation Workshop. Amer Inst Physics, volume 896 of _AIP conference proceedings_ , pp. 81-90.  Hadronic Shower Simulation Workshop, Batavia, IL, SEP 06-08, 2006.
  * 28. Paganetti H, Niemierko A, Ancukiewicz M, Gerweck LE, Goitein M, et al. (2002) Relative biological effectiveness (RBE) values for proton beam therapy.  International journal of radiation oncology, biology, physics 53: 407–421.
  * 29. Chaudhary P, Marshall TI, Perozziello FM, Manti L, Currell FJ, et al. (2014) Relative biological effectiveness variation along monoenergetic and modulated bragg peaks of a 62-mev therapeutic proton beam: A preclinical assessment.  Int J Radiat Oncol Biol Phys 90: 27–35.
  * 30. Goitein M (2008) Radiation Oncology: A Physicist’s-Eye View.  Springer Verlag.
  * 31. Kavanagh JN, Currell FJ, Timson DJ, Savage KI, Richard DJ, et al. (2013) Antiproton induced dna damage: proton like in flight, carbon-ion like near rest.  Sci Rep 3: –.