# Cancelation of dispersion and temporal modulation with non-entangled
frequency-correlated photons

Víctor Torres-Company  Electrical and Computer Engineering Department, Purdue
University, West Lafayette, 47906-IN, USA  Alejandra Valencia  ICFO-Institut
de Ciencies Fotoniques, Mediterranean Technology Park, 08860 Castelldefels
(Barcelona), Spain  Martin Hendrych  ICFO-Institut de Ciencies Fotoniques,
Mediterranean Technology Park, 08860 Castelldefels (Barcelona), Spain  Juan P.
Torres<EMAIL_ADDRESS>ICFO-Institut de Ciencies Fotoniques, Mediterranean
Technology Park, 08860 Castelldefels (Barcelona), Spain  Department of Signal
Theory and Communications, Universitat Politecnica Catalunya, Campus Nord D3,
08034 Barcelona, Spain

###### Abstract

The observation of the so-called dispersion cancelation and temporal phase
modulation of paired photons is generally attributed to the presence of
frequency entanglement between two frequency anticorrelated photons. In this
paper, it is shown that by introducing the appropriate amount of chromatic
dispersion or phase modulation between non-entangled photons, it is also
possible to observe these effects. Indeed, it is found that the relevant
characteristic for the observation of dispersion cancelation or the
cancelation of temporal phase modulation is the presence of certain frequency
correlations between the photons.

###### pacs:

42.25.Kb, 42.50.Dv, 42.50.Ar

##  I Introduction

The theory of quantum coherence glauber1963 ; glauber2006  describes the
temporal and frequency characteristics of a stream of photons through a
hierarchy of correlation functions. In particular, the second-order
correlation function determines the probability to detect a photon at a
certain location and instant time in coincidence with a companion photon at
another location and time. This magnitude is often used to test the quantum or
classical nature of a light source mandel_book . For instance, the normalized
correlations of classical fields (i.e., fields with a positive
$P$-representation) obey certain inequalities, whose violation constitutes an
unequivocal signature of the non-classical properties of the light kimble .

A particular type of correlation between photons is entanglement. Paired
photons that show frequency entanglement can be generated by means of the
process of spontaneous parametric downconversion (SPDC). Interestingly, in an
SPDC process, classical-like features of the fluctuations of the beams coexist
with the strong nonclassical photon-pair correlations loudon .

Two important effects have been attributed to the existence of frequency
entanglement and its demonstration have made use of the correlations existing
between signal and idler photons generated in SPDC. The first effect is
dispersion cancelation franson1992 ; gisin1998 ; kim2009 , which is observed
in the temporal domain. Briefly, if the signal and idler photons generated in
SPDC pumped by a continuous-wave (CW) beam are sent through two separate
dispersive optical elements, such as single-mode fibers, with corresponding
group-delay-dispersion (GDD) coefficients $\Phi_{1}$ and $\Phi_{2}$, the
temporal width of the second-order correlation function increases as
$(\Phi_{1}+\Phi_{2})^{2}$, in a similar way to the broadening of a pulse that
propagates in an optical fiber due to chromatic dispersion valencia2002 . So,
if the GDD parameters of both dispersive media are identical but of opposite
sign, the broadening of the second-order correlation function can be
suppressed. This is in contrast to the case where two identical broadband
coherent light pulses propagate through two dispersive media such as single-
mode optical fibers. In this case, the cross-correlation width broadens as
$(\Phi_{1}^{2}+\Phi_{2}^{2})$ and therefore the cancelation of the dispersion
effects is never possible franson1992 .

The second effect is remote temporal modulation of entangled photons
harris2008 ; harris2009 , which is observed in the frequency domain. For a CW-
pumped SPDC process, the detection of a signal photon at frequency
$\omega_{1}$ would only coincide with the detection of an idler photon at
frequency $\omega_{2}$, given that $\omega_{1}+\omega_{2}=2\omega_{0}$, where
$2\omega_{0}$ is the frequency of the CW pump beam. When synchronously driven
temporal modulators are placed in the signal and idler paths, respectively,
new frequency correlations appear. In a similar manner to dispersion
cancelation, if the two identical modulators are driven in opposite phase,
their global effect is to negate each other and the spectral correlations
appear as those when there are no phase modulators present.

It becomes of fundamental relevance to determine whether these effects are due
to the classical or the quantum-like behavior of the SPDC source. Several
authors have shown that similar dispersion cancelation effects can be obtained
with non-entangled light victor2009 ; shapiro2010 . In particular, the work in
victor2009  considers classical thermal light equally split in two dispersive
arms and demonstrates that the broadening of the second-order correlation
function increases as $(\Phi_{1}-\Phi_{2})^{2}$, so that it remains unaffected
if $\Phi_{1}=\Phi_{2}$. Later, the work in shapiro2010  shows that by
introducing a phase conjugator element in one arm of the intensity
interferometer (producing a Gaussian-state light source), the broadening of
the correlation function increases as $(\Phi_{1}+\Phi_{2})^{2}$, similar to
the signal-idler photon pairs from SPDC, and thus the same dispersion
compensation rules apply.

In this work, we provide new insights into the effects of dispersion
cancelation and temporal modulation applying the Heisenberg picture. This
formalism allows us to include the case of multiphoton pair generation in a
straightforward way, and thus to calculate all relevant coherence functions
(signal-signal and signal-idler correlations). Two important findings are
revealed: First, there is a background term in the second-order correlation
function for all types of correlations, but the signal-to-background ratio for
signal-signal correlations is lower than for the signal-idler correlations
franson2009 ; franson2010 . Second and more importantly, dispersion
cancelation and temporal modulation can also be observed with the correlation
existing among photons from an individual beam of SPDC (i.e., either signal-
signal or idler-idler), which do not show entanglement. This last point
illustrates that it is the existence of certain frequency correlations between
photons, rather than the entanglement, that is the key enabling factor that
allows the observation of dispersion cancelation and remote temporal
modulation.

##  II Quantum description of the light generated in Spontaneous Parametric
Downconversion

The effects studied in this work, i.e., remote dispersion cancelation and
cancelation of temporal modulation are observed in the corresponding second-
order correlation functions. To calculate these magnitudes one needs to
establish the dependence of the creation and annihilation operators with the
physical parameters set by the SPDC process. In this section, we review the
quantum theory of the SPDC process pumped by a CW pump and establish this
dependence.

Let us consider degenerate SPDC produced by a CW plane-wave pump with a
central frequency $2\omega_{0}$ in a crystal of length $L$ and nonlinear
coefficient $\chi^{(2)}$. Within the Heisenberg formalism, the propagation
equations describing the evolution of the signal and idler creation
$\tilde{a}_{s,i}^{{\dagger}}(\omega_{0}+\Omega)$ and annihilation
$\tilde{a}_{s,i}(\omega_{0}+\Omega)$ operators can be written as loudon ;
barak

$\displaystyle\frac{\partial\tilde{a}_{s}\left(z,\Omega\right)}{\partial z}$
$\displaystyle=$
$\displaystyle\sigma\tilde{a}_{i}^{{\dagger}}\left(z,-\Omega\right)\exp\left[i\Delta(\Omega)z\right],$
(1) $\displaystyle\frac{\partial\tilde{a}_{i}\left(z,\Omega\right)}{\partial
z}$ $\displaystyle=$
$\displaystyle\sigma\tilde{a}_{s}^{{\dagger}}\left(z,-\Omega\right)\exp\left[i\Delta(\Omega)z\right],$
(2)

where $\Omega$ is the frequency deviation from the central frequency
$\omega_{0}$,
$\Delta\left(\Omega\right)=k_{p}^{0}-k_{s}\left(\Omega\right)-k_{i}\left(-\Omega\right)$
is the phase matching function, $\sigma$ is a constant parameter proportional
to the number of pump photons traversing the nonlinear crystal,
$k_{j}(\Omega)=(\omega_{0}+\Omega)n_{j}/c$ denotes the signal and idler
wavenumbers, with $n_{j}$ being the refractive index at the corresponding
central frequencies, $k_{p}^{0}=2\omega_{0}n_{p}/c$ the wavenumber of the pump
beam and $c$ the speed of light in a vacuum.

Solving Eqs. (1) and (2) and writing $\tilde{a}_{s}(z,\Omega)$ and
$\tilde{a}_{i}(z,\Omega)$ in terms of the vacuum field operators at the input
face of the nonlinear crystal, $b_{s}(\Omega)$ and $b_{i}(\Omega)$, we obtain
navez ; brambilla

$\displaystyle\tilde{a}_{s}(z,\Omega)$ $\displaystyle=$ $\displaystyle
U(z,\Omega)b_{s}(\Omega)+V(z,\Omega)b_{i}^{\dagger}(-\Omega),$ (3)
$\displaystyle\tilde{a}_{i}(z,\Omega)$ $\displaystyle=$ $\displaystyle
U(z,\Omega)b_{i}(\Omega)+V(z,\Omega)b_{s}^{\dagger}(-\Omega)$ (4)

with

