[a]Gojko Vujanovic

# Multiscale evolution of heavy flavor in the QGP

###### Abstract

Shower development dynamics for a jet traveling through the quark-gluon plasma
(QGP) is a multiscale process, where the heavy flavor mass is an important
scale. During the high virtuality portion of the jet evolution in the QGP,
emission of gluons from a heavy flavor is modified owing to heavy quark mass.
Medium-induced radiation of heavy flavor is sensitive to microscopic processes
(e.g. diffusion), whose virtuality dependence is phenomenologically explored
in this study. In the lower virtuality part of shower evolution, i.e. when the
mass is comparable to the virtuality of the parton, scattering and radiation
processes of heavy quarks differ from light quarks. The effects of these
mechanisms on shower development in heavy flavor tagged showers in the QGP is
explored here. Furthermore, this multiscale study examines dynamical pair
production of heavy flavor (via virtual gluon splittings) and their subsequent
evolution in the QGP, which is not possible otherwise. A realistic event-by-
event simulation is performed using the JETSCAPE framework. Energy-momentum
exchange with the medium proceeds using a weak coupling recoil approach. Using
leading hadron and open heavy flavor observables, differences in heavy versus
light quark energy-loss mechanisms are explored, while the importance of heavy
flavor pair production is highlighted along with future directions to study.

##  1 Introduction

JETSCAPE (Jet Energy-loss Tomography with a Statistically and Computationally
Advanced Program Envelope) framework is a flexible publicly released event-
generator framework providing theoretical modeling of all aspects of heavy ion
collisions. While in this contribution focus is given towards the quenching of
high energy partons, modules are present to cover other aspects of aspect
relativistic heavy-ions, notably relativistic dissipative hydrodynamics and
hadronic transport. Each module can be user-modified, giving any user of the
JETSCAPE framework the ability to focus on a subset of physics simulations of
interest, knowing that other aspects of the simulation are state of the art.
Indeed, the JETSCAPE framework also supplies Bayesian analysis tools needed to
accurately constrain physical processes. These Bayesian tools are actively
being used [1, 2, 3, 4] to improve our understanding of relativistic heavy-ion
physics. The work herein investigates how different physical mechanisms
involved in jet quenching affect the nuclear modification factor $R_{AA}$.
Specifically, this contribution explores (i) how dynamical generation of heavy
quarks inside the QGP effects light and heavy flavor $R_{AA}$ and (ii)
investigates the effects of a virtuality-dependent $\hat{q}$ on light and
heavy flavor $R_{AA}$.

##  2 Numerical simulations

The simulation of parton interaction with the QGP uses a two-step approach for
both light and heavy flavor partons: the high virtuality (and high energy)
portion of parton quenching in the QGP is described using the higher twist
formalism [5, 6] implemented in MATTER (Modular All Twist Transverse-
scattering Elastic-drag and Radiation) [7, 8], while the low virtuality (and
high energy) portion of quenched partons is modeled via Linear Boltzmann
Transport (LBT) [9].

MATTER contains both vacuum showering of highly virtual partons as well as in-
medium corrections to the splitting function, which includes novel coherence
effects giving a virtuality-dependent $\hat{q}$ first proposed in Ref. [10].
The medium-modified splitting function of the process $Q\to Q+g$ — which also
holds for light flavor in the limit $M\to 0$ — is [7]:

$\displaystyle\frac{dN^{\rm vac}}{dzdt}+\frac{dN^{\rm med}}{dzdt}$
$\displaystyle=$ $\displaystyle\frac{\alpha_{s}(t)}{2\pi}\frac{P_{g\leftarrow
Q}(z)}{t}\left\\{1+\int^{\tau^{+}_{Q}}_{0}d\tau^{+}\frac{1}{z(1-z)t(1+\chi)^{2}}\left[2-2\cos\left(\frac{\tau^{+}}{\tau^{+}_{Q}}\right)\right]\right.\times$
$\displaystyle\times$
$\displaystyle\left.\left[\left(\frac{1+z}{2}\right)-\chi+\left(\frac{1+z}{2}\right)\chi^{2}\right]\hat{q}(t)\right\\}.$

where $z$ is the momentum fraction of the daughter heavy quark, $M$ is the
mass of the heavy quark, $\chi=(1-z)^{2}M^{2}/l^{2}_{\perp}$, with
$l^{2}_{\perp}$ being the relative transverse momentum square between the
outgoing daughter partons, determined via $z(1-z)t=l^{2}_{\perp}(1+\chi)$,
while $t$ is the virtuality of the heavy quark and $P_{g\leftarrow
Q}(z)=C_{F}\left(\frac{1+z^{2}}{1-z^{\,\,\,}}\right)$ is the splitting
function. The integral over light-cone time $\tau^{+}$ in Eq. (LABEL:eq:hq_ht)
assumes the medium is in its rest frame, with the upper bound
$\tau^{+}_{Q}=2q^{+}/t$ being given by the ratio of forward light-cone
momentum $q^{+}=\left(q^{0}+{\bf q}\cdot\hat{n}\right)/\sqrt{2}$ (with
$\hat{n}={\bf q}/|{\bf q}|$), and the virtuality $t$. The virtuality-dependent
$\hat{q}(t)$ is parametrized as

$\displaystyle\hat{q}(t)$ $\displaystyle=$
$\displaystyle\hat{q}^{HTL}\frac{c_{0}}{1+c_{1}\ln^{2}(t)+c_{2}\ln^{4}(t)}$

where $c_{1}$ and $c_{2}$ are tunable parameters, $t=E^{2}-{\bf p}^{2}-m^{2}$
is the virtuality of the parton with $p^{\mu}=(E,{\bf p})$, and
$c_{0}=1+c_{1}\ln^{2}(t_{s})+c_{2}\ln^{4}(t_{s})$ is an overall normalization
ensuring that the $t$-dependent contribution is unitless and lies within 0 and
1 as $t$ does not go below $t_{s}$, while the expression for $\hat{q}^{HTL}$
is in Ref. [11].

In addition to medium-modified radiation, MATTER includes $2\to 2$ scatterings
at leading order in the QCD coupling $\alpha_{s}$. An important process that
can only be studied in a multistage model, such as MATTER+LBT, is the
dynamical in-medium heavy flavor production. Thus, heavy flavor pair
production via $g\to Q+\bar{Q}$, where $(\bar{Q})Q$ is the heavy (anti)quark,
is the key process studied here. Indeed, a multiscale study, such as
MATTER+LBT, is required to fully estimate the importance of dynamical in-
medium heavy quark production. Of course, charm pair production is the process
with the largest allowed phase space, and hence D-meson $R_{AA}$, in addition
to charged hadron $R_{AA}$, will be explored below. Unlike recent calculations
[12, 13], the present simulation relies on light-flavor $\hat{q}$ effects when
correcting the $P_{Q\leftarrow g}$ splitting function. However, the medium-
modified $P_{Q\leftarrow g}$ splitting function is taking into account
kinematic cuts associated with heavy quark masses [11]. This approximation for
$P_{Q\leftarrow g}$ will be revisited in the future, as discussed below. In
the low virtuality sector, the LBT simulation is used which also includes
$2\to 2$ scatterings, and employs the radiation kernel of MATTER within a
$2\to 3$ process to describe in-medium emissions (see Ref. [11] for details).
As far as simulations of the QGP and its subsequent hadronic evolution is
concerned, the Bayesian maximum a posteriori parameters are used [14] to
generate the QGP medium for Pb-Pb collisions at $\sqrt{s_{NN}}=5.02$ TeV in
the 0-10% centrality class (cf. details in [11]).

The observable of interest herein is $R_{AA}$:

$\displaystyle
R_{AA}=\frac{\frac{d\sigma_{AA}}{dp_{T}}}{\frac{d\sigma_{pp}}{dp_{T}}}=\frac{\sum_{\ell}\frac{dN_{AA,\ell}}{dp_{T}}\hat{\sigma}_{\ell}}{\sum_{\ell}\frac{dN_{pp,\ell}}{dp_{T}}\hat{\sigma}_{\ell}}$
(3)

where $\frac{dN_{AA}}{dp_{T}}$ and $\frac{dN_{pp}}{dp_{T}}$ are the
multiplicity of either D-mesons or charged hadrons originating from A-A and
p-p collisions in the experimentally given $p_{T}$ bin, respectively. The
cross-section for producing the hard scattering process of the given range
$\ell$ in transverse momentum $\hat{p}_{T}$ is $\hat{\sigma}_{\ell}$
($\hat{p}_{T}$ is the transverse momentum of the exchanged parton in the hard
scattering sampled by PYTHIA). Many $\hat{\sigma}_{\ell}$ are sampled to span
the large kinematic range of the collision [11].

##  3 Results

|
---|---
Figure 1: Nuclear modification factor for D-mesons (a) and charged hadrons (b)
at the $\sqrt{s_{NN}}=5.02$ TeV Pb-Pb collisions at the LHC in the 0-10%
centrality. $c_{1}=10,\ c_{2}=100$ are set within the $\hat{q}(t)$
parametrization [see in Eq. (LABEL:eq:qhat_t)] for the MATTER alone and the
MATTER+LBT curve. The other parameter is the switching virtuality $t_{s}=4$
GeV2 between MATTER and LBT, whose value is found to best describe the
$R_{AA}$ data. Finally, the p-p baseline for the LBT curve is calculated using
PYTHIA whereas the p-p baseline for the MATTER and MATTER+LBT cases are
calculated using MATTER vacuum [15]. Data taken from Ref. [16, 17, 18].

Figure 1 shows the contribution for MATTER and LBT, as well as the combined
result, compared to D-meson and charged hadrons $R_{AA}$, respectively. While
it may be possible to change $\hat{q}$ in the LBT-only portion of the
simulation to obtain a good description of D-meson $R_{AA}$, such a change
would not allow for simultaneous description of both D-meson and charged
hadron $R_{AA}$ without a virtuality-dependent $\hat{q}(t)$. The sensitivity
to $c_{1}$ and $c_{2}$, was well as the importance of having a virtuality-
dependent $\hat{q}$ is shown in Fig. 2, where a similar sensitivity to
$\hat{q}(t)$ is observed between light and heavy flavor. Beyond the
realization that a non-trivial $\hat{q}(t)$ is required, further constraining
$(c_{1},c_{2})$ is difficult and will be explored via Bayesian analysis [4].

|
---|---
Figure 2: (Color online) Nuclear modification factor for D-mesons (a) and
charged hadrons (b) in $\sqrt{s_{NN}}=5.02$ TeV Pb-Pb collisions at the LHC at
0-10% centrality. The parametrization of $\hat{q}(t)$ is varies with $c_{1}$
and $c_{2}$ values indicated. The ratio in the bottom plots are taken with
respect to the $c_{1}=10,c_{2}=100$ case with $\hat{q}(t)$ parametrization
[see Eq. (LABEL:eq:qhat_t)]. The sensitivity of $R_{AA}$ to $t_{s}$ is
explored in Ref. [11].

