# 3D Technologies for Large Area Trackers

Whitepaper Submitted to Snowmass 2013

G. Deptucha, U. Heintzb, M. Johnsona, C. Kenneyd, R. Liptona, M. Narianb, S.
Parkere, A. Shenaia, L. Spiegela, J. Thomc, and Z. Yef
aFermilab, P.O. Box 500, Batavia, Il , USA, bBrown University, Providence, RI
, USA, cCornell University, Ithaca, NY , USA, d SLAC National Accelerator
Laboratory, Menlo Park, CA, USA, e University of Hawaii, Honolulu, HI, USA f
University of Ilinois, Chicago, Chicago, Il, USA

Abstract
We describe technologies which can be developed to produce large area, low
cost pixelated tracking detectors. These utilize wafer-scale 3D electronics
and sensor technologies currently being developed in industry. This can result
in fully active sensor/readout chip tiles which can be assembled into large
area arrays with good yield and minimal dead area. The ability to connect
though the bulk of the device can also provide better electrical performance
and lower mass.

##  1 Introduction

In the next generation of collider experiments detectors will be challenged by
unprecedented luminosity and data volumes. We will need to build large area
arrays (100’s of meter2) of highly pixelated detectors with minimal dead area
and at reasonable cost. Current fine pitch bump bonding technologies require
individual placement of chips on sensors followed by a solder melt cycle.
Technologies which provide wafer-scale interconnect offer prospects of lower
cost, finer pitch and lower mass interconnects between sensors and readout
chips. In addition 3D interconnect technology allows connections through the
bulk of the silicon, providing low inductance and capacitance paths for signal
and power connections. Table 1 summarizes current and projected costs and
yields of various bonding technologies.

In the past pixel sensor arrays have been limited in the module area by:

  1. 1.

Space needed at the edges of the detectors to reduce the field near the
damaged dicing cut regions. These regions act as charge emitters and can cause
unacceptably large currents in the edge strips.

  2. 2.

Die size of the Readout Integrated Circuit (ROIC), determined by the die yield
and reticule area

  3. 3.

The need to provide edge locations for wire bonds for connection to the ROIC.

The development of active edge sensors by C. Kenney, S. Parker and
collaborators has addressed the first problem [1] [2] [3]. Recent work on the
large area FEI4 chip [4] has shown that a large pixel chip can be produced
with good yield. Processes related to 3D electronics can solve the third
problem. The combination of all three can result in a wafer scale fabrication
of tiles or arrays which can be integrated into large modules with high yield
and relatively low cost. Using such tiles, large area pixelated modules can be
assembled with known good integrated sensor/readout die with large pitch
backside bump bond interconnects.

Table 1: Current and projected costs and yields for sensor/readout integration technologies. Component | Current or projected cost | Yield | Comment
---|---|---|---
Readout IC | $8/cm${}^{2}\cite[cite]{[\@@bibref{}{3DCost}{}{}]}$ | 65-70% [6] | Current 3D wafer cost and
|  |  | FEI4 prototype yield
Active Edge Sensors | $53/cm2 |  $\approx$90% | Current cost for prototype
|  |  | 150 mm wafers
Silicon Strip Sensors | $10/cm2 |  $\approx$100% | CMS tracker costs
Bump bonding | $213/cm2 | 98%[7] | CMS forward pixel costs
(2007) |  |  | Yield $\equiv<$20 bad bumps/chip
Bump bonding | $62/cm2 | - | CMS forward pixel upgrade
(2012) |  |  |
Wafer scale | $0.04/cm2 |  $\approx$90% | Projected by Yole Development [8]
DBI bonding |  |  | for high volume production
Target Costs (2020s) | $10/cm2 | 80% | Assuming 200 mm sensor wafers
|  |  | and batch active edge process

##  2 Technologies

###  2.1 3D Circuits

3D circuitry (3DIC) is the generic term for a set of technologies, including
wafer bonding, thinning, and interconnect, which allows vertical
interconnection of multiple layers of CMOS electronics [9] [10]. 3D
interconnects have the advantage of reducing inductance and capacitance while
increasing circuit density and allowing the integration of heterogeneous
device types. The key enabling technology is the Through-Silicon-Via (TSV), a
metal-filled hole etched into the wafer bulk silicon which forms the
conducting path between tiers of a multi-layer assembly. In high energy
physics 3D circuitry would allow us to directly integrate sensors and their
readout electronics without the use of expensive and cumbersome fine pitch
bump bonds.

We have explored three technologies for 3D devices. Our initial studies were
with MIT-Lincoln Labs and used their 0.18 micron three-tier process. This
process utilizes oxide bonding to join the tiers and vias are etched through
the SOI buried oxide after the tiers are bonded. We developed a demonstration
3-tier ILC vertex chip (VIP) in the MIT-LL process [11]. We have used the
Ziptronix oxide bonding process with imbedded metal to mate BTeV FPiX ROIC
wafers to sensors fabricated at MIT-Lincoln Labs [12]. Fermilab has also
sponsored a two-tier 0.13 micron CMOS 3D IC run with Tezzaron/Global Foundries
[13] [14] that features $1.2~{}\mu$ diameter, $6\mu$ deep tungsten filled
TSVs. In this process wafers were bonded face-to-face utilizing either copper
thermocompression or the Ziptronix Direct Bond Interconnect (DBI) oxide
bonding process described below. Bond pitch for these wafers was 4 microns.
After bonding the top wafer is thinned to expose the tungsten through-silicon-
vias and the top is patterned to provide contacts for bonding. The final set
of wafers from Tezzaron have been received. The copper-copper wafers suffered
from alignment problems, lowering the yield of good chips, but the DBI wafers
have both good alignment and yield.

Other 3D processes are becoming commercially available. IBM and Micron
Semiconductor will soon announce a 3D stacked memory product based on IBM’s 32
nm TSV process [15] with copper stud bonding shown in figure 1. This
technology is based on die-to-die bonding, potentially maximizing yield by
utilizing known good die from both sensors and ROICS. Combining this process
with active edge sensors could provide a very appealing path toward active
tiles which could be used to build large area sensors. CEA-LETI is offering
”Open 3D” services for prototyping and/or low volume production. Via-last
processes, such as provided by LETI, can provide TSVs in processed wafers,
with the limitation that the wafers must be thick enough to handle (¿100
microns) and the TSV aspect ratio is generally limited to less than 20:1.

Figure 1: Cross section of a 3D assembly utilizing the IBM 32 nm copper stud
3D interconnect process.

###  2.2 Active Edge Sensors

Active edge sensors are an outgrowth of work done to develop 3D sensors, which
provide good charge collection combined with radiation hardness. The technique
utilizes a deep reactive ion etch of silicon to create a nearly vertical
trench with smooth edges. The high quality of the trench wall avoids charge
generation normally associated with saw-cut edges and allow closer placement
of adjacent sensors. The trenches are filled with doped polycrystalline
silicon. The dopant is diffused into the surrounding single crystal silicon
and annealed for activation. The dopant density gradient will make an electric
field in the collection direction. These steps may be done at the same time
for the other like-type electrodes. The depth of the doped silicon must be
great enough so it is not depleted by the largest applied bias voltage.
Mechanical integrity is maintained by bonding the sensor to a support wafer.
The oxide bond also forms an etch stop for the trenching and singulation
processes. In the case of an oxide bonded handle, in a silicon-on-insulator
structure, the bond also forms an etch stop for subsequent trenching and
singulation processes. An alternative active edge technology based on wafer
cleaving and atomic layer deposition [16] has the prospect of achieving
similar goals without the additional processing needed in the deep trench
process.

###  2.3 Oxide Bonding

Bonding of silicon wafers is a key enabling technology for nanotechnology,
micromachining, and 3D electronics. A variety of techniques have been
developed, including adhesive bonding, metal eutectic bonding, and bonding
based on the silicon oxide surface either grown or deposited on a wafer [17].
Oxide bonding has the advantages of being mechanically robust, chemically
inert, and capable of withstanding the high temperatures typical of silicon
processing.

The direct oxide bond [18] is formed by bringing together silicon wafers which
have been planarized and chemically treated to form a hydrophilic surface.
When the wafers are brought together at room temperature a van der Walls bond
forms between the hydrogen atoms at the wafer surface. Further annealing above
$150^{\circ}{\rm C}$ causes the formation of covalent hydrogen bonds and
provides a substantial increase in bond strength. The Ziptronix DBI process
imbeds nickel or copper in the planarized oxide surface. The metal forms an
interconnect to a seed metal layer in the resulting oxide bonded wafers which
can be used to build 3D interconnect structures [19]. The process requires
good planarity, which is a feature of modern CMOS processes. Dust particles
present during the bonding process can produce local bond voids. These
unbonded areas limit the large area module yield and are a motivation for
active tile development. The DBI process can be used for wafer-to-wafer or
chip-to-wafer bonding. However the wafer to wafer process is least expensive
and most well suited to scaling to large volumes.

###  2.4 Interconnect

Once bonded to a sensor, connections need to be made to the ROIC, which is now
face down. Normally in a bump bonded pixel detector assembly the chip area is
larger than that of the sensor with bond pads extending out of the edges. In
an oxide bonded assembly signals as well as power need to be brought in
through the bulk silicon of the ROIC. If the ROIC has imbedded through silicon
vias the wafers can be thinned and vias exposed using the same technique
utilized in the Tezzaron two-tier 3D integrated circuit run. A non-TSV wafer
can be used by thinning the top silicon to $\approx 10\mu$, etching that
silicon to the I/O pad contacts, and depositing a redistribution metal layer
to provide the final contacts through what normally would be the bottom of the
ROIC. The redistribution layer can provide interconnects between adjacent
chips, route power, and bump bond pads. This layer can also provide
”stitching” between adjacent reticules to provide a large effective area which
can be defined at the dicing stage.

Figure 2: Sketch of an interposer-based large area array.

An alternate interconnect approach is to utilize so-called 2.5D interconnect.
This technology utilizes a silicon or glass interposer which incorporates TSVs
and redistribution of signals (figure 2). This approach decouples the 3D
processing from the sensors and readout chips and offers a substrate with very
fine interconnect capability and coefficient of thermal expansion match to the
silicon ROICs and sensors. An interposer-based solution could partially
decouple sensor and ROIC pitches, allowing for a variety of external
interconnects, and allow for conventional bump bonding of ROICS.

##  3 Large Area Assembly

We now consider techniques for the fabrication of large area arrays utilizing
the tools described in section 2. All are based on wafer-to-wafer bonding of
sensors and ROICs utilizing the DBI process.

Figure 3: Top - 3D sketch showing a view of the wafer stack, including handle
wafer, trenches, and the region between reticules. Bottom- Schematic view of
the final stack with DBI contact layers and sensor and top contacts.