$\displaystyle U(z,\Omega)$ $\displaystyle=$
$\displaystyle\exp\left[\frac{i\Delta(\Omega)z}{2}\right]\left\\{\cosh\left[\Gamma(\Omega)z\right]\right.$
(5)
$\displaystyle\left.-\frac{i\Delta(\Omega)}{2\Gamma(\Omega)}\,\sinh\left[\Gamma(\Omega)z\right]\right\\},$
$\displaystyle V(z,\Omega)$ $\displaystyle=$
$\displaystyle-\frac{i\sigma}{\Gamma(\Omega)}\exp\left[i\frac{\Delta(\Omega)z}{2}\right]$
(6) $\displaystyle\times\sinh\left[\Gamma(\Omega)z\right],$

where $\Gamma(\Omega)=(|\sigma|^{2}-\Delta^{2}(\Omega)/4)^{1/2}$ and $\sigma$
is the nonlinear coefficient of the medium.

From these results, we can calculate the required correlations in the temporal
and frequency domains. We will refer to the second-order coherence function
between signal and idler photons as to _interbeam correlations_. Analogously,
we will refer to the second-order coherence function between signal-signal, or
idler-idler, streams of photons as to _intrabeam correlations_. We will make
use of the commutation rules
$[b_{j}(\Omega_{1}),\,b_{k}(\Omega_{2})]=\delta(\Omega_{1}-\Omega_{2})\delta_{jk}$,
where $\delta_{jk}$ is the Kronecker’s delta and $(j,k)=(s,i)$.

In particular, for the effect of dispersion cancelation, two quantities are of
interest. For the interbeam second-order correlations, we have

$G^{(2)}_{\textrm{inter}}(\tau)=\langle
a_{s}^{\dagger}(t)a_{i}^{\dagger}(t+\tau)a_{i}(t+\tau)a_{s}(t)\rangle.$ (7)

Here, the operators $a_{s,i}(t)$ form a Fourier transform pair with the
creation operators $\tilde{a}_{s,i}(\Omega)$. The above magnitude is related
to the probability of detecting a signal photon at time $t$ in coincidence
with an idler photon at $t+\tau$. For the intrabeam correlations, the
corresponding magnitude is

$G^{(2)}_{\textrm{intra}}(\tau)=\langle
a_{s}^{\dagger}(t)a_{s}^{\dagger}(t+\tau)a_{s}(t+\tau)a_{s}(t)\rangle,$ (8)

where signal photons have been arbitrarily chosen. This function provides the
probability of detecting a signal photon in coincidence with another signal
photon at $t+\tau$.

Recent experiments suggest blauensteiner  that intrabeam correlations in SPDC
can be measured when high-power CW beams are used to pump the nonlinear
crystal. Then the intrabeam higher-order correlations show classical-like
features similar to thermal light yurke . In such a case, the interbeam
correlations must be also adequately accounted for to include the multiphoton
pair effects.

For the effects of cancelation of temporal modulation, the relevant second-
order correlation functions are calculated in the spectral domain. In
particular, for the interbeam photon correlations, one has

$\tilde{G}^{(2)}_{\textrm{inter}}(\Omega_{1},\Omega_{2})=\langle\tilde{a}_{s}^{\dagger}(\Omega_{1})\tilde{a}_{i}^{\dagger}(\Omega_{2})\tilde{a}_{i}(\Omega_{2})\tilde{a}_{s}(\Omega_{1})\rangle,$
(9)

which determines the probability of detecting a signal photon at frequency
$\Omega_{1}$ in coincidence with an idler photon at frequency $\Omega_{2}$.
Analogously, we consider the intrabeam second-order correlation function

$\tilde{G}^{(2)}_{\textrm{intra}}(\Omega_{1},\Omega_{2})=\langle\tilde{a}_{s}^{\dagger}(\Omega_{1})\tilde{a}_{s}^{\dagger}(\Omega_{2})\tilde{a}_{s}(\Omega_{2})\tilde{a}_{s}(\Omega_{1})\rangle,$
(10)

which establishes the probability of detecting in coincidence two signal
photons, one at frequency $\Omega_{1}$ and the other at $\Omega_{2}$.

##  III Interbeam configuration: Cancelation of dispersion and temporal
modulation in the multiphoton regime

In this section, we calculate the change in the interbeam correlations,
$G^{(2)}_{\textrm{inter}}(\tau)$ and
$G^{(2)}_{\textrm{inter}}(\Omega_{1},\Omega_{2})$, when the photons propagate
through the optical systems corresponding to dispersion and temporal
modulation cancelations, respectively. The goal of this study is to extend the
previous theoretical results franson1992 ; harris2008  to the case in which
multiphoton pairs are generated in the SPDC process.

In the interbeam configuration, signal and idler photons are separated (see
Fig. 1), following afterwards paths $1$ and $2$. This arrangement can be
achieved by generating SPDC photons that propagate in different directions
(noncollinear type I) or that show orthogonal polarizations (collinear type
II) and are divided by a polarizing beam splitter.

Figure 1:  (Color online) General configurations for observing interbeam
correlations in (a) the temporal domain and (b) the frequency domain. PBS:
Polarizing beam splitter; SPCM: Single photon counting module; cc: coincidence
counting; Mono: Monochromator. In (a) there are two single-mode optical fibers
in paths $1$ and $2$ with transfer functions $H_{1}(\Omega)$ and
$H_{2}(\Omega)$, respectively. In (b), there are two temporal phase modulators
with transfer functions $m_{1}(t)$ and $m_{2}(t)$. Red and blue lines describe
photons with orthogonal polarizations.

###  III.1 Dispersion cancelation in the interbeam configuration

Let us start by considering remote dispersion cancelation. As depicted in Fig.
1(a), before reaching the detectors, the signal photons traverse path $1$,
with dispersive transfer function $H_{1}(\omega)$, and the idler photons
traverse the dispersive medium $2$, described by $H_{2}(\omega)$. In this
case, the annihilation operators for the signal and idler photons at $D_{1}$
and $D_{2}$ are given by $a^{\prime}_{s}(t)=\int\textrm{d}\Omega
H_{1}(\omega_{0}+\Omega)\tilde{a}_{s}(\Omega)\exp(-i\Omega t)$ and
$a^{\prime}_{i}(t)=\int\textrm{d}\Omega
H_{2}(\omega_{0}+\Omega)\tilde{a}_{i}(\Omega)\exp(-i\Omega t)$. Substituting
these expressions into Eq. (7) and taking into account Eqs. (3) and the
commutation rules, it is easy to show that

$\displaystyle G_{\textrm{inter}}^{(2)}(\tau)=N^{2}+\frac{1}{(2\pi)^{2}}$
$\displaystyle\left|\int\textrm{d}\Omega\exp(i\Omega\tau)R(\Omega)H_{1}(\omega_{0}+\Omega)H_{2}(\omega_{0}-\Omega)\right|^{2},$
(11)

where $R(\Omega)=U(\Omega)V(-\Omega)$. From Eq. (III.1), we observe that
$G_{\textrm{inter}}^{(2)}(\tau)$ consists of two terms: The first one is just
a constant background provided by the flux of photon pairs,
$N=N_{s}(t)=N_{i}(t)=\langle a_{s}^{\dagger}(t)a_{s}(t)\rangle$, and the
second contains the temporal structure of the second-order correlation
function that is indeed affected by the spectral transfer functions of the
elements placed in the photon arms.

In the absence of any dispersive medium in the propagation paths of both
streams of photons, one would get

$G^{(2)}_{\textrm{inter}}(\tau)=N^{2}+\left|1/(2\pi)\,\int\textrm{d}\Omega\exp(i\Omega\tau)R(\Omega)\right|^{2}.$
(12)

Thus, the dispersive media only affect the temporal structure of the second-
order correlation function, not the background term, which is always present.
For the particular case in which the media can be assumed to be first-order
dispersive, such as single-mode fibers, we have
$H_{1}(\Omega)=\exp(i\Phi_{1}\Omega^{2}/2)$ and
$H_{2}(\omega)=\exp(i\Phi_{2}\Omega^{2}/2)$, with the GDD parameters
$\Phi_{k}=\beta_{2k}L_{k}$, $k=(1,2)$. $\beta_{2k}$ is the group velocity
dispersion (GVD) coefficient of fiber $k$ and $L_{k}$ its length. Therefore,
the dispersion effect can be suppressed if $\Phi_{1}=-\Phi_{2}$, i.e., if the
GDD parameters in the two paths have opposite signs franson1992 . The type of
frequency correlation between signal and idler photons, i.e. frequency
anticorrelation $\Omega_{s}+\Omega_{i}=0$, also implies that odd-order
dispersion terms cannot be canceled, but on the contrary, their effects are
added.

The effect of increasing the pair generation rate is to degrade the signal-to-
background ratio. With the aid of Eqs. (5) and (6), it can be shown that this
ratio scales as $1/(NL)$ and tends to infinity as N decreases harris2005 .

###  III.2 Temporal modulation cancelation in the interbeam configuration