Figure 3 shows the importance of in-medium $g\to Q+\bar{Q}$ to both D-meson
and charged hadron $R_{AA}$. According to the ratio plot in Fig. 3, the
D-meson $R_{AA}$ is significantly affected by $g\to Q+\bar{Q}$, on the order
of 20%. Therefore, any future phenomenological studies of heavy flavor energy
loss cannot neglect dynamical $g\to Q+\bar{Q}$ production. The next step for
this calculation is to compute the medium-modification to $P_{Q\leftarrow g}$
splitting function, as well as devise a $\hat{q}$ that has both a $t$\- and
$M$-scale dependence, i.e. $\hat{q}(t,M)$. The formalism to calculate the
medium-correction to the $P_{g\leftarrow Q}$ splitting function is already
established using Soft Collinear Effective Theory (SCET) [7]; that approach
needs to be adapted for $P_{Q\leftarrow g}$. Furthermore, to obtain
$\hat{q}(M,t)$ a modification of the work in Ref. [10] would be needed.

|
---|---
Figure 3: Nuclear modification factor for D-mesons (a) and charged hadrons (b)
in $\sqrt{s_{NN}}=5.02$ TeV Pb-Pb collisions at the LHC at 0-10% centrality.
$c_{1}=10,\ c_{2}=100$ parameters values are employed in Eq.
(LABEL:eq:qhat_t). Ignoring the $g\rightarrow Q+\bar{Q}$ process in MATTER
impacts the D meson $R_{AA}$, while it has a smaller effect on the charged
hadron $R_{AA}$. The dashed line in the ratio subplots divides MATTER+LBT no
$g\to Q+\bar{Q}$ to MATTER+LBT, while the dotted-dashed line divides
MATTER+LBT HTL no $g\to Q+\bar{Q}$ to MATTER+LBT HTL.

##  4 Conclusion

A realistic calculation of heavy quark interaction within the QGP has been
presented, containing both primordial heavy quarks as well as those generated
during the shower of a highly virtual parton inside hot and dense nuclear
matter. The latter can only be fully described using a multiscale numerical
simulation, such as MATTER+LBT. The dynamical production of heavy flavor
during the showering of hard partons, i.e. via the $g\to Q+\bar{Q}$ process,
contributes significantly (at the 20% level) to the overall multiplicity of
heavy quarks, and thus cannot be ignored. Future directions entail improving
the medium-induced $P_{Q\leftarrow g}$ splitting function, and the multiscale
dependence of the transport coefficient $\hat{q}(t,M)$, using SCET formalism
developed in Refs. [7, 10] to include both the virtuality and heavy flavor
mass scales.

The calculation presented herein also shows that $\hat{q}(t)$ is needed to
describe both light flavor $R_{AA}$ as well as D-meson $R_{AA}$. However to
fully constrain the virtuality-dependence of $\hat{q}$, the improvements
mentioned above — i.e. to the $g\to Q+\bar{Q}$ splitting function and
$\hat{q}(t,M)$ — will be needed. Once at hand, they can be applied within a
systematic Bayesian approach to constrain $\hat{q}$. The JETSCAPE
Collaboration has the necessary tools and know-how to carry out such an
endeavor in the future.

## References

  * [1] D. Everett et al.,  Phys. Rev. C 103, 054904 (2021), 2011.01430.
  * [2] D. Everett et al.,  Phys. Rev. Lett. 126, 242301 (2021), 2010.03928.
  * [3] S. Cao et al.,  Phys. Rev. C 104, 024905 (2021), 2102.11337.
  * [4] Y. Chen et al.,  these proceedings.
  * [5] X.-N. Wang and X.-f. Guo,  Nucl. Phys. A 696, 788 (2001), hep-ph/0102230.
  * [6] A. Majumder,  Phys. Rev. D 85, 014023 (2012), 0912.2987.
  * [7] R. Abir and A. Majumder,  Phys. Rev. C 94, 054902 (2016), 1506.08648.
  * [8] A. Majumder,  Phys. Rev. C 88, 014909 (2013), 1301.5323.
  * [9] T. Luo, S. Cao, Y. He, and X.-N. Wang,  Phys. Lett. B 782, 707 (2018), 1803.06785.
  * [10] A. Kumar, A. Majumder, and C. Shen,  Phys. Rev. C 101, 034908 (2020), 1909.03178.
  * [11] W. Fan et al.,  Phys. Rev. C 107, 054901 (2023), 2208.00983.
  * [12] M. Attems et al.,  (2022), 2209.13600.
  * [13] M. Attems et al.,  JHEP 01, 080 (2023), 2203.11241.
  * [14] J. E. Bernhard, J. S. Moreland, and S. A. Bass,  Nature Phys. 15, 1113 (2019).
  * [15] A. Kumar et al.,  Phys. Rev. C 102, 054906 (2020), 1910.05481.
  * [16] A. M. Sirunyan et al.,  Phys. Lett. B 782, 474 (2018), 1708.04962.
  * [17] V. Khachatryan et al.,  JHEP 04, 039 (2017), 1611.01664.
  * [18] S. Acharya et al.,  JHEP 10, 174 (2018), 1804.09083.


# Design and evaluation of a multi-finger skin-stretch tactile interface for
hand rehabilitation robots

Alexandre L. Ratschat1,2, Rubén Martín-Rodríguez1, Yasemin Vardar1, Gerard M.
Ribbers2,
and Laura Marchal-Crespo1,2,3 *This work was funded by the Dutch Research
Council (NWO, VIDI Grant Nr. 18934)1A. L. Ratschat, R. Martín-Rodríguez, Y.
Vardar, and L. Marchal-Crespo are with the Department of Cognitive Robotics,
Delft University of Technology, Delft, The Netherlands2A. L. Ratschat, G. M.
Ribbers, and L. Marchal-Crespo are with the Department of Rehabilitation
Medicine, Erasmus MC, University Medical Center Rotterdam, Rotterdam, The
Netherlands3L. Marchal-Crespo is with the ARTORG Center for Biomedical
Engineering Research, University of Bern, Switzerland Corresponding author:
<EMAIL_ADDRESS>
###### Abstract

Object properties perceived through the tactile sense, such as weight,
friction, and slip, greatly influence motor control during manipulation tasks.
However, the provision of tactile information during robotic training in
neurorehabilitation has not been well explored. Therefore, we designed and
evaluated a tactile interface based on a two-degrees-of-freedom moving
platform mounted on a hand rehabilitation robot that provides skin stretch at
four fingertips, from the index through the little finger. To accurately
control the rendered forces, we included a custom magnetic-based force sensor
to control the tactile interface in a closed loop. The technical evaluation
showed that our custom force sensor achieved measurable shear forces of $\pm$
8 N with accuracies of 95.2–98.4% influenced by hysteresis, viscoelastic
creep, and torsional deformation. The tactile interface accurately rendered
forces with a step response steady-state accuracy of 97.5–99.4% and a
frequency response in the range of most activities of daily living. Our sensor
showed the highest measurement-range-to-size ratio and comparable accuracy to
sensors of its kind. These characteristics enabled the closed-loop force
control of the tactile interface for precise rendering of multi-finger two-
dimensional skin stretch. The proposed system is a first step towards more
realistic and rich haptic feedback during robotic sensorimotor rehabilitation,
potentially improving therapy outcomes.

##  I Introduction

Robotic systems, in combination with virtual reality (VR), can be used in
high-intensity movement training in motivating and safe environments [1, 2],
incorporating games that resemble activities of daily living (ADLs), such as
object manipulation [3]. Yet, the ecological validity of these games is
limited, which might explain the limited gains in ADLs after robotic practice
[4]. A lack of tactile information—e.g., object friction, weight, and
slip—during robotic training may contribute to this lack of ecological
validity [5], and thus the transfer of gained skills to ADLs [3]. The
importance of tactile input in motor performance is illustrated in clinical
syndromes, such as tactile apraxia [6] and sensory ataxia [7].

Tangential skin deformations, commonly called skin stretch, especially arise
during object interactions and could contribute to robotic training. Applying
skin stretch can be exploited to render object friction, as demonstrated by
Provancher and Silvester using a single-degree-of-freedom (DoF) tactor device
on a single finger [8]. In addition to friction, Schorr et al. showed that
skin stretch could render virtual object weight and stiffness with two
wearable three-DoF tactor devices on the user’s index finger and thumb [9].
Quek et al. provided interaction force and torque information via a three-
finger pinch grasp by combining a three-DoF tactor-based tactile interface
with a commercial three-DoF grounded kinesthetic feedback device [10].

A main limitation of most skin stretch tactile interfaces is that they
typically render forces by applying a skin displacement computed through an
estimate of the finger pad stiffness (e.g., [11, 10]). However, considering
the large variability in finger pad stiffness under varying pressure levels
and across individuals [12], the rendered force accuracy is limited. This is
especially relevant since high-fidelity force rendering is required to promote
the transfer of the gained skill in the virtual environment to the real world
[3]. Multi-axis force sensors can overcome this limitation as they enable
closed-loop force control for accurate lateral force rendering at the
fingerpad [13, 14]. However, their high cost and bulky size limit their use in
skin stretch interfaces.

In this work, we present the collected requirements, design, and evaluation of
a novel multi-finger tactile interface that is intended to simulate two-
dimensional skin stretch while training ADLs in virtual environments. The
interface is designed to be integrated into the PRIDE hand rehabilitation
robot developed by Rätz et al. [15], which enables hand setup in a closed
position for patients with hand spasticity. We also describe the design and
evaluation of a custom low-cost and compact magnetic-based force sensor to
render the skin stretch stimuli in a closed-loop fashion. We evaluated our
solution by characterizing the sensor’s static and dynamic errors, hysteresis,
viscoelastic creep, and torsional sensitivity, as well as the tactile
interface’s step and frequency responses. Finally, we discuss our results in
the context of relevant literature.

##  II Methods

###  II-A Tactile interface requirements

We designed our tactile interface with the following requirements in mind: i)
integration in PRIDE (Fig. 1a) while minimally impeding PRIDE’s range of
motion (ROM) and hand donning; ii) two-DoF fingerpad deformations in the
proximodistal and mediolateral directions of the fingers (from index to little
finger) with maximum deformations of 5 $\mathrm{m}\mathrm{m}$ [16]; iii)
maximum renderable forces of 5 $\mathrm{N}$; iv) maximum force rendering error
of 50 $\mathrm{m}\mathrm{N}$ for forces below 1.5 $\mathrm{N}$ and 100
$\mathrm{m}\mathrm{N}$ for larger values based on the just-noticeable
difference (JND) for force perception [17]; v) rise time of 0.35 $\mathrm{s}$
to render object interaction onset [18]; vi) bandwidth of 10
$\mathrm{H}\mathrm{z}$, the highest movement frequency reported during ADLs
[19]; vii) total material costs under 200 € to keep the technology low-cost
and accessible to all.

###  II-B Two-dimensional skin stretch moving platform

Figure 1: Overview of the tactile interface and custom sensor. (a) The tactile
interface is integrated into the PRIDE rehabilitation robot, with the
fingerpads touching the moving platform and pressure applied on the fingers’
dorsal side with the device’s quick release. (b) Custom force sensor with the
magnet embedded in the silicone. (c) Schematic showing the sensor elastomer
deforming under a force F. (d) Front view of the tactile interface with the
moving platform. (e) Back view of the tactile interface with linear rods, DC
motors, and force sensor.

The skin stretch tactile interface was designed to replace PRIDE’s existing
finger end-effector (Fig. 1a). We took advantage of the current quick-release
finger fixation at the back of the fingers to ensure that the platform
movements efficiently transmit into the deformations of the fingerpads.