###  3.1 Tiles with Active Edge Sensors

We are currently exploring a process which utilizes ”standard” active edge
sensors. These are n-on-p devices with 200 micron thick sensors in an SOI
stack with a 500 micron thick handle wafer. The tiles are fabricated by DBI
bonding the sensor wafer stack to a dummy ROIC which rearranges the signal
contact pattern. After bonding contacts on the dummy ROIC wafer are exposed by
grinding and etching away the dummy wafer silicon substrate and first layer of
oxide to expose what would normally be the bottom surface of the readout pads.
This surface is then metalized to provide sites for bump bond placement. The
stack is shown in figure 3. Finally the active die must be singulated. This is
accomplished by etching the 10 microns of dummy wafer silicon and 3 microns of
oxide, followed by an etch of the 200 microns of polysilicon filling the
trench. The SOI oxide forms an etch stop in this process. The wafer is then
attached to a temporary handle wafer and the original SOI handle is ground
away, leaving isolated tiles.

###  3.2 Post-Processed Tiles

A simpler variant on the above process would DBI bond ROIC and untrenched
sensor wafers. The tiles are then defined by etching trenches in the backside
of the sensor wafer. The etching process is low temperature and should not
harm the ROIC. However the resulting edges would not have the doping that
normally defines the side electrode in active edge devices. Post-trench doping
and annealing would not be possible, since the anneal would require higher
temperatures than could be tolerated by the ROIC. In this case a process
similar to that explored for slim-edge fabrication, utilizing atomic layer
deposition to activate the edges should provide acceptable performance. This
process avoids the SOI stack needed for the initial sensor wafer, provides a
more planar surface for oxide bonding, and avoids the complex singulation and
handle wafer removal necessary with the pre-processed active edge devices.

This process is very similar to plasma dicing, which utilizes the same deep
reactive ion etching process to dice wafers [20], so it should be commercially
available for both 200 mm and 300 mm diameter wafers.

###  3.3 Die-to-Die tiles

An alternative which utilizes the die-to-die capability of an IBM-style
process could bond active edge sensor tiles to ROICs with imbedded TSVs. This
probably would not require ROIC and sensor wafers of equal diameter assuming
the copper bond preparation process was available for the smaller diameter
sensor wafers. Although this process might not have the economies of scale of
a wafer level process, improved yields by utilizing tested die might make this
process competitive in cost.

###  3.4 Interposer-based Array

If the density of voids in the DBI bonded assembly is sufficiently small an
appealing assembly would utilize a silicon interposer DBI bonded to the sensor
wafer. Such an assembly could, in principle, utilize the largest sensor that
can be inscribed in the wafer with good yield, typically about 1 cm inset from
the wafer edge. The interposers redistributes the signals to ROIC bump bonds
on the top of the interposer. This works well in a situation where the area of
the ROIC is smaller that that of the associated pixels.

##  4 Prospects

The ability to build large areas of pixelated arrays with complex readout
electronics will be crucial to future detectors. Bonding costs and yield will
define the limitations of these detectors. We have described a plan to
demonstrate pixelated tiles or large area arrays that have the prospect of
making such large area devices affordable. To complete this development sensor
wafers which match the 8” diameter of the ROIC wafers are needed. Such wafers
have been demonstrated in the SOI process by both Lapis/OK [21] and American
Semiconductor. We expect that commercial manufacturers will soon establish 8”
high resistivity sensor production.

## References

  * [1] C. J. Kenney, J. D. Segal, E. Westbrook, S. Parker, J. Hasi, C. Da Via, S. Watts and J. Morse, _Active-edge planar radiation sensors_ , Nucl. Instrum. Meth. A 565, 272 (2006).
  * [2] C. J. Kenney, S. I. Parker, and E. Walckiers, _Results from 3D sensors with wall electrodes: near-cell-edge sensitivity measurements as a preview of active-edge sensors_ , IEEE Trans. Nucl. Sci. 48 2405-2410 (2001) .
  * [3] Cinzia Da Via, Mario Deile, Jasmine Hasi, Christopher Kenney, Angela Kok, Sherwood Parker, Stephen Watts, et al., _3D Active Edge Silicon Detector Tests with 120 GeV Muons_ , IEEE Trans. Nucl. Sci. 56 505 (2009).
  * [4] M.Garcia-Sciveres,et. al.,_The FE-I4 pixel readout integrated circuit_ , Nucl. Instrum. Meth. A 636, 155 (2011).
  * [5] Current cost for 0.18 micron CMOS wafers with through-siicon vias.
  * [6] M. Barbero, Vertex 2011 Workshop, Rust - Austria, June 19th - 24th 2011.
  * [7] P. Merkel et. al. Nucl. Instrum. Meth. A 582 771 (2007).
  * [8] Solid State Technology. http://www.electroiq.com/articles/sst/2008/10/ziptronix-joins-low-cost-quest-for-true-3d-ic.html
  * [9] IBM Journal of Research and Development, Volume 52, No. 6, 2008. Issue devoted to 3D technology.
  * [10] Philip Garrou, Christopher Bower, Peter Ramm, _Handbook of 3D Integration Technology and Applications of 3D Integrated Circuits_ , Wiley-VCH, 2008.
  * [11] G. Deptuch, M. Demarteau, J. Hoff, F. Khalid, R. Lipton, A. Shenai, M. Trimpl and R. Yarema et al., _Vertically integrated pixel readout chip for high energy physics_ , FERMILAB-CONF-11-020-PPD.
  * [12] Z. Ye, in Proceedings of the International Linear Collider Workshop (LCWS08 And ILC08) 2008, edited by M. Barone, Y. Torun, N. Varelas, arXiv:0902.2801v1 .
  * [13] G. Deptuch, M. Demarteau, J. Hoff, R. Lipton, A. Shenai, M. Trimpl, R. Yarema and T. Zimmerman, _Vertically Integrated Circuits at Fermilab_ , IEEE Trans. Nucl. Sci. 57, 2178 (2010).
  * [14] R. Yarema, _The First Multiproject Wafer Run With Chartered/Tezzaron_ , 2011 Front End Electronics Meeting Bergamo Italy, https://indico.cern.ch/getFile.py/access?contribId=34&sessionId=1&resId=1&materialId=slides&confId=122027
  * [15] M. Shapiro, AIDA-Academia Meets Industry, https://indico.cern.ch/contributionDisplay.py?sessionId=5&contribId=38&confId=209454
  * [16] R Mori et. al., _Charge collection measurements on slim-edge microstrip detectors_ , JINST 7 P05002 2012.
  * [17] Q.-Y. Tong and U. G sele _The Electrochemical Society. ed. Semiconductor Wafer Bonding: Science and Technology (1 ed.)_. Wiley-Interscience. ISBN 978-0-471-57481-1 1998.
  * [18] Q.-Y. Tong and E. Schmidt and U. G sele and M. Reiche. _Hydrophobic silicon wafer bonding._ Applied Physics Letters 64 (5): pp. 625-627 1994.
  * [19] P. Enquist, G. Fountain, C. Petteway, A. Hollingsworth, and H. Grady. _Low Cost of Ownership scalable copper Direct Bond Interconnect 3D IC technology for three dimensional integrated circuit applications._ In: 3D System Integration, 3DIC 2009;
  * [20] Wei-Sheng Lei, Ajay Kumar, and Rao Yalamanchili, _Die singulation technologies for advanced packaging: A critical review_ J. Vac. Sci. Technol. B 30, 040801 (2012)
  * [21] Y. Onuki, H. Katsurayama, Y. Ono, H. Yamamoto, Y. Arai, Y. Fujita, R. Ichiyama and Y. Ikegami et al., _SOI detector developments_ , PoS VERTEX 2011, 043 (2011).


††thanks: Present address: Institute of Liberal Arts and Sciences, Nagoya
University, Nagoya, Aichi 464-8601, Japan††thanks: Present address: Center for
General Education, Aichi Institute of Technology, Toyota, Aichi 470-0392,
Japan

#  Current-Induced Metallization and Valence Transition in Black SmS

Shin-ichi Kimura<EMAIL_ADDRESS>Graduate School of
Frontier Biosciences, Osaka University, Suita, Osaka 565-0871, Japan
Department of Physics, Graduate School of Science, Osaka University, Toyonaka,
Osaka 560-0043, Japan  Institute for Molecular Science, Okazaki, Aichi
444-8585, Japan  Hiroshi Watanabe  Graduate School of Frontier Biosciences,
Osaka University, Suita, Osaka 565-0871, Japan  Department of Physics,
Graduate School of Science, Osaka University, Toyonaka, Osaka 560-0043, Japan
Shingo Tatsukawa  Department of Physics, Graduate School of Science, Osaka
University, Toyonaka, Osaka 560-0043, Japan  Takuto Nakamura  Graduate School
of Frontier Biosciences, Osaka University, Suita, Osaka 565-0871, Japan
Department of Physics, Graduate School of Science, Osaka University, Toyonaka,
Osaka 560-0043, Japan  Keiichiro Imura  Department of Physics, Graduate School
of Science, Nagoya University, Nagoya Aichi 464-8602, Japan  Hiroyuki S.
Suzuki  Institute for Solid State Physics, The University of Tokyo, Kashiwa,
Chiba 277-8581, Japan  Noriaki K. Sato  Department of Physics, Graduate School
of Science, Nagoya University, Nagoya Aichi 464-8602, Japan

###### Abstract

A strongly-correlated insulator, samarium mono-sulfide (SmS), presents not
only the pressure-induced insulator-to-metal transition (IMT) with the color
change from black to golden-yellow but also current-induced IMT (CIMT) with
negative resistance. To clarify the origin of the CIMT of SmS, the electronic
structure change has been investigated by optical reflectivity and angle-
integrated photoelectron spectra by applying an electric current. At lower
temperatures than about 100 K, where the nonlinear $V$-$I$ curve has been
observed, the carrier density rapidly increases, accompanied by decreasing
relaxation time of carriers with increasing current. Then, the direct gap size
increases, and the mean valence changes from Sm2+-dominant SmS to the mixed-
valent one with increasing current. These results suggest that the CIMT
originates from increasing the Sm $4f$-$5d$ hybridization intensity induced by
the applied current.

Insulator-to-metal transition (IMT) of solids has been one of the critical
topics of interest for a long time because of the drastic change in electrical
resistivity that is useful for microscopic switching [1, 2]. IMT has been
observed in many transition-metal compounds and organic conductors, in which
they accompany a Jahn-Teller effect [3], a charge/orbital ordering [4], and a
lattice distortion due to the creation of charge/spin-density waves [5]. On
the other hand, in rare-earth compounds, for instance, in SmB6 and YbB12 and
others, metallic property at higher temperatures changes to a Kondo insulator
at low temperatures with a tiny energy gap due to the hybridization between
conduction and $f$ states ($c$-$f$ hybridization) by the Kondo effect [6]. The
nonmagnetic insulating property of EuO at high temperatures changes to a
ferromagnetic metal due to the indirect exchange interaction [7, 8]. These
rare-earth compounds do not show a lattice distortion because of the weak
crystal electric field of $4f$ electrons. Samarium mono-sulfide (SmS), the
title compound, also shows IMT by applying pressure without lattice distortion
as other rare-earth compounds.