The configuration to study remote temporal modulation cancelation is depicted
in Fig. 1(b). The signal photons travel through the modulator $m_{1}(t)$ and
the idler photons through modulator $m_{2}(t)$. In order to measure frequency
correlations, signal and idler photons pass through ideal scanning
monochromators before reaching the photodetectors harris2008 ; harris2009 .

For simplicity, we consider the photons to be phase modulated with sinusoidal
phase modulators with temporal complex transfer function
$m_{1,2}(t)=\exp[i\Delta\theta_{1,2}\sin(\Omega_{m}t)]$, where $\Omega_{m}$ is
the modulation frequency (in the microwave regime) and $\Delta\theta_{1,2}$
are the modulation indexes of the corresponding modulators. In order to
calculate the interbeam second-order correlation function,
$\tilde{G}^{(2)}_{\textrm{inter}}(\Omega_{1},\Omega_{2})$, we need to
calculate the evolution of the operators $\tilde{a}_{s,i}(\Omega)$ after phase
modulation. This is done by making the convolution of the operators with the
frequency response of the corresponding modulation function $m_{1,2}(t)$
harris2008 .

In the frequency domain, the action of a sinusoidal phase modulator can be
written as $M(\Omega)=\sum_{n}J_{n}(\Delta\theta)\delta(\Omega-n\Omega_{m})$,
where $J_{n}$ are the Bessel functions of the first kind. Thus one obtains
that

$\tilde{G}_{\textrm{inter}}^{(2)}(\Omega_{1},\Omega_{2})=N_{s}(\Omega_{1})N_{i}(\Omega_{2})+\left|T(\Omega_{1},\Omega_{2})\right|^{2},$
(13)

with

$\displaystyle
T\left(\Omega_{1},\Omega_{2}\right)=\frac{1}{2\pi}\int\textrm{d}\Omega
R(\Omega)$ $\displaystyle\times
M_{1}(\Omega_{1}+\Omega)M_{2}(\Omega_{2}-\Omega).$ (14)

In an analogous way to the dispersive case, there is a constant background
term proportional to the flux of photons at the considered frequencies and
another intricate non-factorizable term that gets affected by the spectral
transfer functions of the modulators.

State-of-the-art electro-optic phase modulators can provide a maximum optical
bandwidth of tens of GHz, yet much smaller than the bandwidth of most common
SPDC sources with typical values around $10$-$20$ THz. Under these conditions,
we can safely write

$\displaystyle T(\Omega_{1},\Omega_{2})=R(\Omega_{-}/2)$
$\displaystyle\times\sum_{n=-\infty}^{\infty}J_{n}(\Delta\theta_{1}+\Delta\theta_{2})\delta(\Omega_{+}-n\Omega_{m}),$
(15)

where $\Omega_{-}=\Omega_{1}-\Omega_{2}$ and
$\Omega_{+}=\Omega_{1}+\Omega_{2}$.

Before proceeding further, let us write the corresponding expression in the
absence of temporal modulation:

$\displaystyle\tilde{G}_{\textrm{inter}}^{(2)}(\Omega_{1},\Omega_{2})=N_{s}(\Omega_{1})N_{i}(\Omega_{2})+$
$\displaystyle\left|R(\Omega_{1})\right|^{2}\delta(\Omega_{1}+\Omega_{2}).$
(16)

This equation indicates that, besides a background term, there is a strong
correlation between frequencies satisfying $\Omega_{1}=-\Omega_{2}$. From Eq.
(III.2), we conclude that it is possible to recover this result if the
modulation depths of the two modulators fulfill the condition
$\Delta\theta_{1}=-\Delta\theta_{2}$. In this case, the only non-zero Bessel
term is $n=0$, so that
$T(\Omega_{1},\Omega_{2})=R(\Omega_{1})\delta(\Omega_{1}+\Omega_{2})$, i.e.,
the effect of the modulators is canceled harris2008 . Again, the effect of
having multiphoton pairs is only to degrade the signal-to-background ratio,
which decreases with the flux-rate of the generated pairs.

We remark that both dispersion and temporal modulation cancelations effects
are due to the presence of frequency anti-correlation between signal and idler
photons, which can be expressed by the basic relationship
$\langle\tilde{a}_{s}(\Omega_{1})\tilde{a}_{i}(\Omega_{2})\rangle=R(\Omega_{1})\delta(\Omega_{1}+\Omega_{2})$.

##  IV Intrabeam configuration: Remote cancelation of dispersion and temporal
modulation with non-entangled photon pairs

Let us now consider the intrabeam configuration depicted in Fig. 2. Without
loss of generality, one of the beams produced by the source is discarded and
the attention is centered on one of the beams only, for example, the signal.

Figure 2:  (Color online) General configurations for observing intrabeam
correlations in (a) the temporal domain and (b) the frequency domain. PBS:
Polarizing beam splitter; BS: Beam Splitter; SPCM: Single photon counting
module; cc: coincidence counting; Mono: Monochromator. In (a) there are two
single-mode fibers in paths $1$ and $2$ with transfer functions
$H_{1}(\Omega)$ and $H_{2}(\Omega)$, respectively. In (b), there are two
temporal phase modulators with transfer functions $m_{1}(t)$ and $m_{2}(t)$.
Red and blue lines describe photons with orthogonal polarizations. Photons
with a given polarization are discarded.

###  IV.1 Remote dispersion compensation in the intrabeam configuration

Let us consider the situation depicted in Fig. 2(a). After traversing the
dispersive media located in paths $1$ and $2$, the probability to detect a
photon in path $1$ at time $t$ in coincidence with a photon in path $2$ at
time $t+\tau$ is given by the intrabeam second-order correlation function,

$G^{(2)}_{\textrm{intra}}(\tau)=\langle
a^{\dagger}_{s1}(t)a^{\dagger}_{s2}(t+\tau)a_{s2}(t+\tau)a_{s1}(t)\rangle.$
(17)

By analogy to the case of the previous section, the corresponding operators
are calculated as
$a_{s(1,2)}(t)=\int\textrm{d}\Omega\tilde{a}_{s}(\Omega)H_{1(2)}(\Omega)\exp(-i\Omega
t)$, where $H_{1(2)}(\Omega)$ denotes the complex transfer function placed in
the path of photon 1(2). Substituting these expressions and taking into
account Eqs. (5) and (6), it is easy to show that

$\displaystyle
G_{\textrm{intra}}^{(2)}(\tau)=N^{2}+\frac{1}{(2\pi)^{2}}\left|\int\textrm{d}\Omega\,S(\Omega)\right.$
$\displaystyle\left.\times
H_{1}^{*}(\omega_{0}+\Omega)H_{2}(\omega_{0}+\Omega)\exp\left(i\Omega\tau\right)\right|^{2},$
(18)

where $S(\Omega)=|V(\Omega)|^{2}$. We observe that the joint probability
detection contains the same background term as in the interbeam case. Even
more, the functional dependence on $\tau$ depends on the multiplication of the
transfer functions of the dispersive media in such a way, that their phase
difference may alter completely the output. However, now the previous role of
the spectrum $R(\Omega)=U(\Omega)V(-\Omega)$ is replaced by the signal’s
photon spectrum $S(\Omega)$.

For the sake of comparison, we derive the expression of
$G^{(2)}_{\textrm{intra}}(\tau)$ in the absence of dispersive media:

$\displaystyle
G_{\textrm{intra}}^{(2)}(\tau)=N^{2}+\left|\Gamma(\tau)\right|^{2},$ (19)

where $\Gamma(\tau)=1/(2\pi)\,\int\textrm{d}\Omega\exp(i\Omega\tau)S(\Omega)$.
Notice that a similar suppression of the effects of dispersion can be obtained
if both dispersive media are identical, i.e., $H_{1}\equiv H_{2}$. In contrast
to the interbeam case, the suppression of the dispersive effects is not
limited to even-order dispersion terms, but it affects all the terms
victor2009 . On the other hand, the signal-to-background ratio is maximally
bounded to $2$, thus challenging the observation of the effect in an
experiment. The factor of $2$ is a manifestation of the thermal-like character
of the intrabeam correlations. Note that a reduced signal-to-background ratio
could also be attained in the interbeam correlations provided that sufficient
multiple pairs were generated in the SPDC process.

###  IV.2 Remote temporal modulation compensation in the intrabeam
configuration

Finally, we revisit the temporal modulation cancelation scheme with the
photons from only one of the downconverted beams (signal, for example), as
depicted in Fig. 2(b). After modulation, the probability to detect a photon in
path $1$ with frequency $\omega_{0}+\Omega_{1}$ in coincidence with a photon
in path $2$ with frequency $\omega_{0}+\Omega_{2}$ is provided by the
magnitude

$\tilde{G}^{(2)}_{\textrm{intra}}(\Omega_{1},\Omega_{2})=\langle\tilde{a}_{s1}^{\dagger}(\Omega_{1})\tilde{a}_{s2}^{\dagger}(\Omega_{2})\tilde{a}_{s2}(\Omega_{2})\tilde{a}_{s1}(\Omega_{1})\rangle.$
(20)