The platform (Fig. 1d-e) was actuated with two lead screws (M3x30
$\mathrm{m}\mathrm{m}$ stainless steel threaded rod) in the fingers’
proximodistal (X-axis) and the mediolateral (Y-axis) directions connected to
two brushed DC motors (5:1 Micro Metal Gearmotor HPCB 6V with Extended Motor
Shaft, Pololu, USA). Each axis included 3 $\mathrm{m}\mathrm{m}$ diameter
stainless steel linear rods to stabilize the platform movement. To reduce axle
friction, miniature linear bearings (C-LMU3, MISUMI Europe, Germany) and
Teflon tubing were used for the X and Y axes, respectively. Each axis achieved
a total range of $\pm$4.5 $\mathrm{m}\mathrm{m}$. A set of custom 3D-printed
parts made up the device’s housing and provided connections for the actuators
and the mechanical components.

Magnetic rotary encoders (Magnetic Encoder Pair Kit for Micro Metal
Gearmotors, Pololu, USA) were mounted at the rear shaft of the DC motors to
measure the actuators’ positions. Considering the gearbox and lead screw
reductions, the encoders provided a theoretical resolution of 8.3
$\mu\mathrm{m}$. Each motor was controlled using a DRV8876 brushed DC motor
driver (Texas Instruments, USA) using differential pulse-width modulation
(PWM) in slow-decay mode. An ESP32 Devkit-V1 microcontroller (Espressif,
China) integrated and controlled all the components. A custom-made printed
circuit board was designed to host the drivers, regulator, microcontroller,
and sensor and actuator connectors.

###  II-C Custom-made magnetic-field-based force sensor

To create an accurate closed-loop force rendering system, we decided to place
a force sensor between the platform in contact with the fingerpads and the
actuation mechanism (available volume $\sim$1 $\mathrm{c}\mathrm{m}^{3}$). We
opted for a custom three-DoF magnetic-field-based force sensor due to its
simple and cheap manufacturing process and small size (Fig. 1b-c) [20].

The sensor relied on a magnet embedded in an elastomer on top of a three-axis
magnetometer (MLX90393, Melexis, Belgium) located 3.5 $\mathrm{m}\mathrm{m}$
beneath the magnet (air gap) to increase the sensor’s sensitivity. When forces
are applied, the elastomer deforms and the magnet is displaced, varying the
magnetic field at the magnetometer. The magnetic field was generated by a
3.0$\times$1.0 $\mathrm{m}\mathrm{m}$ cylindrical neodymium magnet with N48
remanence (MAGZ-087-P, MagnetPartner, Denmark), equivalent to $\sim$1.4
$\mathrm{T}$. The strong remanence helped to improve the sensor signal-to-
noise ratio, while its small size facilitated its embedding within the
elastomer. Through a calibration procedure (see section II-E2), the sensor
could measure the applied forces based on the measured variations in the
magnetic field.

The elastomer’s shape (cylindrical), size (11 $\mathrm{m}\mathrm{m}$ diameter,
4.5 $\mathrm{m}\mathrm{m}$ height), and material—a silicone with A30 shore
hardness (DragonSkin 30, Smooth-On, USA)—were chosen to achieve the stiffness
needed to cover the range of desired measurable forces. The sensor was
manufactured by mold-casting the elastomer body with a cavity for embedding
the magnet, which was sealed using Sil-poxy (Smooth-On, USA). The elastomer
was attached to the platform and the magnetometer board using cyanoacrylate
adhesive (SUPERGLUE, Loctite, USA). A grooved pattern on both elastomer
surfaces enhanced the adhesive binding strength.

The material costs of the tactile interface, including two sensors, were below
150 €. Although we embedded only one force sensor in the final design, we
manufactured two for evaluation purposes.

###  II-D Control

The custom sensor’s magnetometer was interfaced via Serial Protocol Interface
(SPI) using a 10 $\mathrm{M}\mathrm{H}\mathrm{z}$ clock signal in single-
measurement mode with a 1.55 $\mathrm{m}\mathrm{s}$ reading time. The sensor
was configured to provide data on the 16 least significant bits of the 19-bit
internal analog-to-digital (ADC) converter along the X and Y axes, providing
the highest resolution possible. The resolution was set one bit lower along
the Z axis, i.e., bits 1–17 of the ADC, due to the strong magnetic field along
that axis.

The ESP32 microcontroller was programmed using FreeRTOS and the control loop
for reading sensor data and sending control actions was set to 500
$\mathrm{H}\mathrm{z}$. The encoder readings and custom force sensor readings
were low-pass filtered in real-time at 10 $\mathrm{H}\mathrm{z}$ using Euler-
backward discretization to avoid high-frequency noise. The interface featured
a voltage-driven proportional-integral-derivative (PID) force controller with
manually tuned gains ($k_{px},k_{py}$) = (13.2,
12)($\mathrm{V}\mathrm{/}\mathrm{N}$), ($k_{ix},k_{iy}$) = (3.6,
12)($\mathrm{V}\mathrm{/}\mathrm{N}\mathrm{s}$), ($k_{dx},k_{dy}$) = (0.36,
0.72)($\mathrm{V}\mathrm{s}\mathrm{/}\mathrm{N}$), for the X and Y axes,
respectively.

###  II-E Experimental evaluation

We first calibrated two custom magnetic-based force sensors: a standalone
sensor and a sensor mounted on the tactile interface. The standalone sensor
was employed to fully evaluate and characterize the sensor type, without space
restrictions from the moving platform. The sensor embedded in the tactile
interface was used to force control and evaluate the system’s step and
frequency response.

####  II-E1 Experimental setup

All measurements and evaluations were conducted with the same experimental
setup (Fig. 2). To measure the forces applied to the sensor and, eventually,
exerted by the tactile interface, an ATI Nano 43 force-torque sensor (ATI
Industrial Automation, USA) was used (reference sensor). A three-axis
positioning stage aligned the reference sensor with the standalone custom
sensor or the tactile interface (depending on the evaluation) to apply
controlled displacements. The stage consisted of a two-axis positioning table
(HBM Machines, The Netherlands) connected through a custom mounting plate to a
RK-10 stand (Dino-Lite Europe, The Netherlands). The custom sensor or tactile
interface was mounted on the positioning table, while the reference sensor was
mounted on the stand, allowing relative displacements between the two. Data
from the reference sensor were collected via a data acquisition card
(USB-6351, National Instruments, USA) and sent to a host computer. The
information from the tactile interface, i.e., timestamp, measured forces,
desired forces, and commanded control action, was sent to the host computer
through the USB port of the ESP32 microcontroller. The management of the
communication and the collection of data were performed in Python.

Figure 2: Experimental setup used to calibrate and evaluate the custom sensors
and tactile interface.

####  II-E2 Sensor calibration

To determine the mapping from magnetic field readings to forces (sensor
calibration), we applied controlled deformations to the sensor using the
positioning stage and recorded paired magnetic fields–applied force data from
the reference sensor, at a sampling rate of 100 $\mathrm{H}\mathrm{z}$.

Training and testing data sets were collected to train and test the mapping.
The training set consisted of 100 samples collected over 1 $\mathrm{s}$ at 165
target locations, yielding 16500 samples in total. The target locations were
generated as relative displacements between the custom and reference sensor
with five equally-distributed steps vertically and radially within [0.0, 0.8]
$\mathrm{m}\mathrm{m}$ and [0.0, -2.0] $\mathrm{m}\mathrm{m}$, respectively,
and eight angular steps within [0, 360] $\mathrm{\SIUnitSymbolDegree}$ around
the vertical axis. The test set contained 45 randomly generated target
locations, yielding 4500 samples in total.

The magnetic field data were mapped to forces by combining a third-order
polynomial expansion with a linear regression—i.e., parameter optimization via
stochastic gradient descent with a Huber loss function and L2 regularization
term, similar to [21]. The evaluation of the custom sensor calibration was
performed by calculating the static error as the L2 norm of the mean absolute
errors (MAE) per axis between the reference and custom sensor from the
training and test datasets. We also report the dynamic error, calculated as
the L2 norm of the MAE between the reference and custom sensor readings while
varying the applied forces by the positioning stage over 30 $\mathrm{s}$, with
forces randomly generated within [-5, 4.5] $\mathrm{N}$, [-8, 3] $\mathrm{N}$
and [-9, 1] $\mathrm{N}$ in the X, Y, and Z axes, respectively.

####  II-E3 Custom force sensor characterization

We characterized the custom standalone sensor by evaluating its hysteresis,
viscoelastic creep, and torsional sensitivity. To evaluate the potential
hysteric error, a loading-unloading cycle was applied with the positioning
stage to the custom sensor through an indentation of 0.5
$\mathrm{m}\mathrm{m}$ in the Z axis and of 1.25 $\mathrm{m}\mathrm{m}$ in the
X and Y axes. During the X and Y measurements, an additional 0.2
$\mathrm{m}\mathrm{m}$ indentation in the Z direction was applied to recreate
the working conditions in which the sensor will be used. The hysteresis value
at each axis was then computed as the difference in the custom sensor readings
at loading and unloading conditions at half of the maximum applied force
(measured by the reference sensor) over the maximum applied force (reported in
percentage).

To evaluate the viscoelastic creep, the same indentations as in the hysteresis
test were applied and sustained for a total of 30 $\mathrm{s}$. The
differences between the readouts of the custom and reference sensors were
computed after 5 $\mathrm{s}$ and 30 $\mathrm{s}$ of the force onset. The
creep was then calculated as the difference between those two values.

To evaluate how the custom sensor readouts were affected by torsional
deformations around each axis (torsional sensitivity), we incorporated custom
connectors between the reference sensor and the positioning stage stand. Each
connector aligned with the custom sensor at a 15
$\mathrm{\SIUnitSymbolDegree}$ offset around each of the three axes, i.e.,
three different connectors were employed. We collected data for 20
$\mathrm{s}$, starting with the custom sensor undeformed and ending fully
deformed around the respective axis. The L2 norm of the MAE between the
reference and standard sensor readouts at rest and fully deformed were then
calculated and compared.

####  II-E4 Force control evaluation

A second force sensor unit was manufactured, mounted within the tactile
interface, and calibrated. We performed a step response test and a frequency
bandwidth test to evaluate the force-control performance. In both tests,
measured forces, desired forces, and the commanded control action were
collected at 500 $\mathrm{H}\mathrm{z}$ from the tactile interface, matching
the control frequency of the tactile interface. A normal force of 2
$\mathrm{N}$ was always applied with the positioning stage.

For the step response test, we recorded the system’s response during 10 steps
with an amplitude of 1, 2, and 3 $\mathrm{N}$ along the positive X, Y, and XY
directions (30 steps per direction). Each step lasted 2 $\mathrm{s}$ and was
preceded and followed by 1 $\mathrm{s}$ of zero force. The steady-state error
was computed as the error after 1.8 $\mathrm{s}$ of the step onset, while the
rise time was computed as the time required to rise from 10% to 90% of the
steady-state value of the custom sensor signal. These metrics were computed
for each step and averaged across the 30 steps per direction.

For the bandwidth test, we commanded an exponential chirp signal that swept
from 1 to 100 $\mathrm{H}\mathrm{z}$ at a constant amplitude of 1 $\mathrm{N}$
(2 $\mathrm{N}$ peak-to-peak) in the X and Y directions, for 10 seconds each.
We applied a Hann window to the commanded forces and the forces measured by
the custom sensor to reduce spectral leakage, followed by a Fast Fourier
Transform (FFT). The magnitudes and phases of the Bode plots illustrating the
system’s frequency responses were calculated as the magnitudes and phases of
the measured force FFTs divided by the commanded force FFTs.