SmS with NaCl-type crystal structure has the pressure-induced first-order
phase transition from the black-colored semiconductor ($b$-SmS) to the golden-
yellow-colored semimetal ($g$-SmS) (black-to-golden phase transition: BGT) [9,
10, 11, 12, 13, 14]. BGT can also be produced by applying chemical pressure of
the substituting trivalent yttrium (Y) ions with a smaller ionic radius to
Sm3+ ions [15, 16, 17, 18]. The BGT accompanies the valence transition from
Sm2+ to Sm3+, and the lattice constant decreases by about 5 % due to the
shrinking of the ionic radius without change in the crystal structure [19,
20]. At BGT, the electrons’ contribution of Sm ions changes from Sm2+ $4f^{6}$
in $b$-SmS to Sm3+ $4f^{5}+c$ ($c$: conduction electron) in $g$-SmS by the
BGT, so the appearance of the carriers is considered the origin of the IMT.
Since the valence change and the shrinking of the lattice constant
coincidently occur, however, more than 50 years after its discovery, the
origin of BGT of this material is still under debate [21]. One of the proposed
ideas for the origin of BGT is the BEC–BCS transition of Sm $4f$-$5d$ exciton
[22]. The origin of the BEC–BCS transition is the spontaneous emergence of
electron-hole pairs, namely exciton, due to the decrease of the energy gap by
applying pressure. If so, creating the $4f$-$5d$ exciton due to other methods
is expected to make BGT.

Recently, Ando et al. reported that the voltage ($V$) as a function of current
($I$) changes non-linearly from the high $dV/dI$ at low $I$ to the low $dV/dI$
at high $I$ via a negative $dV/dI$ at lower temperatures than 100 K, namely a
current-induced IMT (CIMT). It is discussed to be related to the BEC–BCS
transition of the $4f$-$5d$ excitons [23]. The change of the $dV/dI$ value
reflects the difference in the parameters of the conduction electrons, the
carrier density ($N$), and/or the mobility ($\mu$), which is related to the
electronic structure change. However, the origin of the CIMT has yet to be
revealed.

In this Letter, to investigate the origin of the CIMT of $b$-SmS from the view
of the electronic structure, we measured reflectivity ($R(\omega)$) spectra in
the broad energy region from the terahertz (THz) to the vacuum-ultraviolet and
angle-integrated photoelectron (PE) spectra by applying electric currents. As
a result, the carrier density rapidly increases with decreasing mobility by
applying current in the negative resistance region. In addition, the energy
gap size increases, and the mean valence of Sm-ions increases with increasing
current. However, the carrier density and the mean valence after CIMT do not
reach those of $g$-SmS, suggesting a different mechanism of BGT or a precursor
of BGT. The obtained results suggest that the localized Sm $4f$ electrons
change to an itinerant character by applying current, similar to melting ice
by flowing water, and the itinerant $4f$ electrons hybridize with the Sm $5d$
conduction band. The observed CIMT of $b$-SmS is concluded as a novel type of
valence transition.

Single crystalline samples of SmS were grown using the vertical Bridgman
method with a high-frequency induction furnace: Pre-reacted materials starting
99.99 % pure (4N) samarium and 6N powdered sulfur were sealed in a vacuum W
crucible [14]. Cleaved samples with the size of about $1\times 1\times 1$ mm3
were mounted on the sample folder with the current flow and cooled down to
about 20 K for the $R(\omega)$ measurement and to 30 K for the PE measurement.
The direct current up to 1.6 A was applied while maintaining a temperature
range of up to 5 K to avoid temperature rise due to Joule heating and to
prevent sample degradation. Near-normal-incident $R(\omega)$ spectra along the
$(001)$ plane were acquired in a wide photon-energy range of 6 meV–10 eV to
obtain optical conductivity ($\sigma_{1}(\omega)$) spectra via the Kramers-
Kronig analysis. Infrared (IR) and THz measurements at the photon energy
$\hbar\omega$ regions of 50 meV–1.5 eV and 6–50 meV have been performed using
IR and THz microscopes with synchrotron radiation at the beamline 6B of UVSOR-
III, Institute for Molecular Science, Japan [24, 25]. In the photon energy
range of 2–10 eV, the $R(\omega)$ spectra were measured using the synchrotron
radiation setup at the beamline 3B of UVSOR-III [26]. PE spectra have been
accumulated using a He-II discharge lamp with a monochromator (VG-Scienta VUV
5000+5040), a hemispherical analyzer (MBScientific A-1), and an own-developed
cold sample finger with a current-flow system. The base pressure for the PE
measurement was less than $1\times 10^{-8}$ Pa.

Figure 1:  (a) $I$-$V$ curves of $b$-SmS at temperatures from 20 to 250 K. (b)
Current-dependent reflectivity ($R(\omega)$) spectra of $b$-SmS in the THz
region at $T=40$ K. (c, d) Temperature and external-current dependence of (c)
the effective carrier density ($N_{eff}$) and (d) the relaxation time ($\tau$)
of $b$-SmS evaluated by the Drude-Lorentz (DL) fitting of the THz $R(\omega)$
spectra shown in (b) and at other temperatures. (e, f) Same as (c) and (d),
but all of the DL fitting data and schematic curved surfaces are also plotted
as guides for the eye.

Firstly, the temperature dependence of the $I$-$V$ behavior has been checked.
Figure 1a shows the $I$-$V$ characteristics of $b$-SmS measured with
increasing electric current at temperatures from 20 to 250 K. The results are
qualitatively consistent with that reported previously [23]. The $I$-$V$ curve
changes to a monotonic increase at temperatures above 80 K. Below 50 K, the
$I$-$V$ curves show a sizable non-linear behavior, and finally, below 40 K,
negative $dV/dI$ appears in the current from 30 to 120 mA. At 40 K, the
$dV/dI$ value of $\sim 10^{2}$ V/A at $I=20$ mA changes to $\sim 2.5$ V/A at 1
A, i.e., 50 times decreases, but at 80 K, it decreases about two times
smaller. This fact suggests that the conduction mechanism and the electronic
structure change between the low-current and high-current regions at
temperatures below 50 K. On the other hand, at higher temperatures than 80 K,
the electronic structure is unchanged. The result at higher than 80 K is
consistent with the recent electric-field-induced IMT results measured at
higher than 77 K, which can be explained by Joule heating effect [27].

To obtain information on carriers at different currents, we derived the
temperature and current dependences of the carrier density and the relaxation
time from the Drude-Lorentz (DL) fitting of the $R(\omega)$ spectra in the THz
region [28, 29]. As an example of the current-dependent THz $R(\omega)$
spectra, spectra at $T=40$ K are shown in Fig. 1b. At $I=0$ A, the spectrum is
a typical insulating one, with a phonon peak at $\sim$25 meV, and the
background intensity is low, suggesting a low carrier density. With increasing
current, the background intensity increases, accompanied by a rapid change in
the phonon shape originating from the interface between the phonon and
background absorptions. After analyzing the spectral change using the DL
fitting, we obtained the strong temperature- and current-dependent effective
carrier density ($N_{eff}$) and the relaxation time ($\tau$) as shown in Figs.
1c and 1d. To show the temperature and current dependences more clearly, the
three-dimensional plots are also shown in Figs. 1e and 1f. On the other hand,
obtained parameters for the phonon peak showed little change with applied
current as shown in Fig. S2 in the Supplementary Material [29], suggesting the
Joule heating effect is not large, consistent with the temperature kept within
5 K. Therefore, the current dependence of $N_{eff}$ and $\tau$ originates from
the applied current itself.

As shown in Fig. 1c, $N_{eff}$ increases with increasing current at all
temperatures. At higher temperatures than 150 K, the rate of increase
($dN_{eff}/dI$) is almost constant at different temperatures ($\sim 1.2\times
10^{16}~{}{\rm cm}^{-3}{\rm A}^{-1}$, dotted lines in Fig. 1c), but the
intercept at 0 A increases with increasing temperature. The $N_{eff}$ at 0 A
originates from the thermally excited carriers because the smallest (indirect)
energy gap size of the material is about 0.1 eV, comparable to the thermal
excitation energy of $k_{\rm B}T$. The similar $dN_{eff}/dI$ value is also
observed at higher currents than 0.6 A of lower temperatures than 100 K,
suggesting the electronic structure in the region is consistent with the high-
temperature range.

On the other hand, at lower currents than 0.6 A and lower temperatures than
100 K, the $dN_{eff}/dI$ value is about $3.2\times 10^{16}$ cm-3A-1 (the
dashed line in Fig. 1c), which is about three times larger than that at the
high-current region. $\tau$ also enhances in the same temperature and current
ranges. Usually, $\tau$, as well as the electron mobility $\mu$, have a
carrier-density dependence due to the electron-electron scattering, but, in
the low-carrier-density region at around $10^{16}$ cm-3, they do not usually
change so much [30, 31]. Therefore, the change of $\tau$ at low temperature
and low current region suggests that the character of carriers and the
electronic structure are modulated from the region of the low temperature and
low current to that of the high temperature and high current.

Figure 2:  Current-dependent optical conductivity($\sigma_{1}(\omega)$)
spectra of $b$-SmS in $\hbar\omega$ below 0.7 eV at $T=50$ K, derived from the
Kramers-Kronig analysis (KKA) of the reflectivity ($R(\omega)$) spectra. In
the low current of $I=0.1$ A, the onset of the reading edge is about 0.28 eV,
whereas it shifts to about 0.32 eV in the high current of 0.9 A. The inset
shows that the absorption edges in the absorption ($\alpha(\omega)$) spectra
obtained from the KKA obey the $\exp(\hbar\omega)$-law, suggesting evidence of
a direct gap from the Urbuch rule of exciton peaks [32]. A sharp peak at about
20 meV originates from a TO phonon that is the same as the peak shown in Fig.
1b.