To calculate the evolution of the operators, one just has to proceed as in
section III(B), i.e., to calculate the convolution of the signal operator with
the modulator transfer function in the spectral domain. Then

$\displaystyle\tilde{G}_{\textrm{intra}}^{(2)}(\Omega_{1},\Omega_{2})$
$\displaystyle=$ $\displaystyle N_{s}(\Omega_{1})N_{s}(\Omega_{2})$ (21)
$\displaystyle+\left|T^{\prime}(\Omega_{1},\Omega_{2})\right|^{2},$

where

$\displaystyle
T^{\prime}(\Omega_{1},\Omega_{2})=\frac{1}{2\pi}\int\textrm{d}\Omega^{\prime}S(\Omega^{\prime})$
$\displaystyle\times
M_{1}^{*}(\Omega_{1}-\Omega^{\prime})M_{2}(\Omega_{2}-\Omega^{\prime}).$ (22)

As before, we obtain a background term proportional to the number of photons
in modes $\omega_{0}+\Omega_{1}$ and $\omega_{0}+\Omega_{2}$. The second term
is a non-separable function in frequency which is affected by the transfer
functions of the modulators. Analogously to the derivation of Eq. (III.2), we
can consider low-bandwidth modulators, so that this function reduces to

$\displaystyle T^{\prime}(\Omega_{1},\Omega_{2})=S(\Omega_{+}/2)$
$\displaystyle\times\sum_{n=-\infty}^{\infty}J_{n}(\Delta\theta_{1}-\Delta\theta_{2})\delta(\Omega_{-}-n\Omega_{m}).$
(23)

To close the loop, we calculate the second-order correlation function that
would be achieved in the absence of modulators,

$\displaystyle\tilde{G}_{\textrm{intra}}^{(2)}(\Omega_{1},\Omega_{2})=N_{s}(\Omega_{1})N_{s}(\Omega_{2})+$
$\displaystyle\left|S(\Omega_{1})\right|^{2}\delta(\Omega_{1}-\Omega_{2}).$
(24)

Similarly to the interbeam case, there is a strong correlation between photons
in paths $1$ and $2$ with the same frequency. Giving a further step, from Eq.
(IV.2), if $\Delta\theta_{1}=\Delta\theta_{2}$, i.e., photons in paths $1$ and
$2$ are equally phase modulated, the effects of the temporal phase modulation
are suppressed. Now, the only non-zero term is $n=0$, and therefore only the
detection of photons with frequencies $\Omega_{1}=\Omega_{2}$ is enhanced.

The origin of these suppression (cancelation) effects when the photons that
follow paths $1$ and $2$ are equally phase-modulated, or traverse equal
dispersive media, is the existence of frequency correlation between the signal
photons, i.e.,
$\langle\tilde{a}_{s}^{\dagger}(\Omega_{1})\tilde{a}_{s}(\Omega_{2})\rangle=S(\Omega_{2})\delta(\Omega_{2}-\Omega_{1})$.
Eqs. (IV.1) and (IV.2), which allow for the remote dispersion and temporal
modulation cancelation, make an explicit use of the existence of these
characteristic frequency correlations. This strong frequency correlation could
be revealed measuring, before any temporal phase modulation takes place, the
number of photons in path $1$ with frequency $\omega_{0}+\Omega_{1}$ in
coincidence with the signal photons with frequency $\omega_{0}+\Omega_{2}$
that traverse path $2$, which is given by Eq. (IV.2).

##  V Summary and conclusions

We have considered the effects of cancelation of dispersion and temporal
modulation in the regime of multiphoton pair generation in the SPDC process.
The two schemes are studied with two kinds of photon correlations, those
arising from the downconverted signal and idler photons (interbeam
correlations), which show entanglement, and those from the individual photons
in a single beam (intrabeam correlations), which show strong correlations but
not entanglement. In the intrabeam regime, it is possible to achieve the
suppression of the effects of either dispersion or temporal modulation at all
orders.

An important point to remark is that the observation of the remote cancelation
of the dispersion and modulation effects in the interbeam configuration
happens even in the high-flux regime ($\sigma L\gg 1$). In this case, an
inequality of normalized second-order correlation functions howell2006 ;
mancini2002  which is fulfilled by classical-like fields, is no longer
violated. This highlights the role of the frequency anticorrelation and
correlation effects as the reason for the observation of any dispersion or
modulation cancelation effects.

The observation in a particular setting of dispersion and modulation
cancelation effects depends on the type of frequency correlations present. In
the interbeam case, where the photons show frequency anticorrelation, the
observation of dispersion and modulation cancelation effects requires that
both photons suffer only even-order dispersion of the opposite sign, or
identical phase modulation but with opposite phases.

On the other hand, if the photons are frequency correlated (intrabeam
correlations with $\Omega_{1}-\Omega_{2}=0$), the observation of dispersion
and modulation cancelation effects requires that both photons suffer equal
dispersion, for all dispersion terms, or identical phase modulation, with the
same phases.

The effects described here bear important similarities with the description of
two-photon imaging experiments with two different types of two-photon sources:
paired photons entangled in the spatial degree of freedom and bunched pairs of
photons coming from a thermal light source scarcelli2004 . In pittman1995 , a
two-photon optical imaging experiment was performed based on the spatial
correlations of the signal and idler photon pairs produced in SPDC. The
possibility to use thermal (or pseudothermal) radiation for two-photon imaging
experiments has also been demonstrated valencia2005 . In both cases, the
important point that enables the observation of similar effects with
dissimilar sources is the presence of spatial correlations between the paired
photons.

Again, the spatial correlations between signal and idler photons show
different characteristics than the spatial correlations of signal pairs. But
the observation of two-photon imaging, as well as the capacity to cancel
diffraction effects when measuring two-photon coincidences, are due to some of
the common characteristics that both types of sources share: the presence of
certain correlations between paired photons.

## Acknowledgements

This work was supported by the Government of Spain (Consolider Ingenio
CSD2006-00019, FIS2010-14831), and was supported in part by FONCICYT project
94142\. The project PHORBITECH acknowledges the financial support of the
Future and Emerging Technologies (FET) programme within the Seventh Framework
Programme for Research of the European Commission, under FET-Open grant
number: 255914.

## References

  * (1) R. J. Glauber, Phys. Rev. 131, 2766 (1963).
  * (2) R. J. Glauber, Rev. Mod. Phys. 78, 1267 (2006).
  * (3) L. Mandel and E. Wolf, Optical coherence and quantum optics, Cambridge University Press, 1995.
  * (4) A. Kuzmich, W. P. Bowen, A. D. Boozer, A. Boca, C. W. Chou, L.-M. Duan, H. J. Kimble, Nature 423, 731 (2003).
  * (5) R. Loudon, Quantum Theory of Light, Oxford University Press, 1st edition 1973, 3rd edition 2000.
  * (6) J. D. Franson, Phys. Rev. A 45, 3126 (1992).
  * (7) J. Brendel, H. Zbinden and N. Gisin, Opt. Comm. 151, 35 (1998).
  * (8) S. Y. Baek, Y. W. Cho and Y. H. Kim, Opt. Express 17, 19244 (2009).
  * (9) A. Valencia, M. V. Chekhova, A. Trifonov and Y. Shih, Phys. Rev. Lett. 88, 183601 (2002).
  * (10) S. E. Harris, Phys. Rev. A 78, 021807 (2008).
  * (11) S. Sensarn, G. Y. Yin, and S. E. Harris, Phys. Rev. Lett. 103, 163601 (2009).
  * (12) V. Torres-Company, H. Lajunen, and A. T. Friberg, New J. Phys. 11, 063041 (2009).
  * (13) J. H. Shapiro, Phys. Rev. A 81 023824 (2010).
  * (14) J. D. Franson, Phys. Rev. A 80, 032119 (2009).
  * (15) J. D. Franson, Phys. Rev. A 81, 023825 (2010).
  * (16) B. Dayan, Phys. Rev. A 76, 043813 (2007).
  * (17) P. Navez, E. Brambilla, A. Gatti, and L. A. Lugiato, Phys. Rev. A 65, 013813 (2001).
  * (18) E. Brambilla, A. Gatti, M. Bache, and L. A. Lugiato, Phys. Rev. A 69, 023802 (2004).
  * (19) B. Blauensteiner, I. Herbauts, S. Bettelli, A. Poppe, and H. Hübel, Phys. Rev. A 79, 063846 (2009).
  * (20) B. Yurke and M. Potasek, Phys. Rev. A 36, 3464 (1987).
  * (21) V. Balic, D. A. Braje, P. Kolchin, G. Y. Yin, and S. E. Harris, Phys. Rev. Lett. 94, 183601 (2005).
  * (22) I. A. Khan and J. C. Howell, Phys. Rev. A 73, 031801(R) (2006).
  * (23) S. Mancini, V. Giovannetti, D. Vitali, and P. Tombesi, Phys. Rev. Lett. 88, 120401 (2002).
  * (24) G. Scarcelli, A. Valencia, and Y. Shih, Phys. Rev. A 70 051802 (2004).
  * (25) T. B. Pittman, Y. Shih, D. V. Strekalov, and A. V. Sergienko, Phys. Rev. A, 52, 3429(R) (1995).
  * (26) A. Valencia, G. Scarcelli, M. D Angelo, and Y. Shih, Phys. Rev. Lett. 94 063601 (2005).