##  III Results

###  III-A Sensor calibration

During the calibration procedure of the custom standalone sensor, we recorded
maximum absolute forces of approximately (8, 8, 16) $\mathrm{N}$ in the X, Y,
and Z axes, respectively. The static error for the training and test datasets
of the standalone sensor, as well as the dynamic error, can be found in Table
I.

TABLE I: Standalone custom sensor static and dynamic error

Evaluation (Dataset) | Error X, Y, Z (MAE $\pm$ SD) / N
---|---
Static Error (Training) |  $0.068\pm 0.102$, $0.056\pm 0.053$, $0.250\pm 0.260$
Static Error (Test) |  $0.090\pm 0.072$, $0.152\pm 0.100$, $0.390\pm 0.280$
Dynamic Error |  $0.108\pm 0.094$, $0.181\pm 0.156$, $1.355\pm 0.627$

###  III-B Force sensor characterization

The hysteresis evaluation of the custom standalone force sensor resulted in
values of 11.57%, 12.39%, and 13.58% in the X, Y, and Z axes, respectively.
Along the Z-axis, we observed an initial offset of 1.168 $\mathrm{N}$ under
no-load conditions, which increased to 2.717 $\mathrm{N}$ under the 12.806
$\mathrm{N}$ maximum load applied during the test. We also observed
viscoelastic creep of 0.090 $\mathrm{N}$, 0.058 $\mathrm{N}$, and 0.508
$\mathrm{N}$ in X, Y, and Z, respectively. Finally, the measurement errors in
the X, Y, and Z directions induced by the torsional distortions were -0.071
$\mathrm{N}$, 1.388 $\mathrm{N}$, and 4.926 $\mathrm{N}$ for rotations around
the X axis; 1.327 $\mathrm{N}$, 0.192 $\mathrm{N}$, and 6.538 $\mathrm{N}$ for
the Y axis; and 0.201 $\mathrm{N}$, -0.137 $\mathrm{N}$, and 0.843
$\mathrm{N}$ for the Z axis.

###  III-C Force control evaluation

####  III-C1 Custom embedded sensor

The calibration of the custom embedded sensor yielded maximum absolute forces
of approximately (6, 6, 6) $\mathrm{N}$ and test set static errors of
($0.214\pm 0.109$, $0.095\pm 0.105$, $1.451\pm 0.337$) $\mathrm{N}$ in the X,
Y, and Z axes, respectively.

####  III-C2 Step response

We commanded a series of step signals to the interface along different axes.
Regarding the steady-state error of the custom sensor w.r.t. the commanded
force, average errors of 9 $\mathrm{m}\mathrm{N}$, 28 $\mathrm{m}\mathrm{N}$,
and 37 $\mathrm{m}\mathrm{N}$ were observed in X, Y, and XY (diagonally),
respectively. The rise times, as measured by the custom sensor, were 0.063
$\mathrm{s}$, 0.065 $\mathrm{s}$, and 0.071 $\mathrm{s}$ in X, Y, and XY,
respectively.

####  III-C3 Frequency response

We commanded a chirp signal of constant amplitude along each system axis. In
the bode plots (Fig. 3), we observed gain-crossover frequencies at
approximately 4 $\mathrm{H}\mathrm{z}$ and 8 $\mathrm{H}\mathrm{z}$ for the X
and Y axes, respectively.

Figure 3: Bode plots of the tactile interface’s frequency response for the X
and Y axes.

##  IV Discussion

Despite the importance of tactile information in motor control, its inclusion
in robotic neurorehabilitation systems remains largely unexplored. Here, we
presented a novel two-dimensional skin stretch tactile interface to extend the
haptic rendering capabilities of the PRIDE hand rehabilitation robot to
simulate ADLs with high ecological validity.

###  IV-A Compact multi-finger tactile interface providing two-dimensional
skin stretch in a hand rehabilitation robot

Our four-finger tactile interface has the potential to enhance the sensory
stimulation capabilities of our hand rehabilitation robot PRIDE by providing
skin stretch on top of the already available kinesthetic feedback. This is an
advantage w.r.t. current technologies that usually only provide single-finger
tactile feedback [9, 11] and do not enable interactions with tangible virtual
objects, such as grasping [14].

The actuation mechanism of our interface achieves the required deformations of
the fingerpads and maximum renderable forces (5 $\mathrm{m}\mathrm{m}$ and 5
$\mathrm{N}$). However, the compliance of the custom force sensor reduced
these maximum values to an effective ROM of $\pm$3 $\mathrm{m}\mathrm{m}$ and
3 $\mathrm{N}$ lateral force, similar to other state-of-the-art devices (e.g.,
[9]). With the added capabilities, the full range of motion of PRIDE with
three DoFs was kept intact, except for small reductions in ROM for the thumb
flexion/extension and circumduction (by approximately 15
$\mathrm{\SIUnitSymbolDegree}$).

###  IV-B Small size magnetic-based custom 3D force sensor with large
measurement ranges

Our low-cost custom magnetic-based sensor recorded the largest shear force
ranges ($\pm$8 $\mathrm{N}$) for a sensor of this type and small size,
compared to those reported in literature (e.g., $\pm$1 $\mathrm{N}$ in [20,
21, 22]). Only one magnetic-based sensor was reported to achieve a similar
measurement range, yet it is ten times larger than ours [23]. The large
measurement-range-to-size ratio of our sensor can partly be attributed to the
stiffer elastomer selected. Furthermore, our experimental setup allowed the
application of forces with reduced slippage.

The static and dynamic errors of our sensors (1.56–4.81 %) were, in general,
higher than those found in the literature (e.g., 0.07–2.2 % in [20, 21, 22]).
This might be due to the limited number of data points employed for
calibration. Due to the time-consuming manual positioning stage in our
calibration procedure, we only recorded samples at 165 target positions for
the training of our magnetic-field-to-force mapping, a number substantially
smaller than those employed in previous studies that used automated
positioning stages, e.g., 4800 in [20] and 13000 in [21].

The observed force offset along the Z axis (approximately 7.3% of the maximum
force) could result from permanent elastomer deformations during calibration
and temperature changes affecting the elastomer’s stiffness. Le Signor et al.
report a 23% increase in measured force in the Z axis for a temperature
increase of 50 $\mathrm{\SIUnitSymbolDegree}\mathrm{C}$ [21]. We chose not to
incorporate temperature compensation as it increases the measurement time by
0.25–1.6 $\mathrm{m}\mathrm{s}$, depending on the oversampling ratio.

We measured a hysteresis of 11–13 % in all axes, in contrast to reported
values in literature, e.g., 3.4 % in [20] and no hysteresis in [22] in the
normal direction. There is, however, no comparable data from similar sensors
in the lateral directions. This difference might originate in the larger
deformations applied to our sensor and to the calibration database that only
contained measurements collected during static, step-wise loading.

The measurement of viscoelastic creep in the lateral directions is a novelty
for this type of sensor. We observed a force change of approximately 2 % of
the applied force in the X and Y axes and 5 % in the Z-axis. This is in line
with the 3 % creep measured in the normal direction by Wang et al., despite
the difference in materials, i.e., shore 00-30 hardness vs. A30 in our sensor
[20]. Notably, we observed that creep effects are less prominent in the shear
vs. the normal direction, likely due to the material properties of the
elastomer.

Finally, the applied torsional deformations around the X and Y axes caused
errors and drifts in the measured normal force. This result highlights the
importance of physically constraining torsional deformations, especially
around the X and Y axes, to achieve accurate force readings.

###  IV-C Closed-loop force rendering in a compact solution

The inclusion of the compact and low-cost force sensor allows for closed-loop
rendering of lateral forces at the fingertips. This addresses the limitations
of current systems that depend on individual calibration of finger pad
stiffness for accurate force rendering [11, 9] and systems requiring costly or
bulky force sensors [14].

The tactile interface’s step-response steady-state errors (0.6–2.5%) indicate
accurate force rendering of the system. The measured errors (9–37
$\mathrm{m}\mathrm{N}$) were smaller than the required maximum errors to be
below the JND for force perception (50–100 $\mathrm{m}\mathrm{N}$) [17].
However, it should be noted that the overall force rendering accuracy of the
system is limited by the accuracy of the custom force sensor. Yet, we argue
that the system’s lateral force rendering quality is an improvement to those
solutions relying on fingerpad stiffness estimations. Fingertip stiffnesses
can be as much as two times larger from one individual to the other, or for
one person between two directions, i.e., mediolateral vs. proximodistal [12].
Thus, our solution has the potential to eliminate the need for tedious
individual fingerpad stiffness calibration procedures.

The observed rise times (63–71$\mathrm{m}\mathrm{s}$) were smaller than those
measured for lifting onset during object manipulation (approximately 350
$\mathrm{m}\mathrm{s}$ [18]), fulfilling our requirement. However, the
effective bandwidth of the interface (4–8 $\mathrm{H}\mathrm{z}$ gain
crossover) was slightly lower than the targeted 10 $\mathrm{H}\mathrm{z}$ and
the effective bandwidth of comparable skin stretch devices (e.g., 10-12
$\mathrm{H}\mathrm{z}$ gain crossover [13]). This might be due to compliance
added from the deformable force sensor, actuation backlash, and the filtering
in the control strategy. The compliance of the force sensor acts as a
mechanical low-pass filter, limiting the effective bandwidth of the system and
adding lag. While the frequency response of the system is sufficient to render
a wide range of object interactions during ADLs [19], especially during
neurorehabilitation, a higher bandwidth would allow for accurate rendering of
extreme cases, such as collisions.

###  IV-D Limitations & Future Work

Most of the limitations of our tactile interface, such as the low accuracy of
the custom sensor and the limited frequency response, stem from our
requirement to create a low-cost and accessible system. Investing in high-
quality force sensors and actuators may result in a more accurate and stiff
system, and alleviate the need for lengthy calibration procedures.

The employed control strategy was a first effort to show the capabilities of
the developed system but still lacked refinement regarding filtering,
parameter tuning, and system modeling. Future work should contribute towards
increasing the system’s force rendering quality, i.e., the force rendering
accuracy and bandwidth, by improving or exchanging the sensor, actuation
mechanism, and control strategy. Finally, the effect of training with the
combined haptic feedback of PRIDE and the tactile interface on motor control
should be investigated with people with motor and sensory deficits.

##  V Conclusion

Existing robotic neurorehabilitation systems largely fail to provide tactile
information, a crucial component for motor control when interacting and
manipulating objects. We presented a novel low-cost tactile interface to
provide skin stretch stimulation in the mediolateral and proximodistal
directions of the little to index fingers to extend the capabilities of hand
rehabilitation robots to render interactions with virtual tangible objects as
a means to enhance the ecological validity of robotic neurorehabilitation
training. We manufactured a custom low-cost three-axis magnetic-based force
sensor with the highest measurement-range-to-size ratio of its kind to achieve
a closed-loop force control scheme. The forces rendered by the tactile
interface showed low steady-state errors and a frequency response sufficient
to simulate a wide range of activities of daily living. We hope that this
study serves as a stepping stone towards including tactile rendering in
robotic neurorehabilitation and, ultimately, enhance the recovery and quality
of life for people after a stroke.