To clarify the electronic structure change by the current application, the
current-dependent $\sigma_{1}(\omega)$ spectra at $T=50$ K are shown in Fig.
2. At the low-current region of $I=0.1$ A, a peak appears at about 0.56 eV
originating from the exciton absorption at the direct gap at the $X$ point
[28], as evidenced in the inset of Fig. 2. By applying a current of $\sim$0.9
A, the peak at $\sim 0.64$ eV shifts to the high-energy side. The onset of the
reading edge also increases by about 40 meV from about 0.28 eV to 0.32 eV with
increasing current. These results indicate the enlargement of the energy gap
size. The obtained $\sigma_{1}(\omega)$ spectrum at high current is completely
different from that of $g$-SmS, in which there is a Drude-like spectral shape
with a plasma edge at about 3 eV [33, 11, 12]. Therefore, the electronic
structure of the high-current phase is not the same as that of $g$-SmS.
Instead, the high-current spectrum in the low-energy region below 0.2 eV is
similar to the spectrum at the pressure of 0.6 GPa, just below the boundary of
BGT [28]. However, the exciton peak shifts to the low-energy side by a
pressure application, which trend is the opposite behavior of the current
application. The pressure-induced energy-gap narrowing below 0.6 GPa indicates
that the material is still in the black phase at the pressure; therefore, the
thermally excited carrier appears. However, in the high-current spectrum, the
energy gap size increases, even though the carrier density increases with
increasing current. This fact suggests that the carriers in the high-current
region do not originate from thermal excitation. One scenario for the energy-
gap enlargement at high currents is the appearance of the bonding and anti-
bonding states of the hybridization between the Sm $4f$ valence band and the
$5d$ conduction band. In the low-current region, a direct gap is expected at
the $X$ point in the Brillouin zone [34]. Then, since the Sm $4f$ states are
isolated, the Sm $5d$ conduction band, which orbital the carriers go through,
is not hybridized. At the high-current region, the isolated $4f$ electrons
change to itinerant; the $4f$ electrons can be hybridized with the carriers,
i.e., the hybridization bands, namely bonding and anti-bonding states of the
$4f$ states and $5d$ conduction bands, are created. The energy gap between the
bonding and anti-bonding states is enlarged with increasing the hybridization
intensity. In addition, the mean valence of Sm ions is expected to be changed
from $2+$ to $2+\delta$ due to the $4f$-$5d$ hybridization. Then, the carrier
density released from the Sm $4f$ states can increase with increasing the
hybridization intensity.

Figure 3:  Current-dependent angle-integrated photoelectron (PE) spectra in
the valence band energy region as a function of binding energy ($E_{B}$). The
PE spectra are accumulated at about 30 K using the He-II line ($h\nu=40.8$ eV)
as a light source. (Inset) Enlarge figure of the Sm $5p$ core level.

To check the hybridization between the Sm $4f$ and $5d$ states at high
currents, we measured the current-dependent valence band PE spectra as shown
in Fig. 3. The figure shows that the PE spectrum has strong current
dependence, i.e., the intensity of the Sm${}^{2+}~{}4f$ multiplet peaks in the
binding energy range below 3 eV decreases, and in contrast, the Sm3+ peak at
about 8 eV increases with increasing current. Similarly, as shown in the
inset, the intensity of the Sm${}^{2+}~{}5p$ peaks decreases, and that of the
Sm${}^{3+}~{}5p$ peaks increases with the current application. These results
strongly support that the applied current increases the Sm mean valence from
$2+$ to $2+\delta$. Therefore, the scenario of the appearance of the
hybridization state between Sm $4f$ and $5d$ at the high-current region is
plausible.

Figure 4:  The scenario of the current-induced metallization and valence
transition of $b$-SmS. (a) Real space description of the current-induced
metallization. Localized $4f$ electrons are induced to be itinerant by applied
currents, like melting ice by flowing water. (b) Reciprocal lattice
description of the current-induced metallization and valence transition. After
changing to the itinerant $4f$ electrons at high currents, they hybridize to
carriers, namely $4f$-$5d$ hybridization. Then the bottom of the conduction
band becomes flat. The relaxation time $\tau$, which corresponds to the
mobility ($\mu=e\tau/m$, where $e$ is the electron charge and $m$ the electron
mass), becomes shorter than that in the low-current case. After the appearance
of the hybridization, the Sm valence changes from $2+$ to $2+\delta$. Then the
Fermi level ($E_{\rm F}$) moves to the higher energy side.

To explain all of the experimental evidence, the lowering $\tau$, the energy
gap enlargement and the valence transition from $2+$ to $2+\delta$ by current
application, we propose a model as shown in Fig. 4. At the ground state
without current application, the Sm${}^{2+}~{}4f$ states are localized, and
the energy gap between the Sm $4f$ and $5d$ states opens [34]. Since the $5d$
conduction band is a free-electron-like band and the Sm ions are localized and
non-magnetic divalent, the carriers’ relaxation time $\tau$ is long because of
the lack of strong interactions. With applying current, the flowing electrons
push the localized $4f$ electrons out from their localized sites by an
electron-electron interaction because the energy level of the localized $4f$
electrons is very close to $E_{\rm F}$, and then the $4f$ electrons become
itinerant. The itinerant $4f$ electrons can hybridize with carriers in the Sm
$5d$ states, so the $4f$-$5d$ hybridization state appears. Carriers in the
hybridization band have a shorter relaxation time because of the strong
interaction.

Similar phenomena to the CIMT and valence transition on $b$-SmS have not been
observed, to the best of our knowledge, even though current-induced phenomena,
such as a Néel vector switching [35], phase transitions in cuprate [36], and a
magnetic thin film [37], has recently been reported. Also, electric-field-
induced metallization on Ca2RuO4 [38, 39] is similar to our observation, but
the electric field in our study is not so high. Our observation of the
current-induced metallization and valence transition completely differs from
such works. This phenomenon may be related to other anomalous phenomena, such
as pressure-induced BGT, alkali-metal-induced surface valence transition [40],
photo-induced band shift [41], and the BEC–BCS transition of Sm $4f$-$5d$
exciton; therefore, a comprehensive understanding of these phenomena will help
to elucidate the fruitful physical properties of SmS.

To summarize, we have investigated the electronic structure change of
insulating $b$-SmS by the current application. By analyzing the THz
reflectivity spectra at temperatures below 100 K, the carrier density rapidly
increased with current application up to 0.8 A and slowly after that,
accompanied by the decreasing relaxation time, suggesting the changing of the
carriers’ character and the electronic structure at high currents. The current
application also increased the direct energy gap size at the $X$ point in the
Brillouin zone and the increase of Sm ions’ mean valence, suggesting that the
applied current induces the Sm $4f$-$5d$ hybridization. The observed current-
induced metallization accompanying the valence transition in $b$-SmS is a
novel phenomenon of insulator-to-metal transitions.

We appreciate Prof. T. Ito for his fruitful discussion and UVSOR Synchrotron
staff members for their support during synchrotron radiation experiments. This
work was partly performed under the Use-of-UVSOR Synchrotron Facility Program
(Proposals Nos. 20-735, 21-677, 22IMS6015, 22IMS6029) of the Institute for
Molecular Science, National Institutes of Natural Sciences. This work was
partly supported by JSPS KAKENHI (Grant No. 20H04453).