# Independent Control of Two Magnetic Robots using External Permanent Magnets:
A Feasibility Study

Joshua Davy, Tomas da Veiga, Giovanni Pittiglio, James H. Chandler,
and Pietro Valdastri Research reported in this article was supported by the
Engineering and Physical Sciences Research Council (EPSRC) under grants number
EP/R045291/1, EP/V009818/1 and EP/V047914/1, and by the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement No 818045). Any opinions, findings and conclusions,
or recommendations expressed in this article are those of the authors and do
not necessarily reflect the views of the EPSRC or the ERC.For the purpose of
open access, the authors have applied a Creative Commons Attribution (CCBY)
license to any Accepted Manuscript version arising. Figures 1 and 2 created
with BioRender.com Joshua Davy, Tomas da Veiga, James H. Chandler, and Pietro
Valdastri are with the STORM Lab, Institute of Autonomous Systems and Sensing
(IRASS), School of Electronic and Electrical Engineering, University of Leeds,
Leeds, UK. Email: {el17jd, eltgdv, j.h.chandler<EMAIL_ADDRESS>Giovanni Pittiglio is with the Department of Cardiovascular Surgery, Boston
Children’s Hospital, Harvard Medical School, Boston, MA 02115, USA. Email:
<EMAIL_ADDRESS>
###### Abstract

The ability to have multiple magnetic robots operate independently in the same
workspace would increase the clinical potential of these systems allowing
collaborative operation. In this work, we investigate the feasibility of
actuating two magnetic robots operating within the same workspace using
external permanent magnets. Unlike actuation systems based on pairs of
electromagnetic coils, the use of multiple permanent magnets comes with the
advantage of a large workspace which better suits the clinical setting. In
this work, we present an optimization routine capable of generating the
required poses for the external magnets in order to control the position and
orientation of two magnetic robots. We show that at a distance of 15cm,
minimal coupling between the magnetic robots can be achieved (3.9% crosstalk)
each embedded with 5mm diameter, 5mm length NdFeB magnets. At smaller
distances, we observe that the ability to independently control the robot
torques decreases, but forces can still achieve independent control even with
alignment of the robots. We test our developed control system in a simulation
of two magnetic robots following pre-planned trajectories in close proximity
(60 mm) showing a mean positional error of 8.7 mm and mean angular error of
16.7∘.

##  I Introduction

Magnetic actuation has been utilized for surgical robotic systems due to its
potential for miniaturization and tetherless control [1]. Magnetic robots can
be miniaturized to the millimeter and, more recently, the micrometer scale
[2], [3] without loosing controllable Degrees of Freedom (DOF). Additionally,
safety concerns with other robotic actuation methods (such as pressurized air)
are not present with magnetic actuation [4]. This ability to guide magnetic
robots through natural orifices allows them to reach deep within the anatomy
without the invasiveness of incisions.

Controlling multiple magnetic robots within a shared workspace may provide new
opportunities for minimally invasive surgical interventions (see Fig. 1).
Magnetic robots have been equipped with surgical tools such as forceps,
cameras, and biopsy needles [2], [5], [6]. The ability to have multiple
magnetic robots work in collaborative operation would provide more options for
surgeons in interventions by allowing multiple tools to be used in parallel.
Furthermore, magnetic robots could be deployed in various parts of the anatomy
(for example, along the gastrointestinal tract) and perform simultaneous
microbiome sampling [7], [8] or drug delivery [9].

Figure 1: Overview of our proposed control scenario. We aim to independently
control the full five degrees of freedom of (each internal robot each
containing a single IPM). Multiple EPMs are manipulated in order to produce
the desired magnetic field. The relative orientation and distance between
robots affects the coupling of the robot control.

Independent control of magnetic robots has been extensively explored using
coil-based electromagnetic systems [10], [11], [12], [13], [14], [15], [16].
The limited workspace of electromagnetic coil systems means the magnetic field
between robots is coupled and therefore cannot be independently controlled
[17]. Instead, various techniques have been employed to control the individual
magnetic wrench on the robots. This has included the design of heterogeneous
robots whose wrench response to the magnetic field differs, the use of
underactuated control where robots are positioned and oriented with limited
manipulability, and the exploitation of inter-agent magnetic interaction [12],
[13], [15].

By varying the external magnetic field, the resultant magnetic wrench on the
robot can be controlled. This external field is usually provided by either a
system of stationary magnetic coils or manipulated permanent magnets [6],
[18]. Actuation systems based on electromagnetic coils have a linear
relationship between applied electric currents and the generated magnetic
field and can be arranged to produce fields with constant gradients across the
workspace leading to a simplified control strategy. However, many designs of
coil systems suffer from limited workspace size and large power supply
requirements to generate strong magnetic fields [19], [20]. Alternatively,
manipulated permanent magnets offer the advantage of low power requirements
and large workspace but produce non-linear magnetic fields. This leads to a
more complex control scenario where the relationship between magnet poses and
induced wrench on the robot is non-trivial [18], [21], [22].

Magnetic agents will experience an aligning torque with the external field and
a force relating to the external magnetic field gradient. Precise control of
this external field allows for a maximum of five Degrees of Freedom (DOFs) for
a magnetic agent (as the torque aligned with the axis of the direction of
magnetism cannot be controlled magnetically) [19]. In a current-free workspace
the number of DOF at a singular point in the workspace is limited to eight
[15], [23]. Our previous work has shown how a single External Permanent Magnet
(EPM) mounted to a serial manipulator is capable of independent control of all
5 DOF of a capsule robot equipped with a Internal Permanent magnet (IPM) [2].
Further, we have shown the independent control of eight DOFs in the case of
orthogonal IPMs assuming a constant magnetic gradient across the agents with
two EPMs [21] [18]. An individual EPM can fully control the five DOF of a
single IPM; in this work we show that at a sufficient distance, where the
interaction between each robot is small, then two EPMs could independently
control the full ten DOF of two IPMs. However, as the IPMs are brought
together, the inter-agent interaction between IPMs and the combined magnetic
field generated by multiple EPMs affects this independence. Understanding this
is key to the feasibility of operating multiple magnetic robots with full
control in the same workspace, such as within the constraints of the human
anatomy.

In this work, we analyse the case of attempting to control two magnetic robots
each containing a single IPM operating in close proximity with EPMs (see Fig.
1). We formulate an optimization strategy capable of calculating the required
poses of the EPMs for a desired wrench on the robots and further analyze the
independence of control between the robots, evaluating our solutions on a
dynamic simulation with a closed-loop controller based on positional and
orientation feedback.

##  II Independent Actuation

Figure 2: The combined magnetic fields of multiple EPMs create highly non-
linear magnetic fields, that vary greatly in direction, magnitude and gradient
across the workspace where the two IPMs are. The separation between IPMs means
the local magnetic field cannot be treated linearly. Our optimization strategy
instead relies on directly solving for the wrench of the IPMs.

###  II-A Magnetics

We consider the control of two medical capsule robots each embedded with a
single IPM actuated under the influence of $N$ EPMs. Via manipulation of the
pose of the EPMs, the local magnetic field acting on the IPMs can be
controlled. This local magnetic field induces a resultant magnetic wrench on
the IPMs.

An IPM with magnetic moment $\mathbf{m}$ under the influence of an external
magnetic field will experience a torque with an external field $\mathbf{B}$

$\bm{\tau}=\mathbf{m}\times\mathbf{B}.$ (1)

Further, in the case of a non-homogeneous magnetic field (i.e. a field where
magnetic field gradients are present) a force is induced

$\textbf{F}=\textbf{B}_{\nabla}^{T}\textbf{m},$ (2)

where $\textbf{B}_{\nabla}$ is the spatial derivative of the magnetic field or
field Jacobian

$\mathbf{B}_{\nabla}=\left(\begin{array}[]{ccc}\frac{\partial B_{x}}{\partial
x}&\frac{\partial B_{x}}{\partial y}&\frac{\partial B_{x}}{\partial z}\\\
\frac{\partial B_{y}}{\partial x}&\frac{\partial B_{y}}{\partial
y}&\frac{\partial B_{y}}{\partial z}\\\ \frac{\partial B_{z}}{\partial
x}&\frac{\partial B_{z}}{\partial y}&\frac{\partial B_{z}}{\partial
z}\end{array}\right).$ (3)

Assuming electro-static contributions to be negligible, and therefore the
workspace to be current-free, Maxwells equations state that the field Jacobian
is traceless and symmetric [10]. This limits the independent components of
this to five

$\mathbf{B}_{\nabla}=\left(\begin{array}[]{ccc}\frac{\partial B_{x}}{\partial
x}&\frac{\partial B_{x}}{\partial y}&\frac{\partial B_{x}}{\partial z}\\\
\frac{\partial B_{x}}{\partial y}&\frac{\partial B_{y}}{\partial
y}&\frac{\partial B_{y}}{\partial z}\\\ \frac{\partial B_{x}}{\partial
z}&\frac{\partial B_{y}}{\partial z}&-\frac{\partial B_{x}}{\partial
x}-\frac{\partial B_{y}}{\partial y}\end{array}\right).$ (4)

This, along with the three independent components of the magnetic field
$\mathbf{B}\in\mathbb{R}^{3}$, forms the eight independent components of the
magnetic field at a singular point in the workspace. In our previous work
[18], [21], we assumed the separation of IPMs to be negligible and therefore
the magnetic field and gradient to be constant across them. With greater
separation of the IPMs, this assumption of linearity no longer holds.

Under the superposition principle and assuming no other source of magnetic
field, the magnetic field at an IPM $j$ at position $\mathbf{p}_{j}$ in the
workspace is the summation of the magnetic field contribution from each
individual EPM and the other IPM $k$. The magnetic field at $\mathbf{p}_{j}$
is

$\mathbf{B}_{j}(\mathbf{p}_{j})=\sum_{i=1}^{N}\mathbf{B}_{i}(\mathbf{m}_{i},\mathbf{p}_{i},\mathbf{p}_{j})+\mathbf{B}_{k}(\mathbf{m}_{k},\mathbf{p}_{k},\mathbf{p}_{j})$
(5)

where $N$ is the number of EPMs. $\mathbf{B}_{i}$ is the magnetic field
contribution from EPM $i$ with magnetic moment $\mathbf{m}_{i}$ and position
$\mathbf{p}_{i}$. (see Fig. 2). Equivalently $\mathbf{B}_{k}$ is the magnetic
field contribution from the other IPM with magnetic moment $\mathbf{m}_{k}$
and position $\mathbf{p}_{k}$, with $k\neq j$.

Equivalently the magnetic field gradient at position $\mathbf{p}_{j}$ is

$\mathbf{B}_{\nabla,j}(\mathbf{p}_{j})=\sum_{i=1}^{N}\mathbf{B}_{\nabla,i}(\mathbf{m}_{i},\mathbf{p}_{i},\mathbf{p}_{j})+\mathbf{B}_{\nabla,k}(\mathbf{m}_{k},\mathbf{p}_{k},\mathbf{p}_{j})$
(6)

where $\mathbf{B}_{\nabla,i}$ is the magnetic field gradient contribution from
EPM $i$ and $\mathbf{B}_{\nabla,k}$ is the contribution from the other IPM.

From (1) and (2) the induced wrench in the global reference on a IPM is

$\mathbf{w}_{j}^{+}=\left(\begin{matrix}\mathbf{m}_{j}\times\mathbf{B}(\mathbf{p}_{j})\\\
\mathbf{B}_{\nabla}(\mathbf{p}_{j})^{T}\mathbf{m}_{j}\end{matrix}\right)$ (7)

where
$\mathbf{w_{j}}^{+}=(\tau_{j,x}^{+}\;\tau_{j,y}^{+}\;\tau_{j,z}^{+}\;F_{j,x}^{+}\;F_{j,y}^{+}\;F_{j,z}^{+})^{T}$,
represented in the global reference frame. The DOF around the axis of IPM
magnetization cannot be controlled therefore, in order to analyze our system
we observe the wrench on the IPMs individually, in a reference frame where the
axis of magnetization aligns with the IPM $x$ axis. As in this rotated
reference frame the torque around the $x$ axis will always be zero the wrench
vector can be represented as the five controllable DOF for a single IPM
$\mathbf{w}_{j}=(\tau_{j,y}\;\tau_{j,z}\;F_{j,x}\;F_{j,y}\;F_{j,z})^{T}$,
where $j\in\\{1,2\\}$ is the index of the IPM.

###  II-B Dipole Modelling

Assuming the distances involved between EPMs and IPMs to be large enough that
neglecting higher order terms is an acceptable assumption, we utilise the
dipole model for modelling the magnetic fields [24]. The field produced by a
magnetic dipole at a point $\mathbf{p}_{b}$ in the workspace is

$\mathbf{B}(\mathbf{m},\mathbf{p}_{a},\mathbf{p}_{b})=\left(\frac{\mu_{0}}{4\pi\left\|\mathbf{r}\right\|^{3}}(3\hat{\mathbf{r}}\hat{\mathbf{r}}^{T}-\mathbb{I})\right)\mathbf{m},$
(8)

where $\mathbf{r}$ = $\mathbf{p}_{a}-\mathbf{p}_{b}$ is the relative
displacement between the dipole position $\mathbf{p}_{a}$ and measurement
point $\mathbf{p}_{b}$. $\hat{\mathbf{r}}$ is the direction vector of
$\mathbf{r}$, $\hat{\mathbf{r}}$ = $\frac{\mathbf{r}}{||\mathbf{r}||}$.
$\mu_{0}$ is the vacuum permeability equal to $4\pi\times 10^{-7}$ Hm-1.
$\mathbb{I}$ is the identity matrix.

The spatial gradient of the magnetic field produced by a dipole is

$\begin{split}\mathbf{B}_{\nabla}(\mathbf{m},\mathbf{p}_{a},\mathbf{p}_{b})=\frac{3\mu_{0}}{4\pi\left\|{\mathbf{r}}\right\|^{4}}(\mathbf{m}\hat{\mathbf{r}}^{T}+\hat{\mathbf{r}}\mathbf{m}^{T}+\\\
\hat{\mathbf{r}}^{T}\mathbf{m}(\mathbb{I}-5\hat{\mathbf{r}}\hat{\mathbf{r}}^{T}))\end{split}$
(9)

[19].

###  II-C Optimization

Our aim is to independently control the five DOFs of both IPMs. We group these
wrenches into a ten dimensional vector $\mathbf{U}$ =
$(\mathbf{w}_{1}\;\mathbf{w}_{2})^{T}$. For a desired $\mathbf{U}_{d}$ the
optimum EPMs pose must be found. We formulate this as an optimization problem
where we aim to minimize the norm error between desired and current wrench

$\underset{\mathbf{x}}{\text{argmin}}||\mathbf{U}_{d}-\mathbf{U}(\mathbf{x})||^{2}$
(10)

where $\mathbf{x}$ is the set of EPM poses which we define in terms of
spherical coordinates for both position and orientation,

$\mathbf{x}=(\alpha_{1}\;\beta_{1}\;r_{1}\;a_{1}\;b_{1}\;...\;\alpha_{N}\;\beta_{N}\;r_{N}\;a_{N}\;b_{N})^{T}.$
(11)

where $\alpha_{i},\beta_{i},r_{i}$ represent the azimuth, elevation and radius
of the spherical position

$\mathbf{p}_{i}=\left(\begin{matrix}x\\\ y\\\
z\end{matrix}\right)=r_{i}\left(\begin{matrix}\cos(\beta_{i})\cos(\alpha_{i})\\\
\cos(\beta_{i})\sin(\alpha_{i})\\\ \sin(\beta_{i})\end{matrix}\right)$ (12)

and $a_{i},b_{i}$ represents the orientation of the EPM

$\mathbf{m}_{i}=||\mathbf{m}||\left(\begin{matrix}cos(b_{i})\cos(a_{i})\\\
\cos(b_{i})\sin(a_{i})\\\ \sin(b_{i})\end{matrix}\right)$ (13)

where $||\mathbf{m}||$ is the magnetic moment norm of the EPM.

Compared to defining magnetic moment and position in the optimization with
Cartesian coordinates, utilising polar representation reduces problem
dimensionality by allowing the magnetic moment to be defined with only two
variables. Further, a minimum radius around the central point of the workspace
$r_{min}$, can be applied linearly rather then via additional non-linear
constraints.

We utilise the interior point algorithm with the implementation made available
in the MATLAB (Mathworks, USA) programming language function fmincon. As
solutions are dependent on the initialization point, we use a random
initialization from within the defined search space of EPM poses and repeat
for a maximum ($M$ = 10) attempts; selecting the most optimum solution with
the lowest error.

##  III System Analysis

(a) Activation of individual DOF with IPMs orthogonal with a separation of 15
cm ($\theta$ = $90^{\circ}$, $d$ = $0.15$ m) with two EPMS.

(b) Activation of individual DOF with IPMs parallel with a separation of 15 cm
($\theta$ = $0^{\circ}$, $d$ = $0.15$ m) with two EPMs.

(c) Mean Crosstalk across all independent activations of DOF with various
seperations with increasing relative angle $\theta$.

(d) Activation of individual DOF with IPMs orthogonal with a separation of 5
cm ($\theta$ = $90^{\circ}$, $d$ = $0.05$ m) with two EPMs.

(e) Activation of individual DOF with IPMs parallel with a separation of 5 cm
($\theta$ = $0^{\circ}$, $d$ = $0.05$ m) with two EPMs.

(f) Mean Crosstalk across all independent activations of DOF with various
relative orientations $\theta$ with increasing separation $d$.

Figure 3:

A variety of magnet sizes and shapes are utilized within magnetic robots that
navigate the gastrointestinal tract. In our analysis, we utilise cylindrical
magnets made from N52 grade NdFeB of $5$ mm diameter and $5$ mm length to
represent the IPMs of the magnetic robots. By varying the distance $d$ and
relative orientation $\theta$ between the IPMs we expect to see the
independence of the their DOFs vary. We also consider the number of EPMs to
control the system and in order to make a fair comparison we utilize a fixed
volume of magnetic material split between them. We choose this to be the
equivalent of the volume of the two EPMs used in our previous work ($100$ mm
diameter $100$ mm length cylindrical N52 grade NdFeB magnets) [22]. In our
previous experience, these magnets represent a rough upper limit on the size
of EPMs that can be manipulated with standard collaborative manipulators
without causing electromagnetic interface with the manipulator joint encoders.
The distance between EPMs is constrained in order to prevent their mutual
attraction overpowering the manipulators. To do so, we set the minimal
distance between EPMs to ensure their mutual attraction does not exceed $10g$
N, which represents a typical maximum payload of standard collaborative
manipulators $\approx 10$ kg ($g$ = $9.81$ ms-2). We also impose a minimum
radius $r_{min}$ of $0.15$ m between EPMs and the center of the IPMs (see Fig.
2). We aim to independently actuate the forces and torques on each IPM to
values of $0.5$ N and $0.05$ Nm-1 respectively whilst minimizing the
activation of other DOFs. The required force and torques are application
dependent but these represent typical values that could manipulate a magnetic
robot while overcoming external forces such as gravity or tether drag [25],
[26].

(a) Activation of individual DOF with IPMs parallel with a separation of 15 cm
($\theta$ = $90^{\circ}$, $d$ = $0.15$ m) with three EPMS.

(b) Activation of individual DOF with IPMs parallel with a separation of 5 cm
($\theta$ = $0^{\circ}$, $d$ = $0.15$ m) with three EPMs.

(c) Mean crosstalk across all independent activations of DOF for parallel IPMs
compared between dual and three EPMs.

Figure 4:

###  III-A Two EPMs

We first consider the manipulation of two IPMs using two EPMs. Here, the
number of control variables (i.e. the DOF of EPM poses, 10) matches the number
of output variables. As stated above, with significant separation of the IPMs,
their control can be considered separately with an individual EPM manipulating
a single IPM. This scenario would lead to full control. Once this separation
is decreased, then this independence decreases. For independence of control,
the best case scenario for relative orientation of IPMs is if they are
orthogonal ($\theta$ = 90∘). From our previous works, we have shown this
provides a minimum of eight DOF with zero separation [21]. If this condition
was applied for the case of parallel IPMs (where the resultant magnetic field
and gradient vector is the same across them), this would only lead to 5
independent DOFs with the actuation between IPMs fully coupled. Instead, we
make no assumption of constant gradient leading to this greater DOF.

To observe this independence, we utilize a metric to define the ability to
independently actuate a single DOF. We define crosstalk as

$C(\mathbf{U},i)=\max_{j\in\\{1,10\\}}\frac{|U_{j}|}{|U_{i}|},\text{ }i\neq
j.$ (14)

where $i$ is the index of the desired independently actuated component of
$\mathbf{U}$. A crosstalk of 0% would show that the configuration of the EPMs
can individually activate $U_{i}$ without any activation of any other
component.

To evaluate the ability to independently actuate all ten DOF of the two IPMs
using two EPMs, we utilize our presented optimization strategy to find, if
given the constraints on EPM poses, each DOF that can be activated with
minimal activation of the others. Fig. 3(a) and 3(b) show the relative
activation of each DOF for the case of orthogonal and parallel IPMs with a
separation of 15 cm ($d$ = $0.15$ m) and result in small mean crosstalks of
3.9% and 15.0% respectively. Fig. 3(d) and 3(e) show the same case of
orthogonal and parallel IPMs but with a reduced separation of 5 cm ($d$ =
$0.05$ m). As expected, the independence of control is lost with a smaller
distance between IPMs. In particular, it can be observed in the parallel case,
tight coupling with the torques of the IPMs with crosstalk reaching 100%. This
suggests that two EPMs cannot create largely different directions of magnetic
fields over short distances. The forces however, (which are proportional to
the field gradients), have much less coupling. Fig. 3(c) shows how with
greater IPM separation, the mean crosstalk decreases. Fig. 3(f) shows how the
relative IPM angle affects the maximum crosstalk when independently activating
each DOF. It can be observed, that the greater angle leads to greater
independence of activation.

###  III-B Three EPMs

The use of more but smaller EPMs comes with several advantages for independent
control. The first is the ability to create more complex magnetic fields
because of the greater number of sources of magnetic fields. Secondly, the
mutual attraction between EPMs is decreased allowing them to be operating in
closer proximity safely. However, with more EPMs the number of control
variables increases. This complicates the control, but may come with the
advantage of redundancy due to multiple solutions. Two EPMs are enough to
completely control the 8 magnetic DOF of a singular point in the workspace,
but with greater separation, the use of three EPMs allows us to control more
DOF across multiple points in the workspace. Fig. 4(a) shows the case of
parallel IPMs with a separation of 15 cm ($\theta$ = 0∘ $d$ = $0.15$ m) for
the three EPM case achieving a crosstalk of just 1.6% compared to 15.0% of the
dual EPM case. However, at small separations, the same issue of torque
coupling is apparent (see Fig. 4(b)). Fig. 4(c) compares the separation with
crosstalk for both the dual and three EPM cases.

##  IV Dynamic Simulation

Figure 5: Overview of our simulation setup. Two magnetic robots each embedded
with a singular IPM are controlled through pre-planned trajectories within the
workspace. The controller positions the EPMs in order to give the desired
wrench on the robots.

(a) Robot trajectories with a minimal separation of 60mm. The controller
maintains stability throughout the entire trajectory. Red crosses correspond
to robot 1’s path. Blue crosses correspond to robot 2’s path.

(b) Robot trajectories with a minimal separation of 30mm. As the robots enter
close proximity, the controller becomes instable and the robots diverge from
the path. Yellow crosses correspond to robot 1’s path. Purple crosses
correspond to robot 2’s path.

(c) Positional error for both robots in both scenarios. In the second
scenario, at $t$ = $51$ s, control is lost when the robots enter close
proximity ($33$ mm).

Figure 6:

In order to validate the ability to control two magnetic robots in close
proximity with multiple EPMs, we utilize a dynamic planar simulation of two
robots each embedded with a singular IPM. By manipulating the pose of the EPMs
in the simulation, the magnetic wrench on IPMs can be controlled. The IPMs
embedded in the robots are the same cylindrical NdFeB magnets used in our
above analysis ($5$ mm diameter, $5$ mm length). Although it was observed that
three EPMs lead to better independence of control, we choose to focus this
simulation analysis on the case of two EPMs due to the simpler control
scenario and being more applicable to our previous work based on real-world
dual EPM actuation [18]. As before, the EPMs are cylindrical NdFeB magnets
($100$ mm diameter, $100$ mm length) and the robots operate in a spherical
workspace with radius $r_{min}$ = $15$ cm.

Pre-planned trajectories for the robots are formulated within the central
workspace. These trajectories are formed of sets of target poses for the robot
at a given time-step. A PID controller calculates the required wrench on each
robot in order to reduce the error between the current pose and target pose at
a frequency of 2Hz. In our simulation, we presume the EPMs can be moved
instantaneously between poses. This simplification would not be applicable to
real world application, where the transition of EPMs will affect the wrench on
the robots, but for the purpose of this feasibility study suffices for
understanding if control of multiple magnetic robots is possible with EPMs.
Fig. 5 visualizes the robots and the pre-planned trajectories. The inter-agent
wrench between IPMs is modelled in the simulation. The pre-planned
trajectories are chosen in order to have the robots start with a significant
distance between them (in order to ensure their control can be separated) but
come to a minimal separation $s$ (see Fig. 5). Videos showing the robot
trajectories and the corresponding EPM poses can be found in the Supplementary
Video.

To compare the effect of separation on the control of the robots, we alter the
trajectories to control the separation $s$. In the first scenario (Fig. 6(a)),
the minimal separation is $60$ mm. The controller is capable of maintaining
control of the robot poses throughout the trajectory with a mean positional
error of 8.7mm and mean angular error of $16.7^{\circ}$. Although from our
above analysis, at these small separations, there is large crosstalk between
applied wrench, the smooth trajectories used here require small wrenches on
the robot to follow and therefore control is maintained.

The second scenario (Fig. 6(b)) shows a reduced minimal separation of $30$ mm.
At $t$ = $51$ s, the robots enter into close proximity ($33$ mm), and the
controller fails (Fig. 6(c)). This is due, to the inability of the optimizer
to find suitable EPM poses that can generate the required wrench on the
robots. At this short distance the inter-agent wrench is large, and the
inability to create large change in magnetic field properties over the short
distance means that the wrench between robots cannot be sufficiently
controlled independently. In reality, this would lead to a instability where
the robots would attract each other and collide.

##  V Conclusion

In this work, we have considered the ability to control two magnetic robots
within a shared workspace with multiple manipulated EPMs. We have developed an
optimization scheme capable of generating the optimum pose for EPMs in order
to induce a desired magnetic wrench onto the robots. It can be observed that
with separation between robots ($15$ cm), full independence of control ($3.9$%
crosstalk) can be achieved; however, with smaller distances and alignment, the
coupling between the torques on both robots increases. We have analyzed how,
the relative orientation and separation affects this coupling and how it could
be potentially countered by using more EPMs. Finally, we have utilized our
optimization strategy as part of a dynamic controller in order to control two
simulated magnetic robots in varying proximity on pre-planned trajectories.
Our simulation confirms our analysis, showing that at small separations the
coupling between robots leads to instability in our control. Future work will
consider how these robots can be utilised to operate together to perform
collaborative operation or operate separately in various parts of the
gastrointestinal tract.

In our control scenario, we have assumed that EPMs can be moved instantly
between time steps. In reality, the transition between poses will affect the
resultant wrench on IPMs and must be considered. The stability of our scenario
has not been considered and a better understanding of how accumulating errors
will affect the control must be considered for real world feasibility. Future
work, would consider how collaborative manipulators could be utilized for
controlling the EPMs. This would impose tighter and more complex constraints
on potential EPM poses due to the reach of the manipulators. Currently, there
is no constraint between sequential EPM poses, this should be constrained in
order to produce realistic EPM trajectories, this may limit controllability.
The tight coupling between magnetic torques at small distances limits the
controllability with small separations. Differences in geometric or magnetic
design between magnetic robots could be utilized to better separate their
response to similar magnetic fields. The ability to control the robots is
dependent on receiving good feedback of their relative position to the EPMs;
future work should consider how this feedback can be obtained using magnetic
field or imaging sensors.

## References

  * [1] Y. Kim and X. Zhao, “Magnetic soft materials and robots,” _Chemical Reviews_ , vol. 122, no. 5, pp. 5317–5364, 2022.
  * [2] J. W. Martin, B. Scaglioni, J. C. Norton, V. Subramanian, A. Arezzo, K. L. Obstein, and P. Valdastri, “Enabling the future of colonoscopy with intelligent and autonomous magnetic manipulation,” _Nature machine intelligence_ , vol. 2, no. 10, pp. 595–606, 2020.
  * [3] E. Diller and M. Sitti, “Micro-scale mobile robotics,” _Foundations and Trends® in Robotics_ , vol. 2, no. 3, pp. 143–259, 2013. [Online]. Available: http://dx.doi.org/10.1561/2300000023
  * [4] M. Van Damme, P. Beyl, B. Vanderborght, R. Versluys, R. Van Ham, I. Vanderniepen, F. Daerden, and D. Lefeber, “The safety of a robot actuated by pneumatic muscles—a case study,” _International Journal of Social Robotics_ , vol. 2, no. 3, pp. 289–303, 2010.
  * [5] C. Forbrigger, A. Lim, O. Onaizah, S. Salmanipour, T. Looi, J. Drake, and E. D. Diller, “Cable-less, magnetically driven forceps for minimally invasive surgery,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 2, pp. 1202–1207, 2019.
  * [6] D. Son, H. Gilbert, and M. Sitti, “Magnetically actuated soft capsule endoscope for fine-needle biopsy,” _Soft robotics_ , vol. 7, no. 1, pp. 10–21, 2020.
  * [7] J. Nam, Y. P. Lai, L. Gauthier, G. Jang, and E. Diller, “Resonance-based design of wireless magnetic capsule for effective sampling of microbiome in gastrointestinal tract,” _Sensors and Actuators A: Physical_ , p. 113654, 2022.
  * [8] P. Shokrollahi, Y. P. Lai, S. Rash-Ahmadi, V. Stewart, M. Mohammadigheisar, L.-A. Huber, N. Matsuura, A. E. Zavodni, J. Parkinson, and E. Diller, “Blindly controlled magnetically actuated capsule for noninvasive sampling of the gastrointestinal microbiome,” _IEEE/ASME Transactions on Mechatronics_ , vol. 26, no. 5, pp. 2616–2628, 2020.
  * [9] W. Chen, X. Chen, M. Yang, S. Li, X. Fan, H. Zhang, and H. Xie, “Triple-configurational magnetic robot for targeted drug delivery and sustained release,” _ACS Applied Materials & Interfaces_, vol. 13, no. 38, pp. 45 315–45 324, 2021.
  * [10] E. Diller, J. Giltinan, and M. Sitti, “Independent control of multiple magnetic microrobots in three dimensions,” _The International Journal of Robotics Research_ , vol. 32, no. 5, pp. 614–631, 2013.
  * [11] D. Wong, E. B. Steager, and V. Kumar, “Independent control of identical magnetic robots in a plane,” _IEEE Robotics and Automation Letters_ , vol. 1, no. 1, pp. 554–561, 2016.
  * [12] S. Salmanipour and E. Diller, “Eight-degrees-of-freedom remote actuation of small magnetic mechanisms,” in _2018 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2018, pp. 3608–3613.
  * [13] S. Tottori, L. Zhang, K. E. Peyer, and B. J. Nelson, “Assembly, disassembly, and anomalous propulsion of microscopic helices,” _Nano letters_ , vol. 13, no. 9, pp. 4263–4268, 2013.
  * [14] M. Salehizadeh and E. Diller, “Two-agent formation control of magnetic microrobots in two dimensions,” _Journal of Micro-Bio Robotics_ , vol. 12, pp. 9–19, 2017.
  * [15] ——, “Three-dimensional independent control of multiple magnetic microrobots via inter-agent forces,” _The International Journal of Robotics Research_ , vol. 39, no. 12, pp. 1377–1396, 2020.
  * [16] F. Ongaro, S. Pane, S. Scheggi, and S. Misra, “Design of an electromagnetic setup for independent three-dimensional control of pairs of identical and nonidentical microrobots,” _IEEE transactions on robotics_ , vol. 35, no. 1, pp. 174–183, 2018.
  * [17] Q. Boehler, S. Gervasoni, S. L. Charreyron, C. Chautems, and B. J. Nelson, “On the workspace of electromagnetic navigation systems,” _IEEE Transactions on Robotics_ , pp. 1–17, 2022.
  * [18] G. Pittiglio, M. Brockdorff, T. da Veiga, J. Davy, J. H. Chandler, and P. Valdastri, “Collaborative magnetic manipulation via two robotically actuated permanent magnets,” _IEEE Transactions on Robotics_ , 2022.
  * [19] J. J. Abbott, E. Diller, and A. J. Petruska, “Magnetic methods in robotics,” _Annual Review of Control, Robotics, and Autonomous Systems_ , vol. 3, no. 1, pp. 57–90, 2020.
  * [20] M. Yousefi and H. Nejat Pishkenari, “Independent position control of two identical magnetic microrobots in a plane using rotating permanent magnets,” _Journal of Micro-Bio Robotics_ , vol. 17, no. 1, pp. 59–67, 2021.
  * [21] G. Pittiglio, J. H. Chandler, M. Richter, V. K. Venkiteswaran, S. Misra, and P. Valdastri, “Dual-arm control for enhanced magnetic manipulation,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2020, pp. 7211–7218.
  * [22] G. Pittiglio, P. Lloyd, T. da Veiga, O. Onaizah, C. Pompili, J. H. Chandler, and P. Valdastri, “Patient-specific magnetic catheters for atraumatic autonomous endoscopy,” _Soft Robotics_ , 2022.
  * [23] S. Salmanipour, O. Youssefi, and E. D. Diller, “Design of multi-degrees-of-freedom microrobots driven by homogeneous quasi-static magnetic fields,” _IEEE Transactions on Robotics_ , vol. 37, no. 1, pp. 246–256, 2020.
  * [24] A. J. Petruska and J. J. Abbott, “Optimal permanent-magnet geometries for dipole field approximation,” _IEEE transactions on magnetics_ , vol. 49, no. 2, pp. 811–819, 2012.
  * [25] B. Scaglioni, L. Previtera, J. Martin, J. Norton, K. L. Obstein, and P. Valdastri, “Explicit model predictive control of a magnetic flexible endoscope,” _IEEE robotics and automation letters_ , vol. 4, no. 2, pp. 716–723, 2019.
  * [26] G. Pittiglio, L. Barducci, J. W. Martin, J. C. Norton, C. A. Avizzano, K. L. Obstein, and P. Valdastri, “Magnetic levitation for soft-tethered capsule colonoscopy actuated with a single permanent magnet: A dynamic control approach,” _IEEE robotics and automation letters_ , vol. 4, no. 2, pp. 1224–1231, 2019.