## References

  * [1] R. Gassert and V. Dietz, “Rehabilitation robots for the treatment of sensorimotor deficits: a neurophysiological perspective,” _Journal of NeuroEngineering and Rehabilitation_ , vol. 15, no. 1, p. 46, Jun. 2018. [Online]. Available: https://doi.org/10.1186/s12984-018-0383-x
  * [2] E. Basalp, P. Wolf, and L. Marchal-Crespo, “Haptic training: Which types facilitate (re)learning of which motor task and for whom? answers by a review,” _IEEE Transactions on Haptics_ , vol. 14, pp. 722–739, 10 2021. [Online]. Available: https://ieeexplore.ieee.org/document/9513580/
  * [3] D. E. Levac, M. E. Huber, and D. Sternad, “Learning and transfer of complex motor skills in virtual reality: a perspective review,” _Journal of NeuroEngineering and Rehabilitation_ , vol. 16, no. 1, p. 121, Oct. 2019. [Online]. Available: https://doi.org/10.1186/s12984-019-0587-8
  * [4] K. Lo, M. Stephenson, and C. Lockwood, “Effectiveness of robotic assisted rehabilitation for mobility and functional ability in adult stroke patients: a systematic review,” _JBI database of systematic reviews and implementation reports_ , vol. 15, pp. 3049–3091, 12 2017. [Online]. Available: https://journals.lww.com/jbisrir/fulltext/2017/12000/effectiveness_of_robotic_assisted_rehabilitation.15.aspx
  * [5] J. M. Walker, A. A. Blank, P. A. Shewokis, and M. K. O’Malley, “Tactile feedback of object slip facilitates virtual object manipulation,” _IEEE Transactions on Haptics_ , vol. 8, pp. 454–466, 10 2015. [Online]. Available: https://ieeexplore.ieee.org/document/7080902/
  * [6] F. Binkofski, E. Kunesch, J. Classen, R. J. Seitz, and H. J. Freund, “Tactile apraxia: Unimodal apractic disorder of tactile object exploration associated with parietal lobe lesions,” _Brain_ , vol. 124, pp. 132–144, 1 2001. [Online]. Available: https://dx.doi.org/10.1093/brain/124.1.132
  * [7] A. Grundmann, M. Hardwick, D. Ledingham, and J. Miller, “Sensory ataxia with cranial nerve palsies,” _Practical Neurology_ , vol. 22, no. 1, pp. 85–89, 2022. [Online]. Available: https://pn.bmj.com/content/22/1/85
  * [8] W. Provancher and N. Sylvester, “Fingerpad skin stretch increases the perception of virtual friction,” _IEEE Transactions on Haptics_ , vol. 2, pp. 212–223, 10 2009. [Online]. Available: http://ieeexplore.ieee.org/document/5196670/
  * [9] S. B. Schorr and A. M. Okamura, “Fingertip Tactile Devices for Virtual Object Manipulation and Exploration,” in _Proceedings of the 2017 CHI Conference on Human Factors in Computing Systems_ , ser. CHI ’17. New York, NY, USA: Association for Computing Machinery, May 2017, pp. 3115–3119. [Online]. Available: https://doi.org/10.1145/3025453.3025744
  * [10] Z. F. Quek, S. B. Schorr, I. Nisky, W. R. Provancher, and A. M. Okamura, “Sensory substitution of force and torque using 6-dof tangential and normal skin deformation feedback,” in _2015 IEEE International Conference on Robotics and Automation (ICRA)_ , vol. 2015-June. IEEE, 5 2015, pp. 264–271. [Online]. Available: http://ieeexplore.ieee.org/document/7139010/
  * [11] D. Leonardis, M. Solazzi, I. Bortone, and A. Frisoli, “A 3-rsr haptic wearable device for rendering fingertip contact forces,” _IEEE Transactions on Haptics_ , vol. 10, pp. 305–316, 7 2017. [Online]. Available: https://ieeexplore.ieee.org/document/7784835/
  * [12] M. Wiertlewski and V. Hayward, “Mechanical behavior of the fingertip in the range of frequencies and displacements relevant to touch,” _Journal of Biomechanics_ , vol. 45, no. 11, pp. 1869–1874, Jul. 2012. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0021929012003260
  * [13] Y. Kamikawa and A. M. Okamura, “Comparison between force-controlled skin deformation feedback and hand-grounded kinesthetic force feedback for sensory substitution,” _IEEE Robotics and Automation Letters_ , vol. 3, pp. 2174–2181, 7 2018. [Online]. Available: http://ieeexplore.ieee.org/document/8304766/
  * [14] F. E. van Beek, R. J. King, C. Brown, M. D. Luca, and S. Keller, “Static weight perception through skin stretch and kinesthetic information: Detection thresholds, jnds, and pses,” _IEEE Transactions on Haptics_ , vol. 14, pp. 20–31, 1 2021. [Online]. Available: https://ieeexplore.ieee.org/document/9141483/
  * [15] R. Ratz, R. M. Muri, and L. Marchal-Crespo, “Design of a haptic palmar device with thumb flexion and circumduction movements for sensorimotor stroke rehabilitation,” in _2022 44th Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC)_, vol. 2022-July. IEEE, 7 2022, pp. 2644–2647. [Online]. Available: https://ieeexplore.ieee.org/document/9871303/
  * [16] N. Nakazawa, R. Ikeura, and H. Inooka, “Characteristics of human fingertips in the shearing direction,” _Biological Cybernetics_ , vol. 82, no. 3, pp. 207–214, Feb. 2000. [Online]. Available: https://doi.org/10.1007/s004220050020
  * [17] E. E. Brodie and H. E. Ross, “Sensorimotor mechanisms in weight discrimination,” _Perception & Psychophysics_, vol. 36, no. 5, pp. 477–481, Sep. 1984. [Online]. Available: https://doi.org/10.3758/BF03207502
  * [18] R. S. Johansson and G. Westling, “Roles of glabrous skin receptors and sensorimotor memory in automatic control of precision grip when lifting rougher or more slippery objects,” _Experimental Brain Research_ , vol. 56, no. 3, pp. 550–564, Oct. 1984. [Online]. Available: https://doi.org/10.1007/BF00237997
  * [19] K. A. Mann, F. W. Wernere, and A. K. Palmer, “Frequency spectrum analysis of wrist motion for activities of daily living,” _Journal of Orthopaedic Research_ , vol. 7, pp. 304–306, 3 1989. [Online]. Available: https://onlinelibrary.wiley.com/doi/10.1002/jor.1100070219
  * [20] H. Wang, G. De Boer, J. Kow, A. Alazmani, M. Ghajari, R. Hewson, and P. Culmer, “Design Methodology for Magnetic Field-Based Soft Tri-Axis Tactile Sensors,” _Sensors_ , vol. 16, no. 9, p. 1356, Sep. 2016. [Online]. Available: https://www.mdpi.com/1424-8220/16/9/1356
  * [21] T. L. Signor, N. Dupre, and G. F. Close, “A gradiometric magnetic force sensor immune to stray magnetic fields for robotic hands and grippers,” _IEEE Robotics and Automation Letters_ , vol. 7, pp. 3070–3076, 4 2022. [Online]. Available: https://ieeexplore.ieee.org/document/9695271/
  * [22] A. Dwivedi, A. Ramakrishnan, A. Reddy, K. Patel, S. Ozel, and C. D. Onal, “Design, Modeling, and Validation of a Soft Magnetic 3-D Force Sensor,” _IEEE Sensors Journal_ , vol. 18, no. 9, pp. 3852–3863, May 2018. [Online]. Available: https://ieeexplore.ieee.org/abstract/document/8314107
  * [23] Q. Nie and F. C. Sup, “A Soft Four Degree-of-Freedom Load Cell Based on the Hall Effect,” _IEEE Sensors Journal_ , vol. 17, no. 22, pp. 7355–7363, Nov. 2017. [Online]. Available: https://ieeexplore.ieee.org/abstract/document/8047967


# Spatial cage solitons — taming light bullets

Chao Mei  Max Born Institute for Nonlinear Optics and Short Pulse
Spectroscopy, 12489 Berlin, Germany
and School of Computer and Communication Engineering, University of Science
and Technology Beijing (USTB), 100083 Beijing, China  Ihar Babushkin
Institute of Quantum Optics, Leibniz University Hannover, 30167 Hannover,
Germany  Tamas Nagy  Max Born Institute for Nonlinear Optics and Short Pulse
Spectroscopy, 12489 Berlin, Germany  Günter Steinmeyer<EMAIL_ADDRESS>
Max Born Institute for Nonlinear Optics and Short Pulse Spectroscopy, 12489
Berlin, Germany
and Institut für Physik, Humboldt Universität zu Berlin, 12489 Berlin, Germany

###### Abstract

Multimode nonlinear optics offers to overcome a long-standing limitation of
fiber optics, tightly phase locking several spatial modes and enabling the
coherent transport of a wavepacket through a multimode fiber. A similar
problem is encountered in the temporal compression of multi-mJ pulses to few-
cycle duration in hollow gas-filled fibers. Scaling the fiber length to up to
six meters, hollow fibers have recently reached 1 TW of peak power. Despite
the remarkable utility of the hollow fiber compressor and its widespread
application, however, no analytical model exists to enable insight into the
scaling behavior of maximum compressibility and peak power. Here we extend a
recently introduced formalism for describing mode-locking to the spatially
analogue scenario of locking spatial fiber modes together. Our formalism
unveils the coexistence of two soliton branches for anomalous modal dispersion
and indicates the formation of stable spatio-temporal light bullets that would
be unstable in free space, similar to the temporal cage solitons in mode-
locking theory. Our model enables deeper understanding of the physical
processes behind the formation of such light bullets and predict the existence
of multimode solitons in a much wider range of fiber types than previously
considered possible.

Spatial solitons have fascinated researchers since the early days of nonlinear
optics [1, 2]. In combination with self-phase modulation, the self-focusing
effect offers the possibility for three-dimensional contraction of an optical
wavepacket and concomitant intensity increase. While there exist numerous
reports on such light bullets [3, 4, 5], this intriguing nonlinear mechanism
found little application, probably because of the limiting action of a spatial
modulation instability [6, 7]. Starting from small imperfection in the beam
profile, this process induces a rapid small-scale breakup of the beam profile
into filaments when the critical power $P_{\rm crit}$ is exceeded [8, 9], thus
limiting the obtainable nonlinear interaction length. Given this severe
constraint, high power pulse compression and nonlinear conversion techniques
have resorted to hollow capillaries for extended nonlinear interaction length
[10, 11, 12, 13]. While several other competing techniques [14, 15] have been
discussed for the compression of pulses with gigawatt peak powers, the hollow
fiber is currently the most established compression technique and found
widespread application in attosecond pulse generation and other high-field
experiments, see, e.g., [16, 17]. Utilizing the advanced stretched-fiber
technique [18, 19] recently enabled record-breaking continuous powers above
300 W [20] and peak powers exceeding 1 TW for the first time [21]. Despite the
widespread utility of this technique, however, there exist only relatively few
analytical approaches [22, 23] for modeling the nonlinear broadening processes
inside the hollow fiber. Numerical simulations often resorted to a simplified
one-dimensional approximation as full modal expansions [24, 25, 26] are
numerically cumbersome. In the following, we present a completely analytical
approach for determining spatial soliton solutions in nonlinear multimode
fiber geometries. Similar spatio-temporal solitons have previously been
observed in numerical simulations [27, 28]. Moreover, our approach is
mathematically similar to the cage soliton solutions of the Haus master
equation of mode-locking [29, 30], and it is also applicable for the thriving
field of multimode fiber nonlinear optics [31, 32, 33]. Assuming adiabatic
pulse shaping, the results of this analysis enable the derivation of universal
scaling laws for the design of nonlinear multimode waveguides.