## References

  * Mott [1968] N. F. Mott, Rev. Mod. Phys. 40, 677 (1968).
  * Mott [2004] N. Mott,  _Metal-Insulator Transitions_ (CRC Press, 2004).
  * Han  _et al._ [2000] J. E. Han, E. Koch, and O. Gunnarsson, Phys. Rev. Lett. 84, 1276 (2000).
  * Moritomo  _et al._ [1997] Y. Moritomo, H. Kuwahara, Y. Tomioka, and Y. Tokura, Phys. Rev. B 55, 7549 (1997).
  * Grüner [2018] G. Grüner,  _Density Waves in Solids_ (CRC Press, 2018).
  * Kasuya [1994] T. Kasuya, Europhys. Lett. 26, 277 (1994).
  * Lee and Liu [1984] V.-C. Lee and L. Liu, Phys. Rev. B 30, 2026 (1984).
  * Miyazaki  _et al._ [2009] H. Miyazaki, T. Ito, H. J. Im, S. Yagi, M. Kato, K. Soda, and S. Kimura, Phys. Rev. Lett. 102, 227203 (2009).
  * Jayaraman  _et al._ [1970] A. Jayaraman, V. Narayanamurti, E. Bucher, and R. G. Maines, Phys. Rev. Lett. 25, 1430 (1970).
  * Batlogg  _et al._ [1974] B. Batlogg, J. Schoenes, and P. Wachter, Phys. Lett. A 49, 13 (1974).
  * Batlogg  _et al._ [1976] B. Batlogg, E. Kaldis, A. Schlegel, and P. Wachter, Phys. Rev. B 14, 5503 (1976).
  * Travaglini and Wachter [1984] G. Travaglini and P. Wachter, Phys. Rev. B 30, 5877 (1984).
  * Matsubayashi  _et al._ [2007a] K. Matsubayashi, K. Imura, H. S. Suzuki, G. Chen, N. Mori, T. Nishioka, K. Deguchi, and N. K. Sato, J. Phys. Soc. Japan 76, 033602 (2007a).
  * Matsubayashi  _et al._ [2007b] K. Matsubayashi, K. Imura, H. S. Suzuki, T. Mizuno, S.-i. Kimura, T. Nishioka, K. Kodama, and N. K. Sato, J. Phys. Soc. Japan 76, 064601 (2007b).
  * Penney and Holtzberg [1975] T. Penney and F. Holtzberg, Phys. Rev. Lett. 34, 322 (1975).
  * Guntherodt and Holtzberg [1976] G. Guntherodt and F. Holtzberg, Solid State Commun. 18, 181 (1976).
  * Kaneko  _et al._ [2014] M. Kaneko, M. Saito, T. Ito, K. Imura, T. Hajiri, M. Matsunami, S. Kimura, H. S. Suzuki, and N. K. Sato, JPS Conf. Proc. 3, 011080 (2014).
  * Yokoyama  _et al._ [2019] Y. Yokoyama, H. Hasegawa, Y. Mizuno, D. Asai, Y. Okamoto, H. S. Suzuki, K. Takehana, Y. Imanaka, and K. Takenaka, Phys. Rev. B 100, 245143 (2019).
  * Jarrige  _et al._ [2013] I. Jarrige, H. Yamaoka, J.-P. Rueff, J.-F. Lin, M. Taguchi, N. Hiraoka, H. Ishii, K. D. Tsuei, K. Imura, T. Matsumura, A. Ochiai, H. S. Suzuki, and A. Kotani, Phys. Rev. B 87, 115107 (2013).
  * Imura  _et al._ [2020] K. Imura, Y. Suzuki, N. Kawade, K. Sakamoto, K. Deguchi, H. Yamaoka, Y. Yamamoto, J. Mizuki, N. Hiraoka, H. Ishii, H. S. Suzuki, and N. K. Sato, JPS Conf. Proc. 30, 011131 (2020).
  * Watanabe [2021] S. Watanabe, J. Phys. Soc. Japan 90, 023706 (2021), 2101.02484 .
  * Varma [1976] C. M. Varma, Rev. Mod. Phys. 48, 219 (1976).
  * Ando  _et al._ [2020] H. Ando, K. Sakamoto, K. Imura, K. Deguchi, H. S. Suzuki, and N. K. Sato, JPS Conf. Proc. 30, 011132 (2020).
  * Kimura  _et al._ [2006] S. Kimura, E. Nakamura, T. Nishi, Y. Sakurai, K. Hayashi, J. Yamazaki, and M. Katoh, Infrared Phys. Technol. 49, 147 (2006).
  * Kimura  _et al._ [2007] S.-i. Kimura, Y. Sakurai, E. Nakamura, and T. Mizuno, in  _AIP Conf. Proc._, Vol. 879 (AIP, 2007) pp. 595–598.
  * Fukui  _et al._ [2014] K. Fukui, R.-i. Ikematsu, Y. Imoto, M. Kitaura, K. Nakagawa, T. Ejima, E. Nakamura, M. Sakai, M. Hasumoto, and S.-i. Kimura, J. Synchrotron Radiat. 21, 452 (2014).
  * Kishida  _et al._ [2022] K. Kishida, J. Shibutani, T. Hatano, Y. Yokoyama, N. Katayama, Y. Okamoto, Y. Nakamura, H. Kishida, H. S. Suzuki, K. Takehana, Y. Imanaka, and K. Takenaka, Appl. Phys. Lett. 121, 121901 (2022).
  * Mizuno  _et al._ [2008] T. Mizuno, T. Iizuka, S.-i. Kimura, K. Matsubayashi, K. Imura, H. S. Suzuki, and N. K. Sato, J. Phys. Soc. Japan 77, 113704 (2008).
  * [29] See Supplemental Material at http://?????? for the Drude-Lorentz fitting method and the fitting data of phonons.
  * Masetti  _et al._ [1983] G. Masetti, M. Severi, and S. Solmi, IEEE Trans. Electron Devices 30, 764 (1983).
  * Reggiani  _et al._ [2000] S. Reggiani, M. Valdinoci, L. Colalongo, M. Rudan, and G. Baccarani, VLSI Des. 10, 467 (2000).
  * Urbach [1953] F. Urbach, Phys. Rev. 92, 1324 (1953).
  * Kirk  _et al._ [1972] J. L. Kirk, K. Vedam, V. Narayanamurti, A. Jayaraman, and E. Bucher, Phys. Rev. B 6, 3023 (1972).
  * Kimura  _et al._ [2008] S.-i. Kimura, T. Mizuno, K. Matsubayashi, K. Imura, H. S. Suzuki, and N. K. Sato, Phys. B Condens. Matter 403, 805 (2008).
  * Wu  _et al._ [2022] H. Wu, H. Zhang, B. Wang, F. Groß, C.-Y. Yang, G. Li, C. Guo, H. He, K. Wong, D. Wu, X. Han, C.-H. Lai, J. Gräfe, R. Cheng, and K. L. Wang, Nat. Commun. 13, 1629 (2022).
  * Naito  _et al._ [2021] T. Naito, T. Nishizaki, and N. Kobayashi, J. Appl. Phys. 129, 043901 (2021).
  * Kim  _et al._ [2008] D.-k. Kim, A. Domenicucci, and S. S. Iyer, J. Appl. Phys. 103, 073708 (2008).
  * Nakamura  _et al._ [2013] F. Nakamura, M. Sakaki, Y. Yamanaka, S. Tamaru, T. Suzuki, and Y. Maeno, Sci. Rep. 3, 2536 (2013).
  * Zhang  _et al._ [2019] J. Zhang, A. S. McLeod, Q. Han, X. Chen, H. A. Bechtel, Z. Yao, S. N. Gilbert Corder, T. Ciavatti, T. H. Tao, M. Aronson, G. L. Carr, M. C. Martin, C. Sow, S. Yonezawa, F. Nakamura, I. Terasaki, D. N. Basov, A. J. Millis, Y. Maeno, and M. Liu, Phys. Rev. X 9, 011032 (2019).
  * Nakamura  _et al._ [2023] T. Nakamura, T. Nakaya, Y. Ohtsubo, H. Sugihara, K. Tanaka, R. Yukawa, M. Kitamura, H. Kumigashira, K. Imura, H. S. Suzuki, N. K. Sato, and S.-i. Kimura, Phys. Rev. B 107, L041102 (2023).
  * Chen  _et al._ [2022] Y. Chen, T. Nakamura, H. Watanabe, T. Suzuki, Q. Ren, K. Liu, Y. Zhong, T. Kanai, J. Itatani, S. Shin, K. Okazaki, K. Imura, H. S. Suzuki, N. K. Sato, and S.-i. Kimura https://doi.org/10.48550/arXiv.2212.00983 (2022), arXiv:2212.00983 .


Started document

Finished title page

^^JFinished preliminary forms^^J^^J

Modified page settings

[1]
L. Ambrosio, N. Fusco, and D. Pallara.
Functions of Bounded Variation and Free Discontinuity Problems.
Oxford Science Publications. Clarendon Press, 2000.

[2]
N. Aronszajn.
A unique continuation theorem for solutions of elliptic partial
differential equations or inequalities of second order.
J. Math Pures Appl, 36(9):235–249, 1957.

[3]
P. H. Bérard.
Spectral Geometry: Direct and Inverse Problems.
Lecture Notes in Mathematics. Springer, Berlin, Heidelberg, 1986.

[4]
V. Bobkov and E. Parini.
On the higher Cheeger problem.
Journal of the London Mathematical Society, 97(3):575–600,

[5]
V. Bobkov and E. Parini.
On the Cheeger problem for rotationally invariant domains.
manuscripta mathematica, 166(3):503–522, Nov 2021.

[6]
S. Breen.
Uniform upper and lower bounds on the zeros of Bessel functions of
the first kind.
Journal of Mathematical Analysis and Applications,
196(1):1–17, 1995.

[7]
Y. D. Burago and V. A. Zalgaller.
Geometric Inequalities.
Grundlehren der mathematischen Wissenschaften. Springer, Berlin,
Heidelberg, 1988.

[8]
P. Buser.
On Cheeger's inequality $\lambda_1 \ge h^2/4$.
In R. Osserman and A. Weinstein, editors, Geometry of the
Laplace Operator, pages 29–77. American Mathematical Society, Providence,
Rhode Island, 1980.

[9]
P. Buser.
A note on the isoperimetric constant.
Annales scientifiques de l'École Normale Supérieure, Ser.
4, 15(2):213–230, 1982.

[10]
G. Buttazzo, G. Carlier, and M. Comte.
On the selection of maximal Cheeger sets.
Differential and Integral Equations, 20(9):991 – 1004, 2007.

[11]
A. Caboussat, R. Glowinski, and V. Pons.
An augmented Lagrangian approach to the numerical solution of a
non-smooth eigenvalue problem.
Journal of Numerical Mathematics, 17(1):3–26, 2009.

[12]
G. Carlier, M. Comte, and G. Peyré.
Approximation of maximal Cheeger sets by projection.
ESAIM: Mathematical Modelling and Numerical Analysis -
Modélisation Mathématique et Analyse Numérique, 43(1):139–150, 2009.

[13]
M. Caroccia and S. Ciani.
Dimensional lower bounds for contact surfaces of Cheeger sets.
Journal de Mathématiques Pures et Appliquées, 157:1–44,

[14]
V. Caselles, G. Facciolo, and E. Meinhardt.
Anisotropic Cheeger sets and applications.
SIAM Journal on Imaging Sciences, 2(4):1211–1254, 2009.

[15]
I. Chavel.
Eigenvalues in Riemannian Geometry.
Academic Press, 1984.

[16]
I. Chavel.
Riemannian Geometry: A Modern Introduction.
Cambridge Studies in Advanced Mathematics. Cambridge University
Press, second edition, 2006.

[17]
J. Cheeger.
A lower bound for the smallest eigenvalue of the Laplacian.
In Problems in Analysis: A Symposium in Honor of Salomon
Bochner, pages 195–199. Princeton University Press, 1970.

[18]
R. R. Coifman and S. Lafon.
Diffusion maps.
Applied and Computational Harmonic Analysis, 21(1):5–30, 2006.
Special Issue: Diffusion Maps and Wavelets.

[19]
Y. Colin de Verdière.
Sur la multiplicité de la première valeur propre non nulle du
Commentarii Mathematici Helvetici, 61:254–270, 1986.

[20]
R. Courant and D. Hilbert.
Methods of Mathematical Physics, volume 1.
Interscience Publishers, New York, 1953.

[21]
X. Dai, G. Wei, and Z. Zhang.
Neumann isoperimetric constant estimate for convex domains.
Proc. Amer. Math. Soc., 146:3509–3514, 2018.

[22]
A. Daneshgar, H. Hajiabolhassan, and R. Javadi.
On the isoperimetric spectrum of graphs and its approximations.
Journal of Combinatorial Theory, Series B, 100(4):390–412,

[23]
E. B. Davies.
Metastable states of symmetric Markov semigroups II.
Journal of the London Mathematical Society, s2-26(3):541–556,

[24]
N. de Ponti and A. Mondino.
Sharp Cheeger–Buser type inequalities in $\mathsf{RCD}(k,
\infty)$ spaces.
The Journal of Geometric Analysis, 31:2416–2438, 2021.

[25]
G. Froyland.
Dynamic isoperimetry and the geometry of Lagrangian coherent
Nonlinearity, 28(10):3587–3622, 2015.

[26]
G. Froyland and O. Junge.
Robust FEM-based extraction of finite-time coherent sets using
scattered, sparse, and incomplete trajectories.
SIAM Journal on Applied Dynamical Systems, 17(2):1891–1924,

[27]
G. Froyland and E. Kwok.
A dynamic Laplacian for identifying Lagrangian coherent
structures on weighted Riemannian manifolds.
Journal of Nonlinear Science, Jun 2017.

[28]
G. Froyland, C. P. Rock, and K. Sakellariou.
Sparse eigenbasis approximation: Multiple feature extraction across
spatiotemporal scales with application to coherent set identification.
Communications in Nonlinear Science and Numerical Simulation,
77:81–107, 2019.

[29]
K. Funano.
Eigenvalues of Laplacian and multi-way isoperimetric constants on
weighted Riemannian manifolds.
arXiv e-prints, page arXiv:1307.3919, Jul 2013.

[30]
D. Gilbarg and N. S. Trudinger.
Elliptic partial differential equations of second order.
Springer-Verlag Berlin Heidelberg, 1977.

[31]
D. S. Grebenkov and B.-T. Nguyen.
Geometrical structure of Laplacian eigenfunctions.
SIAM Review, 55(4):601–667, 2013.