Figure 1: (a) Ray-optical representation of hollow-fiber transmission. The
wavevector $k_{0}$ can be decomposed into a transverse component $k_{\perp}$
and a longitudinal component $\beta$, which are connected by Pythagoras’
Theorem [34]. (b) The propagation constants $\beta_{n}$ of the individual EH1n
modes follow an approximate $n^{2}$ dependence whereas the $k_{n}^{\perp}$
underly a linear relationship with $n$. (c) Mode fields of the first four EH1n
modes considered in this study. Intensities are depicted by colors and
electric fields are represented by arrows.

The linear optical properties of a cylindrical hollow dielectric waveguide
with radius $a$ have been first modeled by Marcatili and Schmeltzer [10].
Assuming linear polarization, one finds hybrid solutions of the wave equation,
which have originally been designated as ${\rm EH}_{mn}$ modes. In solid core
multimode (SCM) fibers, widely similar solutions are referred to as ${\rm
LP}_{(m-1)n}$ modes. Strictly speaking, both designations are not identical
because of the Goos-Hänchen effect [35], which causes the ${\rm LP}_{(m-1)n}$
modes to extend into the cladding. In contrast, the ${\rm EH}_{mn}$ exhibit a
node of the electric field at the dielectric interface. Assuming $m=1$, i.e.,
azimuthal homogeneity, the radial field profile of the ${\rm EH}_{0n}$ mode is
given by

$E_{n}(r)\propto J_{0}\left(\frac{u_{n}r}{a}\right),$ (1)

with the zero-order Bessel function $J_{0}$ and its $n$th zero $u_{n}$
($n>0$). The complex-valued propagation constant of these modes is given by

$\kappa_{n}=\beta_{n}+i\alpha_{n}=k_{0}-\frac{u_{n}^{2}}{2k_{0}a^{2}}\left(1-i\frac{\epsilon+1}{k_{0}a\sqrt{\epsilon-1}}\right),$
(2)

where $k_{0}=2\pi/\lambda$ is the wavenumber and $\sqrt{\epsilon}=n_{\rm
clad}$ the refractive index of the cladding material. Exploiting the
trigonometric relationship between propagation constant and wave number
displayed in Fig. 1, one finds an approximate parabolic dependence [34, 36]

$\displaystyle\beta_{n}=k_{0}\cos\left(\frac{u_{n}}{k_{0}a}\right)$
$\displaystyle\approx$
$\displaystyle\beta_{1}-\frac{\pi^{2}(n-1)^{2}}{2k_{0}a^{2}}$ (3)
$\displaystyle=$ $\displaystyle\beta_{1}+\mathscr{B}(n-1)^{2}.$

This approximation may appear crude for small $n$ and mostly serves to point
out the close analogy to the role of dispersion in the propagation of
femtosecond pulses. Therefore, all subsequent computations have been repeated
with and without the approximation, yet with little resulting difference in
the outcome. Moreover, as $\beta_{n}<k_{0}$ in a hollow waveguide, the modes
$E_{n}$ propagate at superluminal phase velocity inside, which can be
explained by $-\pi$ phase jumps due to Fresnel reflection at the interface in
the geometrical picture of Refs. [34, 36]. A similar parabolic approximation
can be made for SCM fibers [37], yet with positive curvature $\mathscr{B}$.
Exploiting the relation

$\beta^{2}+k_{\perp}^{2}=k_{0}^{2}\sin^{2}\theta+k_{0}^{2}\cos^{2}\theta=k_{0}^{2}$
(4)

depicted in Fig. 1(a), one also finds an approximate linear relationship for
the transverse wave number

$k_{n}^{(\perp)}=k_{0}\sin\left(\frac{u_{n}}{k_{0}a}\right)\approx\left(n-\frac{1}{4}\right)\frac{\pi}{a}.$
(5)

In the following, we describe the evolution of the spatial beam profile
$E(z,r)$ upon propagation along the coordinate $z$ with the transverse wave
equation

$\partial_{z}E=\frac{i\beta}{r}\partial_{r}r\partial_{r}E+i\,\Gamma|E|^{2}E,$
(6)

due to competing diffractive effects and the action of self-focusing with
$\Gamma=k_{0}n_{2}$ and the nonlinear refractive index $n_{2}$. Equation (6)
is widely similar to the Nonlinear Schrödinger Equation (NLSE), yet with the
temporal coordinate $t$ replaced by the radial coordinate $r$. While it is
customary to numerically solve the NLSE by a split-step Fourier method [38],
the equation can be reformulated in the frequency domain [30], which leads to
a system of nonlinearly coupled ordinary differential equations. Rewriting Eq.
(6) in the $k^{(\perp)}$ domain leads to a similar system of coupled equations
[27]

$\partial_{z}\tilde{E}_{n}=i\beta_{n}\tilde{E}_{n}+i\,\Gamma\mkern-18.0mu\sum_{j+k-\ell=n}\mkern-18.0mu\tilde{E}_{j}\tilde{E}_{k}\tilde{E}^{*}_{\ell}.$
(7)

Here the summation over the four-wave mixing terms is restricted to those
products that satisfy conservation of the transverse wavenumber $k_{\perp}$,
cf. Fig.1(b). The partially degenerate case $j=k$ is also known as self-
diffraction [39]. As was noted by DeLong et al. [40], the self-diffraction
process is not exactly phase matched. This slight mismatch is accounted for by
the intermodal dispersion term $\beta_{n}\tilde{E}_{n}$, which takes the part
of group-delay dispersion with its quadratic dependence in the frequency
representation of the Haus master equation [30]. As the parabolic dependence
on mode number $n$ is convex for hollow fibers one can associate this case
with anomalous dispersion whereas SCM fibers display normal modal dispersion.
Moreover, the slight phase mismatch can be understood with the neglection of
the linear term in Eq. (3) and is otherwise similar to higher-order dispersion
contributions in traditional mode-locking theory [29], i.e., these
contributions give rise to deviations from equidistance of the cold cavity
modes. Solving the mode-locking version of Eqs. (LABEL:eq:ODE2) and (10)
requires restrictive assumptions on the number of coupled longitudinal modes
[30]. In the transverse case, it is often considered sufficient to include
only a few spatial modes (e.g., $N=3$ [26]) for treating nonlinear propagation
through a hollow waveguide. We therefore write out Eq. (7) for 4 modes, using
excessive content in $\tilde{E}_{4}$ as an indicator for the breakdown of our
simplifying assumptions

$\displaystyle\partial_{z}\tilde{E}_{1}$ $\displaystyle=$ $\displaystyle
i\beta_{1}\tilde{E}_{1}+i\Gamma\left[\eta_{1234}\tilde{E}_{2}\tilde{E}_{3}\tilde{E}_{4}^{*}+\eta_{123}\tilde{E}_{2}^{2}\tilde{E}_{3}^{*}+\left(|\tilde{E}_{1}|^{2}+\eta_{12}|\tilde{E}_{2}|^{2}+\eta_{13}|\tilde{E}_{3}|^{2}+\eta_{14}|\tilde{E}_{4}|^{2}\right)\tilde{E}_{1}\right]$
$\displaystyle\partial_{z}\tilde{E}_{2}$ $\displaystyle=$ $\displaystyle
i\beta_{2}\tilde{E}_{2}+i\Gamma\left[\eta_{1234}\tilde{E}_{1}\tilde{E}_{3}^{*}\tilde{E}_{4}+\eta_{123}\tilde{E}_{1}\tilde{E}_{2}^{*}\tilde{E}_{3}+\eta_{234}\tilde{E}_{3}^{2}\tilde{E}_{4}^{*}\right.$
$\displaystyle+\left.\left(\eta_{12}|\tilde{E}_{1}|^{2}+|\tilde{E}_{2}|^{2}+\eta_{23}|\tilde{E}_{3}|^{2}+\eta_{24}|\tilde{E}_{4}|^{2}\right)\tilde{E}_{2}\right]$
$\displaystyle\partial_{z}\tilde{E}_{3}$ $\displaystyle=$ $\displaystyle
i\beta_{3}\tilde{E}_{3}+i\Gamma\left[\eta_{1234}\tilde{E}_{1}\tilde{E}_{2}^{*}\tilde{E}_{4}+\eta_{123}\tilde{E}_{1}^{*}\tilde{E}_{2}^{2}+\eta_{234}\tilde{E}_{2}\tilde{E}_{3}^{*}\tilde{E}_{4}\right.$
$\displaystyle+\left.\left(\eta_{13}|\tilde{E}_{1}|^{2}+\eta_{23}|\tilde{E}_{2}|^{2}+|\tilde{E}_{3}|^{2}+\eta_{34}|\tilde{E}_{4}|^{2}\right)\tilde{E}_{3}\right]$
$\displaystyle\partial_{z}\tilde{E}_{4}$ $\displaystyle=$ $\displaystyle
i\beta_{4}\tilde{E}_{4}+i\Gamma\left[\eta_{1234}\tilde{E}_{1}^{*}\tilde{E}_{2}\tilde{E}_{3}+\eta_{234}\tilde{E}_{2}^{*}\tilde{E}_{3}^{2}\left(\eta_{14}|\tilde{E}_{1}|^{2}+\eta_{24}|\tilde{E}_{2}|^{2}+\eta_{34}|\tilde{E}_{3}|^{2}+|\tilde{E}_{4}|^{2}\right)\tilde{E}_{4}\right].$

Here the $\beta_{n}$ terms describe intermodal dispersion, i.e., the phase
velocity differences between the individual modes. Self-diffraction is
accounted for by partially degenerate four-wave mixing terms
$\propto\eta_{k\ell}\tilde{E}^{2}_{k}E^{*}_{\ell}$. Self-focusing of the
individual modes is described by the fully degenerate terms
$\propto|E_{k}|^{2}E_{k}$. In addition, cross-phase modulation (XPM) terms
appear, cf. Table 1. In Eq. (LABEL:eq:ODE2), we introduced modal overlap
factors $\eta_{nj}$, $\eta_{njk}$ and $\eta_{njk\ell}$ for non-collinear four-
wave mixing processes as they have been previously discussed, e.g., by Chapman
et al. [41]. For the fully non-degenerate process we define

$\eta_{njk\ell}=\frac{\int\limits_{0}^{a}\tilde{E}_{n}\tilde{E}_{j}\tilde{E}_{k}\tilde{E}_{\ell}r{\rm
d}r}{\prod\limits_{m=\\{n,j,k,\ell\\}}\sqrt[4]{\int_{0}^{a}\tilde{E}_{m}^{4}r{\rm
d}r}}.$ (9)

For degenerate mixing processes, we use a shorthand notation, e.g.,
$\eta_{1221}=\eta_{2112}=\eta_{12}/2$ and include the degeneracy factor in the
respective overlap factors. Values for the various $\eta$ are listed in Table
1. Provided that the nonlinear length $L_{\rm NL}$ is much shorter than the
dispersive length $L_{\rm D}$ [38], Eq. (LABEL:eq:ODE2) can be used as a
highly efficient tool for simulating the propagation via solving a set of
coupled ordinary differential equations. As will be further discussed below,
such adiabaticity can be assumed in typical hollow-fiber compression
scenarios. Moreover, propagation losses can be accounted for by a complex-
valued redefinition of the $\beta_{n}$.