[32]
A. Grigor'yan.
Isoperimetric inequalities and capacities on Riemannian manifolds.
In J. Rossmann, P. Takáč, and G. Wildenhain, editors, The Maz'ya Anniversary Collection, pages 139–153. Birkhäuser Basel,
Basel, 1999.

[33]
A. Grigor’yan.
Heat kernels on weighted manifolds and applications.
In J. Jorgenson and L. Walling, editors, The Ubiquitous Heat
Kernel, pages 93–192. American Mathematical Society, Providence, Rhode
Island, 2006.

[34]
A. Hassannezhad, G. Kokarev, and I. Polterovich.
Eigenvalue inequalities on Riemannian manifolds with a lower
Ricci curvature bound.
Journal of Spectral Theory, 6(4):807–835, December 2016.

[35]
S. Hofmann, M. Mitrea, and M. Taylor.
Singular integrals and elliptic boundary problems on regular
Semmes-Kenig-Toro domains.
International Mathematics Research Notices,
2010(14):2567–2865, 12 2009.

[36]
A. J. Izenman.
Introduction to manifold learning.
WIREs Computational Statistics, 4(5):439–446, 2012.

[37]
J. Jost.
Partial Differential Equations.
Graduate Texts in Mathematics, 214. Springer New York, New York,
second edition, 2007.

[38]
D. Karrasch and J. Keller.
A geometric heat-flow theory of Lagrangian coherent structures.
Journal of Nonlinear Science, 30:1849–1888, 2020.

[39]
D. Karrasch and N. Schilling.
A Lagrangian perspective on nonautonomous advection-diffusion
processes in the low-diffusivity limit.
arXiv e-prints, 2021.

[40]
T. C. Kwok, L. C. Lau, Y. T. Lee, S. Oveis Gharan, and L. Trevisan.
Improved Cheeger's inequality: Analysis of spectral partitioning
algorithms through higher order spectral gap.
In Proceedings of the Forty-Fifth Annual ACM Symposium on Theory
of Computing, STOC '13, page 11–20. Association for Computing Machinery,
New York, NY, USA, 2013.

[41]
O. Lablée.
Spectral Theory in Riemannian Geometry.
EMS Textbooks in Mathematics. European Mathematical Society, 2015.

[42]
T. Lachand-Robert and E. Oudet.
Minimizing within convex bodies using a convex hull method.
SIAM Journal on Optimization, 16(2):368–379, Jan. 2005.

[43]
M. Ledoux.
Spectral gap, logarithmic Sobolev constant, and geometric bounds.
Surveys in Differential Geometry, 9:219–240, 2004.

[44]
J. M. Lee.
Introduction to Smooth Manifolds.
Graduate Texts in Mathematics. Springer-Verlag New York, second
edition, 2012.

[45]
J. R. Lee, S. Oveis Gharan, and L. Trevisan.
Multiway spectral partitioning and higher-order Cheeger
J. ACM, 61(6):37:1–37:30, Dec. 2014.

[46]
G. P. Leonardi.
An overview on the Cheeger problem.
In A. Pratelli and G. Leugering, editors, New Trends in Shape
Optimization, pages 117–139. Springer International Publishing, Cham, 2015.

[47]
G. P. Leonardi, R. Neumayer, and G. Saracco.
The Cheeger constant of a Jordan domain without necks.
Calculus of Variations and Partial Differential Equations,
56(164), 2017.

[48]
G. P. Leonardi and A. Pratelli.
On the Cheeger sets in strips and non-convex domains.
Calculus of Variations and Partial Differential Equations,
55:15:1, 2016.

[49]
J. Leydold.
On the number of nodal domains of spherical harmonics.
Topology, 35(2):301–321, 1996.

[50]
B. Madafiglio.
Spectral geometry of the Laplacian and isoperimetric constants:
Cheeger's inequality.
Unpublished Honours thesis, 2019.

[51]
V. G. Maz'ya.
The negative spectrum of the $n$-dimensional Schrödinger
Dokl. Akad. Nauk SSSR, 144:721–722, 1962.

[52]
L. Miclo.
On eigenfunctions of Markov processes on trees.
Probability Theory and Related Fields, 142(3):561–594, Nov

[53]
L. Miclo.
On hyperboundedness and spectrum of Markov operators.
Inventiones mathematicae, 200(1):311–343, Apr 2015.

[54]
E. Milman.
On the role of convexity in isoperimetry, spectral gap and
Inventiones mathematicae, 177(1):1–43, Jul 2009.

[55]
A. Mukherjee.
Differential Topology.
Springer International Publishing, second edition, 2015.

[56]
E. Parini.
The second eigenvalue of the $p$-Laplacian as $p$ goes to 1.
International Journal of Differential Equations, Article ID
985671, 2010.

[57]
E. Parini.
An introduction to the Cheeger problem.
Surveys in Mathematics and its Applications, 6:9–21, 2011.

[58]
M. H. Protter and H. F. Weinberger.
Maximum Principles in Differential Equations.
Springer-Verlag New York, 1984.

[59]
A. Volkmann.
Regularity of isoperimetric hypersurfaces with obstacles in
Riemannian manifolds.
Diploma thesis, Albert-Ludwigs-Universitat Freiburg, 2010.

[60]
S.-T. Yau.
Isoperimetric constants and the first eigenvalue of a compact
Riemannian manifold.
Annales scientifiques de l'École Normale Supérieure, Ser.
4, 8(4):487–507, 1975.


# Subdiffusion with particle immobilization process described by differential
equation with Riemann–Liouville type fractional time derivative

Tadeusz Kosztołowicz<EMAIL_ADDRESS>Institute of Physics,
Jan Kochanowski University,
Uniwersytecka 7, 25-406 Kielce, Poland

###### Abstract

An equation describing subdiffusion with possible immobilization of particles
is derived by means of the continuous time random walk model. The equation
contains a fractional time derivative of Riemann–Liouville type which is a
differential-integral operator with the kernel defined by the Laplace
transform. We propose the method for calculating the inverse Laplace transform
providing the kernel in the time domain. In the long time limit the
subdiffusion–immobilization process reaches a stationary state in which the
probability density of a particle distribution is an exponential function.

##  I Introduction

In diffusion process particles can be eliminated from further diffusion in
different ways. There may be a particle decay due to a reaction when it meets
other molecules. Since the particle disappears, the probability density
$P(x,t)$ that the particle is at a point $x$ in time $t$ is not normalized,

$\int_{-\infty}^{\infty}P(x,t)dx<1.$ (1)

Another process that eliminates a particle from further diffusion is the
permanent immobilization of the particle. Both processes mentioned above can
occur in the diffusion of antibiotic molecules in a bacterial biofilm. One of
defense mechanisms is to disintegrate the antibiotic molecules, the process
can be described by diffusion-reaction equations. In the other one bacteria
can thicken the biofilm immobilizing antibiotic molecules aot ; mot , see also
km ; kmwa  and the references cited therein. The immobilized molecules have
not disappeared, they can further interact with the environment. In this case,
the probability of finding a molecule in the system is equal to one at any
time. We call the process subdiffusion with particle immobilization. It is
obvious that this process cannot be described by a diffusion–reaction
equation.

The immobilization of molecules can occur in a medium in which the movement of
particles is very hindered, as in the biofilm mentioned above, subdiffusion
may occur in such a system, see for example Refs. mk ; mk1 ; mks ; bg ; ks ;
klages2008 ; skb ; sk ; barkai2000 ; barkai2012 . We derive an equation
describing subdiffusion with particle immobilization in a one-dimensional
homogeneous system. We assume that after each jump a particle can be
immobilized with the same probability which does not change with time and is
independent of the particle position.

##  II Model

To derive the subdiffusion–immobilization equation we use the continuous time
random walk (CTRW) model mk ; mk1 ; mks ; skb ; barkai2000 ; ks ; montroll1965
; compte ; hilferanton ; chechkin . Within the model, when the average length
of a single particle jump $\epsilon$ is finite the form of the subdiffusion
equation is determined by the probability density $\psi$ of the waiting time
for the particle to jump. In terms of the Laplace transform,
$\mathcal{L}[f(t)](s)=\int_{0}^{\infty}{\rm e}^{-st}f(t)dt\equiv\hat{f}(s)$,
the equation generated by the function $\psi$ is as follows

$s\hat{P}(x,s)-P(x,0)=\frac{\epsilon^{2}s\hat{\psi}(s)}{2[1-\hat{\psi}(s)]}\frac{\partial^{2}\hat{P}(x,s)}{\partial
x^{2}},$ (2)

the derivation of this equation is described in Appendix.

We make the following assumptions:

  1. 1.

the probability of finding a particle in the system is equal to one at any
time,

$\int_{-\infty}^{\infty}P(x,t)dx=1,$ (3)

  2. 2.

since the particle can be permanently stopped, the probability that the
particle will make a jump is less than one,

$\int_{0}^{\infty}\psi(t)dt<1.$ (4)

###  II.1 Subdiffusion equation

To obtain the subdiffusion equation we assume

$\hat{\psi}(s)=\frac{1}{1+\tau s^{\alpha}},$ (5)

$0<\alpha<1$, where $\tau$ is a parameter with the units of
$\text{s}^{\alpha}$. This function satisfies the normalization condition

$\int_{0}^{\infty}\psi(t)dt\equiv\hat{\psi}(0)=1.$ (6)

This condition means that the particle cannot be stopped permanently with non-
zero probability. From Eqs. (2) and (5) we get

$s\hat{P}(x,s)-P(x,0)=Ds^{1-\alpha}\frac{\partial^{2}\hat{P}(x,s)}{\partial
x^{2}},$ (7)

where $D=\epsilon^{2}/2\tau$ is a subdiffusion coefficient given in the units
of $\text{m}^{2}/\text{s}^{\alpha}$. Due to the relations

$\mathcal{L}^{-1}\left[s\hat{f}(s)-f(0)\right](t)=\frac{df(t)}{dt},$ (8)
$\mathcal{L}^{-1}\left[s^{\beta}\hat{f}(s)\right](t)=\frac{{}^{RL}d^{\beta}f(t)}{dt^{\beta}},$
(9)

$0<\beta<1$, where

$\frac{{}^{RL}d^{\beta}f(t)}{dt^{\beta}}=\frac{1}{\Gamma(1-\beta)}\frac{d}{dt}\int_{0}^{t}(t-u)^{-\beta}f(u)du$
(10)

is the Riemann–Liouville time fractional derivative of the order
$\beta\in(0,1)$. From Eqs. (7)–(9) we get the subdiffusion equation

$\frac{\partial P(x,t)}{\partial
t}=D\frac{{}^{RL}\partial^{1-\alpha}}{\partial
t^{1-\alpha}}\frac{\partial^{2}P(x,t)}{\partial x^{2}}.$ (11)

###  II.2 Subdiffusion–immobilization equation

In order to find a function $\psi(t)$ that satisfies Eq. (4), i.e.
$\hat{\psi}(0)<1$, we assume that the Laplace transform of the function is

$\hat{\psi}(s)=\frac{1}{1+\tau\gamma+\tau s^{\alpha}},$ (12)

$0<\alpha<1$, the parameter $\gamma$, which controls molecule immobilization,
is given in the units of $1/\text{s}^{\alpha}$. The probability $p_{s}$ of
stopping the molecule permanently is
$p_{s}=1-\hat{\psi}(0)=\tau\gamma/(1+\tau\gamma)$. From Eqs. (2) and (12) we
get

$s\hat{P}(x,s)-P(x,0)=D\frac{s^{1-\alpha}}{1+\gamma
s^{-\alpha}}\frac{\partial^{2}\hat{P}(x,s)}{\partial x^{2}}$ (13)

The inverse Laplace transform of the right-hand side of Eq. (13) is calculated
using the formula

$\mathcal{L}^{-1}\left[\frac{s^{1-\alpha}}{1+\gamma
s^{-\alpha}}\hat{f}(s)\right](t)=\frac{{}^{RL}_{\;\;F}d^{1-\alpha}f(t)}{dt^{1-\alpha}},$
(14)

where

$\frac{{}^{RL}_{\;\;F}d^{1-\alpha}f(t)}{dt^{1-\alpha}}=\frac{d}{dt}\int_{0}^{t}F_{\alpha}(t-t^{\prime};\gamma)f(t^{\prime})dt^{\prime}$
(15)

is the Riemann–Liouville type fractional derivative with the kernel
$F_{\alpha}$ which is defined by its Laplace transform

$\hat{F}_{\alpha}(s;\gamma)=\frac{1}{\gamma+s^{\alpha}}.$ (16)

For $\gamma=0$, this derivative is the Riemann–Liouville derivative Eq. (10)
of the order $1-\alpha$. Eqs. (13)–(16) provide the following
subdiffusion–immobilization equation

$\frac{\partial P(x,t)}{\partial
t}=D\frac{{}^{RL}_{\;\;F}\partial^{1-\alpha}}{\partial
t^{1-\alpha}}\frac{\partial^{2}P(x,t)}{\partial x^{2}}.$ (17)

Calculation of the inverse transform of Eq. (16) is usually done by power
series expansion of the function when $\gamma/s^{\alpha}<1$, and then
inverting the transform term by term using the formula
$\mathcal{L}^{-1}[1/s^{\beta}](t)=t^{\beta-1}/\Gamma(\beta)$, $\beta>0$. The
result is the Mittag-Leffler function mainardi1 ; mainardi2 . However, this
procedure is valid for relatively large values of the parameter $s$, which
correspond to small values of time variable. To get the inverse Laplace
transform over the whole time domain we propose to use the following method:
(1) instead of $\hat{F}_{\alpha}$ Eq. (16) find the inverse transform of
$\hat{F}_{\alpha}(s,\gamma){\rm e}^{-as^{\mu}}$, $a,\mu>0$, (2) expand
$\hat{F}_{\alpha}$ in a power series of $s$ considering both cases
$s^{\alpha}>\gamma$ and $s^{\alpha}<\gamma$ separately, (3) use the formula
tkoszt2004

$\displaystyle\mathcal{L}^{-1}\left[s^{\nu}{\rm e}^{-as^{\mu}}\right](t)\equiv
f_{\nu,\mu}(t;a)$ (18)
$\displaystyle=\frac{1}{t^{\nu+1}}\sum_{n=0}^{\infty}\frac{1}{n!\Gamma(-n\mu-\nu)}\left(-\frac{a}{t^{\mu}}\right)^{n}$

$a,\mu>0$, (4) calculate the limit of $a\rightarrow 0^{+}$ in the obtained
functions. We note that

$f_{\nu,\mu}(t;0^{+})=\frac{1}{t^{\nu+1}\Gamma(-\nu)},$ (19)

and the result is independent of the parameter $\mu$.

From the formula