Assuming propagation of the $\tilde{E}_{j}$ at identical phase velocity, i.e.,
as solitonic wavepacket, we find a wavenumber offset $\psi$ relative to the
fundamental mode $\beta_{1}$. We further renormalize real-valued electric
field amplitudes $a_{n}=|\tilde{E}_{n}/\tilde{E}_{1}|$ to yield $a_{1}=1$ and
redefine an effective nonlinearity
$\gamma=\Gamma\mathscr{B}^{-1}|\tilde{E}_{1}|^{-2}$. Using these
simplifications, we can extract an algebraic discriminant for the resulting
spatial soliton

$\displaystyle\psi$ $\displaystyle=$
$\displaystyle\gamma\left(1+\eta_{1234}a_{2}a_{3}a_{4}+\eta_{123}a_{2}^{2}a_{3}+\eta_{12}a_{2}^{2}+\eta_{13}a_{3}^{2}+\eta_{14}a_{4}^{2}\right)$
(10) $\displaystyle=$
$\displaystyle\gamma\left[\eta_{1234}a_{3}a_{4}+\eta_{123}a_{2}a_{3}+\eta_{234}a_{3}^{2}a_{4}+\left(\eta_{12}+a_{2}^{2}+\eta_{23}a_{3}^{2}+\eta_{24}a_{4}^{2}\right)a_{2}\right]+a_{2}$
$\displaystyle=$
$\displaystyle\gamma\left[\eta_{1234}a_{2}a_{4}+\eta_{123}a_{2}^{2}+\eta_{234}a_{2}a_{3}a_{4}+\left(\eta_{13}+\eta_{23}a_{2}^{2}+a_{3}^{2}+\eta_{34}a_{4}^{2}\right)a_{3}\right]+4a_{3}$
$\displaystyle=$
$\displaystyle\gamma\left[\eta_{1234}a_{2}a_{3}+\eta_{234}a_{2}a_{3}^{2}+\left(\eta_{14}+\eta_{24}a_{2}^{2}+\eta_{34}a_{3}^{2}+a_{4}^{2}\right)a_{4}\right]+9a_{4}.$

Table 1: Nonlinear mode coupling factors. XPM | $\eta_{12}$ | = | 1.709 | $\eta_{13}$ | = | 1.529 | $\eta_{14}$ | = | 1.408
---|---|---|---|---|---|---|---|---|---
XPM | $\eta_{23}$ | = | 1.847 | $\eta_{24}$ | = | 1.725 | $\eta_{34}$ | = | 1.897
FWM | $\eta_{123}$ | = | 0.788 | $\eta_{234}$ | = | 0.874 | $\eta_{1234}$ | = | 1.462

One can now retrieve all real-valued roots of Eq. (10) and compute beam
diameter $w_{\rm eff}$, loss $\alpha_{\rm eff}$, and, most importantly, the
resulting “soliton phase” $\psi$ [42], which is more correctly defined as the
total propagation constant involving both linear and nonlinear effects
relative to the phase velocity frame. $\psi$ vanishes for $\gamma\rightarrow
0$, as we already accounted for linear propagation effects by subtracting the
propagation constant $\beta_{1}$ in Eq. (10). This effective subtraction of
the linear propagation phase now gives immediate access to the nonlinear phase
via $\partial\varphi_{\rm nl}/\partial z=\psi$. Using the relation
$\varphi_{\rm nl}=k_{0}zn_{2}PA_{\rm eff}^{-1}$ with the nonlinear refractive
index $n_{2}$ and the effective area [38] of the fundamental EH11 mode

$A_{\rm eff}=\frac{1}{2\pi}\frac{\left(\int|E(r)|^{2}r\rm
dr\right)^{2}}{\int|E(r)|^{4}r\rm dr}\approx 1.50a^{2},$ (11)

one can then relate the soliton phase to the power $P$ in the nonlinear
waveguide. Inserting the definition of the critical power for self-focusing
$P_{\rm cr}=0.147\lambda^{2}/n_{2}$ [9, 8], we then finally yield

$P\approx 10.2P_{\rm cr}\frac{\psi a^{2}}{k_{0}\lambda^{2}}.$ (12)

Figure 2: Spatial soliton solution branches of Eq. (LABEL:eq:ODE2). For normal
modal dispersion ($n_{\rm core}>n_{\rm clad}$), a single solution branch
exists (green). In hollow fibers, two branches coexist (blue and red). The red
branch is considered unstable (see discussion in text). (a) Radially
integrated intensity $\int E(r)^{2}r{\rm d}r=P$ of the spatial cage soliton
solutions vs. effective nonlinearity. Powers have been normalized to the
critical power of self-focusing $P_{\rm cr}$ in free space [9]. Insets show
$E(r)$ for parameters indicated by symbols. (b) Root mean square mode diameter
of spatial solitons normalized to the HE11 mode. Insets show spatial intensity
profiles $|E(r)|^{2}$.

The resulting peak power vs. effective nonlinearity is depicted in Fig. 2(a).
In the anomalous regime (hollow fibers), the higher-order modes ${\rm
EH}_{1n},n>1$ of the spatial soliton solutions carry the opposite sign of the
fundamental ${\rm EH}_{11}$ mode. This can be understood by a partial
cancellation of nonlinearity and modal dispersion, similar to the formation of
Schrödinger solitons [2]. In turn, a depression of the peak intensity and an
increased effective area $A_{\rm eff}$ of the beam profile result. [Fig.
2(b)]. Consequently, the hollow waveguide can host beams with a peak power of
$\approx 1.4P_{\rm cr}$, yet at a characteristic donut profile, cf. Fig. 3. In
the normal dispersion regime (SCM fibers), peak intensities are enhanced as
nonlinearity and modal dispersion add up, leading to cusp-like beam profiles
for $\kappa>0$. Comparing both waveguide dispersion regimes, it is striking
that the asymmetry in the $P(\gamma)$ relation is only caused by beam profile
variations. For the SCM case, this effect leads to a maximum peak power
hosting of 0.5 $P_{\rm cr}$. Apart from these two fundamental soliton
branches, our investigations identifies a second solution branch in the
anomalous dispersion regime (red curves in Fig. 2). In this branch, the limit-
value donut solution converges toward a “Mexican hat” solution upon subsequent
reduction of power. As the latter solutions involve an increasing amount of
higher-order modes, such beam profiles would experience rapidly increasing
losses, i.e., this solution branch is considered unstable. In contrast, the
fundamental soliton branches (blue and green in Fig. 2) self-stabilize upon
propagation as their higher-order contents reduces with decreasing peak power.

Figure 3: Three-dimensional visualization of spatio-temporal soliton structure
at the stability limit ($\approx 1.4P_{\rm crit}$) in the anomalous modal
dispersion regime. Equi-intensity surfaces are shown with colors red ($80\%$
peak int.) to blue ($10\%$ peak int.). In the center, a donut structure
dominates which evolves into an ellipsoidal shape in the temporal wings. The
glass-gas interface of the hollow fiber is depicted in light gray for
comparison.

To this end, it appears illustrative to compute the relevant interaction
lengths [38] inside a hollow fiber. Given the rather low group-velocity
dispersion of noble gases like argon, which may be additionally cancelled out
by the waveguide (group-velocity) dispersion of the hollow fiber, dispersion
lengths $L_{D}\approx\tau^{2}/|\beta_{2}|$ for pulse durations
$\tau_{0}>20\,$fs generally exceed the fiber length by a large factor for
experimental conditions in [19, 20, 21]. The absorption length is commonly
also chosen longer than the actual fiber length. This is contrasted by the
nonlinear length $L_{\rm NL}=\Gamma^{-1}P^{-1}$, which amounts to only a few
millimeters for powers $P>0.1P_{\rm cr}$ at identical experimental conditions.
One can therefore conclude that spatial soliton effects strongly dominate the
nonlinear dynamics inside the hollow fiber, causing an adiabatic reshaping of
the solitons as a reaction to the comparatively slow waveguide losses. Non-
solitonic contents is stripped off into linearly propagating higher-order
modes, which travel at reduced group velocities. For example, for
$a=200\,\mu$m and $\lambda=1\,\mu$m, the linear group delay of the ${\rm
EH}_{12}$ relative to the fundamental mode amounts to about 100 fs per meter
propagation length, i.e., the non-solitonic contents will lead to the
formation of a temporal continuum background after recompression. Ignoring
reshaping effects due to group-velocity dispersion, one can now compute the
structure of the emerging spatio-temporal light bullets, see Fig. 3. Here we
have chosen the highest possible peak power in Fig. 2, i.e., $P=1.4\,P_{\rm
cr}$, which lead to the formation of a donut spatial structure at pulse
center. At lower intensities this structure goes over into the more common
ellipsoidal shape of conventional light bullets.

Figure 4: Comparison of model results with measured data. (a) Total losses
(linear and nonlinear) vs. ratio of $a^{3}$ and $\lambda^{2}$ (curve and red
dots). Early measurements with relatively short hollow fibers [19] indicated
significantly higher losses whereas more recent measurements [20, 21] showed
excellent agreement (blue triangles). (b) Maximum beneficial length (solid
curve and solid symbols) and maximum compressibility (dashed line and hollow
symbols), cf. Eq. (13). This analysis confirms that superior compression can
be reached with longer hollow fibers and larger core diameters.

Utilizing the spatial soliton solutions of Eq. (LABEL:eq:ODE2), one can now
derive a few design rules for hollow fiber compressors. As nonlinearly
broadened spectra typically exhibit near-perfect spectral symmetry, one can
employ the simple relation [38]

$\chi(z)=\frac{\Delta\nu_{\rm
rms}}{\Delta\nu_{0}}=\sqrt{1+\frac{4}{3\sqrt{3}}\varphi_{\rm nl}^{2}(z)}$ (13)

to estimate the compressibility of the input pulse. In view of applications,
increase of peak power is typically considered more important than ultimate
shortness of the pulse. Let us therefore define the figure of merit

$M(a,\lambda)=\left.{\rm max}\left[\chi(z)\exp(-2\alpha_{\rm
tot}z)\right]\right|_{z}$ (14)

as the criterion for the maximum beneficial propagation length $z_{\rm max}$
inside the hollow fiber. Neglecting any possible loss in the subsequent
compression process (e.g., by chirped mirrors), peak powers will not further
increase upon additional propagation. Pertinent computations are shown in Fig.
4, indicating the dependence of the total (linear and nonlinear) loss
$\alpha_{\rm tot}$ as a function of $a^{3}/\lambda^{2}$, see Fig. 4(a).
Comparing with experimental data [blue symbols in Fig. 4(a)], it appears
striking that early work with relatively short fibers [19] reported losses
that deviate from predictions of our model whereas more recent reports with
longer fibers [20, 21] appear to completely agree with our model. Moreover, we
also compare the maximum beneficial length and the observed spectral
broadening between theory and previous experimental findings in Fig. 4(b).
Here it is not overly surprising that our projections are too optimistic, with
maximum compressibilities and beneficial fiber lengths that are about two
thirds of our predictions. Nevertheless, the relations derived from our
spatial soliton model clearly explain the trends observed in previous
experimentation, confirming that long hollow fibers promise superior
performance compared to the single-meter-long segments of early
experimentation.