$\displaystyle\frac{{\rm
e}^{-as^{\mu}}}{\gamma+s^{\alpha}}=\left\\{\begin{array}[]{c}{\rm
e}^{-as^{\mu}}\sum_{n=0}^{\infty}(-\gamma)^{n}s^{-(n+1)\alpha},\;s>\gamma^{1/\alpha},\\\
\\\ \frac{{\rm
e}^{-as^{\mu}}}{\gamma}\sum_{n=0}^{\infty}\left(-\frac{1}{\gamma}\right)^{n}s^{n\alpha},\;s<\gamma^{1/\alpha},\end{array}\right.$
(23)

and Eqs. (18) and (19) we obtain

$\displaystyle
F_{\alpha}(t;\gamma)=\left\\{\begin{array}[]{c}\frac{1}{t^{1-\alpha}}E_{\alpha,\alpha}(-\gamma
t^{\alpha}),\;t<t_{b},\\\ \\\
-\frac{1}{\gamma^{2}t^{1+\alpha}}\tilde{E}_{\alpha,\alpha}\left(-\frac{1}{\gamma
t^{\alpha}}\right),\;t>t_{b},\end{array}\right.$ (27)

where $E_{\alpha,\beta}(u)=\sum_{n=0}^{\infty}\frac{u^{n}}{\Gamma(\alpha
n+\beta)}$, $\alpha,\beta>0$, is the two–parameter Mittag–Leffler (ML)
function,
$\tilde{E}_{\alpha,\beta}(u)=\sum_{n=0}^{\infty}\frac{u^{n}}{\Gamma(-\alpha
n-\beta)}$ is a generalization of the ML function for negative parameters. We
note that conditions $s>\gamma^{1/\alpha}$ and $s<\gamma^{1/\alpha}$ do not
determine the parameter $t_{b}$. For example, the condition
$s>\gamma^{1/\alpha}$ is equivalent to
$1/s^{\beta+1}<1/(s^{\beta}\gamma^{1/\alpha})$ for $\beta>0$ (assuming that
$s$ is a real positive parameter). The inverse Laplace transform of the
inequality provides $t<\beta/\gamma^{1/\alpha}$ where $\beta$ is a positive
number. Thus, the above inequality does not determine $t_{b}$. Here we define
the parameter $t_{b}$ as the shorter time at which the upper and the lower
functions in Eq. (27) are matched, see Fig. 1.

In terms of the Laplace transform the solution to Eq. (16) (the Green’s
function) for the initial condition $P(x,0)=\delta(x)$, where $\delta$ is the
Dirac–delta function, and boundary conditions $P(\pm\infty,t)=0$ is

$\hat{P}(x,s)=\frac{\sqrt{\gamma+s^{\alpha}}}{2s\sqrt{D}}\;{\rm
e}^{-|x|\frac{\sqrt{\gamma+s^{\alpha}}}{\sqrt{D}}}$ (28)

The solution fulfils the condition $\int_{-\infty}^{\infty}\hat{P}(x,s)dx=1/s$
what provides the normalization of the function $P$ Eq. (3).

Let $\gamma\neq 0$. We calculate the inverse Laplace transform of the function
(28) for small and large values of $s$ separately. In calculation, we use the
formulas $\sqrt{1+u}\approx 1+u/2-u^{2}/8$ and ${\rm e}^{-u}\approx
1-u+u^{2}/2$, $u\rightarrow 0$, and keep the leading terms in the obtained
series. When $s^{\alpha}>\gamma$ we obtain

$\hat{P}(x,s)=\frac{1}{2\sqrt{D}s^{1-\alpha/2}}\left(1-\frac{b_{1}}{s^{\alpha/2}}+\frac{b_{2}}{s^{\alpha}}\right){\rm
e}^{-\frac{|x|}{\sqrt{D}}s^{\alpha/2}},$ (29)

where $b_{1}=\gamma|x|/2\sqrt{D}$ and
$b_{2}=(\gamma/2)(1+|x|^{2}\gamma/2\sqrt{D})$. If $s^{\alpha}<\gamma$, we get

$\hat{P}(x,s)=\frac{\sqrt{\gamma}}{2s\sqrt{D}}{\rm
e}^{-\sqrt{\frac{\gamma}{D}}|x|(1+\frac{s^{\alpha}}{2\gamma})}\left[1+\frac{s^{\alpha}}{2\gamma}-b\frac{s^{2\alpha}}{\gamma^{2}}\right],$
(30)

where $b=\sqrt{\gamma/D}|x|+1/8$. Eqs. (18) and (29) provide the Green’s
functions in the limit of short time

$\displaystyle
P(x,t)=\frac{1}{2\sqrt{D}}\Big{[}f_{-1+\alpha/2,\alpha/2}(t;\eta)$ (31)
$\displaystyle-
b_{1}f_{-1,\alpha/2}(t;\eta)+b_{2}f_{-1-\alpha/2,\alpha,2}(t;\eta)\Big{]},$

where $\eta=|x|/\sqrt{D}$. From Eqs. (18) and (30) we get the Green’s function
in the long time limit

$\displaystyle P(x,t)=\frac{1}{2}\sqrt{\frac{\gamma}{D}}{\rm
e}^{-\sqrt{\frac{\gamma}{D}}|x|}\Big{[}f_{-1,\alpha}(t;\xi)$ (32)
$\displaystyle+\frac{1}{2\gamma}f_{\alpha-1,\alpha}(t;\xi)-\frac{b}{\gamma^{2}}f_{2\alpha-1,\alpha}(t;\xi)\Big{]},$

where $\xi=|x|/2\sqrt{D\gamma}$.

Since the mean particle position equals zero, in terms of the Laplace
transform the mean square displacement of the particle is

$\mathcal{L}\left[\left\langle\right(\Delta
x)^{2}(t)\rangle\right](s)=\int_{-\infty}^{\infty}x^{2}\hat{P}(x,s)dx=\frac{2D}{s(\gamma+s^{\alpha})}.$
(33)

When $\gamma\neq 0$, for small $s$ we have
$\mathcal{L}\left[\left\langle\right(\Delta
x)^{2}(t)\rangle\right](s)=2D/[1/s-1/(\gamma s^{1-\alpha})]$. Thus, in the
limit of long time we get

$\left\langle(\Delta
x)^{2}(t)\right\rangle=\frac{2D}{\gamma}\left[1-\frac{1}{\gamma\Gamma(1-\alpha)t^{\alpha}}\right].$
(34)

In the limit $t\rightarrow\infty$, the stationary state described by the
following function is reached,

$P(x,t\rightarrow\infty)\equiv
P_{st}(x)=\frac{1}{2}\sqrt{\frac{\gamma}{D}}\;{\rm
e}^{-\sqrt{\frac{\gamma}{D}}|x|}.$ (35)

Figure 1: Plot of the function $F_{\alpha}$. The dashed vertical line shows
the location of the parameter $t_{b}=11.5$. The solid line with squares is the
plot of the upper function in Eq. (27) which describes $F_{\alpha}$ for
$t<t_{b}$, the solid line with circles is the plot of the lower function in
Eq. (27) which represents $F_{\alpha}$ for $t>t_{b}$. In the numerical
calculations, the leading 20 terms in the series appearing in the functions
$E_{\alpha,\alpha}$ and $\tilde{E}_{\alpha,\alpha}$ have been taken into
account. Figure 2: Plots of Green’s functions for times given in the legend.
The plots represent the function Eq. (29) for $t=0.1,0.5$ and Eq. (30) for
$t=15,50,100$. Figure 3: Plots of the function $P_{st}$ Eq. (35) for different
values of the ratio $\gamma/D$ given in the legend.

For illustration, plots of functions $F_{\alpha}$ and $P$ are shown in Figs. 1
and 2, respectively. The parameters are $\alpha=0.7$, $\gamma=0.6$, and
$D=10$, all parameters are given in arbitrarily chosen units. In Fig. 3 the
Green’s functions for the stationary state are presented.

##  III Final remarks

The process of subdiffusion with particle immobilization can be described by
an equation with a fractional time derivative of the Riemann–Liouville type,
which is a differential-integral operator with the kernel $F_{\alpha}$ defined
by its Laplace transform Eq. (16). Normal diffusion and subdiffusion have a
different stochastic interpretation. However, the normal diffusion-
immobilization equation can be obtained from Eq. (17) by substituting
$\alpha=1$. We have proposed a method for determining the inverse Laplace
transform of the kernel. In our opinion, this method can be widely used for
calculating inverse Laplace transforms $\mathcal{L}^{-1}[\hat{f}(s)](t)$ for a
wide class of functions $f$.

In a homogeneous unbounded system the subdiffusion-immobilization process
reaches a stationary state which is described by $P_{st}(x)$ Eq. (35). This
distribution depends only on the quotient $\gamma/D$ expressed in the units of
$1/\text{m}^{2}$ and it does not explicitly depend on the parameter $\alpha$.
The achievement of the steady state is suggested by Fig. 2, where the Green’s
functions for relatively long times differ very little from each other. In the
stationary state there is $\left\langle(\Delta
x)^{2}(t\rightarrow\infty)\right\rangle=\frac{2D}{\gamma}$, the particle is
finally immobilized with probability equal to one.

The subdiffusion–immobilization process is described by Eq. (17) that can be
obtained in practice by replacing the time fractional Riemann-Liouville
derivative Eq. (10) with the more general Riemann–Liouville type derivative
with the kernel $F_{\alpha}$ Eq. (15) in the “ordinary” subdiffusion equation
Eq. (11), orders of both derivatives are the same. There is a different
situation than in the subdiffusion–reaction equation in which the reaction
term is involved in the “ordinary” subdiffusion equation, see for example
Refs. mendez ; yuste ; kl2014 . We mention that the Riemann–Liouville type
fractional derivatives with different kernels have been considered in Ref.
hl1998 ; lh1999 ; lh2000 ; yang .

## Acknowledgment

The author wishes to express his thanks to Aldona Dutkiewicz for fruitful
discussions.

## Appendix. Derivation of Eq. (2)

To derive the subdiffusion equation we use a simple model of a particle random
walk along a one–dimensional homogeneous lattice. Usually, in the CTRW model
both a particle jump length and waiting time for a particle to jump are random
variables. We assume that the jump length distribution $\lambda$ has the form
$\lambda(x)=\frac{1}{2}[\delta(x-\epsilon)+\delta(x+\epsilon)]$. Random walk
with discrete time $n$ is described by the equation
$P_{n+1}(m)=\frac{1}{2}P_{n}(m+1)+\frac{1}{2}P_{n}(m-1)$, where $P_{n}(m)$ is
a probability that a diffusing particle is at the position $m$ after $n$-th
step. Let the initial particle position be $m=0$. Moving from discrete $m$ to
continuous $x$ spatial variable we assume $x=m\epsilon$ and
$P_{n}(x)=P_{n}(m)/\epsilon$, where $\epsilon$ is a distance between discrete
sites. The above equations and the relation
$[P_{n}(x+\epsilon)+P_{n}(x-\epsilon)-2P_{n}(x)]/\epsilon^{2}=\partial^{2}P_{n}(x)/\partial
x^{2}$, $\epsilon\rightarrow 0$, provide the following equation in the limit
of small $\epsilon$

$P_{n+1}(x)-P_{n}(x)=\frac{\epsilon^{2}}{2}\frac{\partial^{2}P_{n}(x)}{\partial
x^{2}}.$ (36)

To move from discrete to continuous time we use the formula
$P(x,t)=\sum_{n=0}^{\infty}Q_{n}(t)P_{n}(x)$ montroll1965 , where $Q_{n}(t)$
is the probability that a diffusing particle takes $n$ step in the time
interval $(0,t)$. The function $Q_{n}$ is a convolution of $n$ distributions
$\psi$ of a waiting time for a particle to jump and a function
$U(t)=1-\int_{0}^{t}\psi(t^{\prime})dt^{\prime}$ which is the probability that
a particle does not change its position after $n$-th step,
$\hat{U}(s)=[1-\hat{\psi}(s)]/s$,
$Q_{n}(t)=(\underbrace{\psi\ast\psi\ast\ldots\ast\psi}_{n\;times}\ast U)(t)$,
where $(f\ast h)(t)=\int_{0}^{t}f(u)h(t-u)du$. Due to the following property
$\mathcal{L}[(f\ast h)(t)](s)=\hat{f}(s)\hat{h}(s)$ we obtain

$\hat{P}(x,s)=\frac{1-\hat{\psi}(s)}{s}\sum_{n=0}^{\infty}\hat{\psi}^{n}(s)P_{n}(x).$
(37)

Combining Eqs. (36) and (37) we get Eq. (2).

## References

  * (1) G. G. Anderson and G. A. O’Toole, Bacterial Biofilms, Current Topics in Microbiology and Immunology 322 (Berlin, Springer, 2008).
  * (2) T. F. C. Mah and G. A. O’Toole, Mechanisms of biofilm resistance to antimicrobial agents, Trends Microbiol. 9, 34 (2001).
  * (3) T. Kosztołowicz and R. Metzler, Diffusion of antibiotics through a biofilm in the presence of diffusion and absorption barriers, Phys. Rev. E 102, 032408 (2020).
  * (4) T. Kosztołowicz, R. Metzler, S. Wa̧sik, and M. Arabski, Modelling experimentally measured of ciprofloxacin antibiotic diffusion in Pseudomonas aeruginosa biofilm formed in artificial sputum medium, PLoS ONE 15(12), e0243003 (2020).
  * (5) R. Metzler and J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Phys. Rep. 339, 1 (2000).
  * (6) R. Metzler and J. Klafter, The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics, J. Phys. A 37, R161 (2004).
  * (7) R. Metzler, J. Klafter, and I. M. Sokolov, Anomalous transport in external fields: Continuous time random walks and fractional diffusion equations extended, Phys. Rev. E 58, 1621 (1998).
  * (8) J. P. Bouchaud and A. Georgies, Anomalous diffusion in disordered media: statistical mechanisms, models and physical applications, Phys. Rep. 195, 127 (1990).
  * (9) J. Klafter and I. M. Sokolov, First Step in Random Walks. From Tools to Applications (Oxford UP, New York, 2011).
  * (10) R. Klages, G. Radons, and I. M. Sokolov, Anomalous Transport: Foundations and Applications (Wiley, New York, 2008).
  * (11) I. M. Sokolov, J. Klafter, and A. Blumen, Fractional kinetics, Phys. Today 55, 11, 48-54 (2002).
  * (12) I. M. Sokolov and J. Klafter, From diffusion to anomalous diffusion: a century after Einstein’s Brownian motion, Chaos 15, 026103 (2005).
  * (13) E. Barkai, R. Metzler, and J. Klafter, From continuous time random walks to the fractional Fokker-Planck equation, Phys. Rev. E 61, 132 (2000).
  * (14) E. Barkai, Y. Garini, and R. Metlzer, Strange kinetics of single molecules in living cells, Phys. Today 65, 29 (2012).
  * (15) E. W. Montroll and G. H. Weiss, Random walks on lattices. II, J. Math. Phys. 6, 167 (1965).
  * (16) A. Compte, Stochastic foundations of fractional dynamics, Phys. Rev. E 53, 4191 (1996).
  * (17) R. Hilfer and L. Anton, Fractional master equations and fractal time random walks, Phys. Rev. E 51, R848 (1995).
  * (18) A. V. Chechkin, M. Hofmann, and I. M. Sokolov, Continuous-time random walk with correlated waiting times, Phys. Rev. E 80, 031112 (2009).
  * (19) F. Mainardi, A tutorial on the basic special functions of fractional calculus, WSEAS Trans. Math. 19, 74 (2020).
  * (20) F. Mainardi, Why the Mittag-Leffler function can be considered the Queen function of the fractional calculus?, Entropy 22, 1359 (2020).
  * (21) T. Kosztołowicz, From the solutions of diffusion equation to the solutions of subdiffusive one, J. Phys. A: Math. Gen. 37, 10779 (2004).
  * (22) V. Méndez, S. Fedotov, and W. Horsthemke, Reaction–Transport Systems: Mesoscopic Foundations, Fronts, and Spatial Instabilities (Springer, Berlin, 2010).
  * (23) S.B. Yuste, L. Acedo, and K. Lindenberg, Reaction front in an $A+B\rightarrow C$ reaction-subdiffusion process, Phys. Rev. E 69, 036126 (2004).
  * (24) T. Kosztołowicz and K.D. Lewandowska, Subdiffusion-reaction processes with $A\rightarrow B$ reactions versus subdiffusion-reaction processes with $A+B\rightarrow B$ reactions, Phys. Rev. E 90, 032136 (2014).
  * (25) T.T. Hartley and C.F. Lorenzo, A solution to the fundamental linear fractional order differential equation, NASA/TP-1998-208693 (1998).
  * (26) C.F. Lorenzo and T.T. Hartley, Generalized functions for the fractional calculus, NASA/TP-1999-209424/REV1 (1999).
  * (27) C.F. Lorenzo and T.T. Hartley, R–function relationships for application in the fractional calculus, NASA/TM-2000-210361 (2000).
  * (28) X.-J. Yang, General Fractional Derivatives. Theory. Methods and Applications, CRC Press, Taylor and Francis Group, Boca Raton (2019), p.177.