In conclusion, the field of multimode nonlinear optics bears a number of
appealing applications, which are mostly ruled by a $a^{3}/\lambda^{2}$
relationship. In particular, in hollow fibers, where group-velocity dispersion
plays an inferior role, spatial soliton formation appears to take a previously
unrecognized lead role at large diameters and lengths. Given that dispersive
and absorptive lengths are orders of magnitude larger than the soliton length,
adiabatic reshaping dominates the nonlinear dynamics of pulse propagation
through the hollow fiber. As the critical power $P_{\rm cr}$ plays a decisive
limiting role for the performance of a nonlinear multimode waveguide, further
upscaling of compressible peak powers require usage of lower pressures or
gases with lower refractive index than the commonly used argon. As the
currently demonstrated highest peak powers already used fibers with several
meter length, hosting even higher powers requires the use of significantly
longer fibers of tens or even hundreds of meters length to accumulate
sufficient nonlinear phase for the broadening process. Such dimensions appear
to be out of range for universities but could certainly be implemented in
large scale facilities, in particular linear accelerators. Using SCM fibers
near their zero-dispersion wavelength instead, spectral broadening can be
accomplished at much higher pulse energies than in single-mode fibers. In
contrast to previous demonstrations of nonlinear multimode optics, our
theoretical investigations suggest that the exact refractive index profile
plays only a minor rule, enabling the use of simple step-index architectures
rather than relying only on parabolic profiles.

Acknowledgments. GS gratefully acknowledges fruitful discussions with Pavel
Sidorenko and Frank Wise (Cornell University) as well as with Howard Milchberg
(UMD).

## References

  * Chiao  _et al._ [1965] R. Y. Chiao, E. Garmire, and C. H. Townes, Self-trapping of optical beams, Phys. Rev. Lett. 14, 479 (1965).
  * Zakharov and Shabat [1972] V. E. Zakharov and A. B. Shabat, Exact theory of two-dimensional self-focusing and one-dimensional self-modulation of waves in nonlinear media, Sov. Phys. JETP 34, 62 (1972).
  * Silberberg [1990] Y. Silberberg, Collapse of optical pulses, Opt. Lett. 15, 1282 (1990).
  * Moll  _et al._ [2003] K. D. Moll, A. L. Gaeta, and G. Fibich, Self-similar optical wave collapse: observation of the Townes profile, Phys. Rev. Lett. 90, 203902 (2003).
  * Majus  _et al._ [2014] D. Majus, G. Tamošauskas, I. Gražulevičiūtė, N. Garejev, A. Lotti, A. Couairon, D. Faccio, and A. Dubietis, Nature of spatiotemporal light bullets in bulk Kerr media, Phys. Rev. Lett. 112, 193901 (2014).
  * Bespalov and Talanov [1966] V. I. Bespalov and V. I. Talanov, Filamentary structure of light beams in nonlinear liquids, JETP Lett. 3, 307 (1966).
  * Bliss  _et al._ [1974] E. S. Bliss, D. R. Speck, J. F. Holzrichter, J. H. Erkkila, and A. J. Glass, Propagation of a high-intensity laser pulsewith small-scale intensity modulation, Appl. Phys. Lett. 25, 448 (1974).
  * Marburger [1975] J. Marburger, Self-focusing: Theory, Progr. Quantum Electron. 4, 35 (1975).
  * Boyd [2020] R. W. Boyd,  _Nonlinear Optics_ , 4th ed. (Academic Press, London, UK, 2020).
  * Marcatili and Schmeltzer [1964] E. Marcatili and R. Schmeltzer, Hollow metallic and dielectric waveguides for long distance optical transmission and lasers, Bell System Tech. J. 43, 1783 (1964).
  * Nisoli  _et al._ [1996] M. Nisoli, S. De Silvestri, and O. Svelto, Generation of high energy 10 fs pulses by a new pulse compression technique, Appl. Phys. Lett. 68, 2793 (1996).
  * Nisoli  _et al._ [1997] M. Nisoli, S. De Silvestri, O. Svelto, R. Szipöcs, K. Ferenz, C. Spielmann, S. Sartania, and F. Krausz, Compression of high-energy laser pulses below 5 fs, Opt. Lett. 22, 522 (1997).
  * Durfee  _et al._ [1999] C. G. Durfee, A. R. Rundquist, S. Backus, C. Herne, M. M. Murnane, and H. C. Kapteyn, Phase matching of high-order harmonics in hollow waveguides, Phys. Rev. Lett. 83, 2187 (1999).
  * Skupin  _et al._ [2006] S. Skupin, G. Stibenz, L. Bergé, F. Lederer, T. Sokollik, M. Schnürer, N. Zhavoronkov, and G. Steinmeyer, Self-compression by femtosecond pulse filamentation: Experiments versus numerical simulations, Phys. Rev. E 74, 056604 (2006).
  * Jargot  _et al._ [2018] G. Jargot, N. Daher, L. Lavenu, X. Delen, M. H. N. Forget, and P. Georges, Self-compression in a multipass cell, Opt. Lett. 43, 5643 (2018).
  * Sansone  _et al._ [2006] G. Sansone, E. Benedetti, F. Calegari, C. Vozzi, L. Avaldi, R. Flammini, L. Poletto, P. Villoresi, C. Altucci, R. Velotta, S. Stagira, S. De Silvestri, and M. Nisoli, Isolated single-cycle attosecond pulses, Science 314, 443 (2006).
  * Ma  _et al._ [2016] X. Ma, J. Dostál, and T. Brixner, Broadband 7-fs diffractive-optic-based 2D electronic spectroscopy using hollow-core fiber compression, Opt. Express 18, 20781 (2016).
  * Nagy  _et al._ [2008] T. Nagy, M. Forster, and P. Simon, Flexible hollow fiber for pulse compressors, Appl. Opt. 47, 3264 (2008).
  * Nagy  _et al._ [2011] T. Nagy, V. Pervak, and P. Simon, Optimal pulse compression in long hollow fibers, Opt. Lett. 36, 4422 (2011).
  * Nagy  _et al._ [2019] T. Nagy, S. Hädrich, P. Simon, A. Blumenstein, N. Walther, R. Klas, J. Buldt, H. Stark, S. Breitkopf, P. Jójárt, I. Seres, Z. Várallyay, T. Eidam, and J. Limpert, Generation of three-cycle multi-millijoule laser pulses at 318 W average power, Optica 6, 1423 (2019).
  * Nagy  _et al._ [2020] T. Nagy, M. Kretschmar, M. J. J. Vrakking, and A. Rouzée, Generation of above-terawatt 1.5-cycle visible pulses at 1 kHz by post-compression in a hollow fiber, Opt. Lett. 45, 3313 (2020).
  * Tempea and Brabec [1998] G. Tempea and T. Brabec, Theory of self-focusing in hollow waveguide, Opt. Lett. 23, 762 (1998).
  * Fibich and Gaeta [2000] G. Fibich and A. Gaeta, Critical power for self-focusing in bulk media and hollow waveguides, Opt. Lett. 25, 335 (2000).
  * Nurhuda  _et al._ [2003] M. Nurhuda, A. Suda, K. Midorikawa, M. Hatayama, and K. Nagasaka, Propagation dynamics of femtosecond laser pulses in a hollow fiber filled with argon: constant gas pressure versus differential gas pressure, J. Opt. Soc. Am. B 20, 2002 (2003).
  * Andreasen and Kolesik [2013] J. Andreasen and M. Kolesik, Midinfrared femtosecond laser pulse filamentation in hollow waveguides: A comparison of simulation methods, Phys. Rev. E 87, 053303 (2013).
  * Safaei  _et al._ [2020] R. Safaei, G. Fan, O. Kwon, K. Légaré, P. Lassonde, B. E. Schmidt, H. Ibrahim, and F. Légaré, High-energy multidimensional solitary states in hollow-core fibres, Nature Photonics 14, 732 (2020).
  * Crego  _et al._ [2019] A. Crego, E. C. Jarque, and J. S. Roman, Influence of the spatial confinement on the self-focusing of ultrashort pulses in hollow-core fibers, Sci. Rep. 9, 9546 (2019).
  * López-Zubieta  _et al._ [2018] B. A. López-Zubieta, E. C. Jarque, Í. J. Sola, and J. S. Roman, Spatiotemporal-dressed optical solitons in hollow-core capillaries, OSA Continuum 1, 930 (2018).
  * Haus [2000] H. A. Haus, Mode-locking of lasers, IEEE J. Sel. Top. Quantum Electron. 6, 1173 (2000).
  * Escoto  _et al._ [2021] E. Escoto, A. Demircan, and G. Steinmeyer, Cage solitons, IEEE J. Quantum Electron. 57, 1300106 (2021).
  * Wise [2010] F. W. Wise, Generation of light bullets, Physics 3, 107 (2010).
  * Renninger and Wise [2013] W. H. Renninger and F. W. Wise, Optical solitons in graded-index multimode fibres, Nature Communications 4, 1719 (2013).
  * Lopez-Galmiche  _et al._ [2016] G. Lopez-Galmiche, Z. S. Eznaveh, M. A. Eftekhar, J. A. Lopez, L. G. Wright, F. Wise, D. Christodoulides, and R. A. Correa, Visible supercontinuum generation in a graded index multimode fiber pumped at 1064 nm, Opt. Lett. 41, 2553 (2016).
  * Crenn [1984] J. P. Crenn, Optical study of the EH11 mode in a hollow circular oversized waveguide and Gaussian approximation of the far-field pattern, Appl. Opt. 23, 3428 (1984).
  * Goos and Hänchen [1947] F. Goos and H. Hänchen, Ein neuer und fundamentaler Versuch zur Totalreflexion, Annalen der Physik 436, 333 (1947).
  * Crenn [1993] J. P. Crenn, Optical propagation of the HE11 mode and Gaussian beams in hollow circular waveguides, Int. J. Infrared Millimeter Waves 14, 1947 (1993).
  * Snyder and Love [1983] A. W. Snyder and J. D. Love,  _Optical Waveguide Theory_ (Chapman and Hall, London, UK, 1983).
  * Agrawal [2019] G. P. Agrawal,  _Nonlinear Fiber Optics_ , 6th ed. (Academic Press, London, UK, 2019).
  * Vinetskiǐ  _et al._ [1979] V. L. Vinetskiǐ, N. V. Kukhtarev, S. G. Odulov, and M. S. Soskin, Dynamic self-diffraction of coherent light beams, Sov. Phys. Usp. 22, 742 (1979).
  * DeLong  _et al._ [1995] K. W. DeLong, R. Trebino, and D. Kane, Comparison of ultrashort-pulse frequency-resolved-optical-gating traces for three common beam geometries, J. Opt. Soc. Am. B 11, 1595 (1995).
  * Chapman  _et al._ [2010] R. T. Chapman, T. J. Butcher, P. Horak, F. Poletti, J. G. Frey, and W. S. Brocklesby, Modal effects on pump-pulse propagation in an Ar-filled capillary, Opt. Express 18, 13279 (2010).
  * Blow  _et al._ [1992] K. J. Blow, N. J. Doran, and S. J. D. Phoenix, The soliton phase, Opt. Commun. 88, 137 (1992).