# Influence of Orbital Angular Momentum of light on Random Spin-split modes in
Disordered Anisotropic Optical media

Anwesha Panda  Department of Physical Sciences
Indian Institute of Science Education and Research Kolkata
Mohanpur, India - 741246  Sneha Dey  Department of Physics, University of
Calcutta, 92, APC Rd, Kolkata, India - 700009  Yogishree Arabinda Panda
Cluster of Excellence Physics of Life, Dresden University of Technology,
Dresden, Germany - 01307  Aditya Anurag Dash  Department of Physical Sciences
Indian Institute of Science Education and Research Kolkata
Mohanpur, India - 741246  Aloke Jana  Department of Physics, Colorado School
of Mines, 1500 Illinois St, Golden, CO 80401, USA  Nirmalya Ghosh  Department
of Physical Sciences
Indian Institute of Science Education and Research Kolkata
Mohanpur, India - 741246<EMAIL_ADDRESS>
###### Abstract

Spin orbit interaction of light in a disordered anisotropic medium is known to
yield spin split modes in the momentum domain because of the random spatial
gradient of the geometric phase of light. Here, we have studied the statistics
of such spin split modes for beams carrying intrinsic orbital angular momentum
through the quantification of momentum domain entropy and investigated its
dependence on various beam parameters. The influence of the spatial structure
of the beam and the phase vortex on the statistics of the spin split modes
were separately investigated using input Laguerre-Gaussian and Perfect Vortex
beams passing through disordered anisotropic medium with controlled input
disorder parameter, which was realized by modulating the pixels of a liquid
crystal-based spatial light modulator. The results of systematic
investigations on the impact of beam waist, spot size and topological charge
of the vortex beam shows that the influence of the spot size on the emergence
of the random spin split modes is much more significant as compared to the
other beam parameters.
Key words: Spin-orbit interaction of light, Vortex beam, Spatial phase
gradient.

##  1 Introduction

Spin-orbit coupling, also known as spin-orbit interaction (SOI), is an
universal concept in physics that involves the coupling of spin and orbital
degrees of freedom in both particles with mass (e.g., electrons) and massless
particles (e.g., photons) due to relativistic effects. It occurs in various
systems, ranging from atomic, condensed-matter systems to optical
technologies, giving rise to interesting phenomena and potential applications.
However, it is intriguing to note that this phenomenon, although typically
discussed in the context of quantum particles, can also be manifested in
classical light. This follows that classical light beams can carry both spin
angular momentum (SAM) related to circular polarization and orbital angular
momentum (OAM) associated with helical wavefronts of light, both of which
interconvert into each other to create rich physics associated with SOI.

The SOI of light gives rise to two closely intertwined phenomena. The first
phenomenon involves the influence of the light’s trajectory on its state of
polarization, resulting in the emergence of spin-dependent optical vortices
which is commonly observed in systems with cylindrical or spherical symmetry.
The second phenomenon is the reciprocal effect where polarization affects the
trajectory of light. This effect is termed the Spin Hall effect (SHE) of light
and is typically observed when symmetry is broken[1, 2, 3, 4, 5, 6, 7, 8, 5,
9]. SOI of light are basically of two types - one is geometric phase mediated
and other is transverse angular momentum mediated. Geometric phases and their
gradients, along with the preservation of the total angular momentum of light,
are closely linked to optical Spin-Orbit Interaction (SOI) phenomena. Within
this context, two distinct forms of geometric phase play a role—the spin
redirection Berry phase and the Pancharatnam-Berry (PB) geometric phase. Spin
Hall effect (SHE) originating from the geometric phase gradient, leads to
either spatial domain or momentum domain shift. The other type of SHE, which
is completely independent of the geometrical phase of light, originates from
the transverse spin-angular momentum of light and is observed in the case of
surface waves, evanescent waves, and waveguide modes, is not investigated in
this paper [2, 3, 10, 11, 4, 12, 13].

These advancements have resulted in a variety of fundamental effects in the
realm of photonic SOI across diverse light-matter interactions. Remarkable
phenomena like spin-to-vortex transformation, orbital Hall effect, optical
Rashba effect, plasmonic Aharonov–Bohm effect, spin-dependent transverse
momentum, transverse SAM, spin-momentum locking, and spin-controlled
directional waveguiding. These breakthroughs have paved the way for novel
insights into universal SOI principles, offering new possibilities for
designing spin-orbit photonic-devices.[14, 15, 16, 17, 18, 19, 13]

Most of the scenarios described previously in the context of spin orbit
interactions specifically those dealing with geometrical phases are for
ordered inhomogeneous anisotropic medium. One of such realization is the
metasurface which is fabricated by spatially structuring anisotropic media in
the nanometer length scale. However, a recent discovery showcased the
possibility of obtaining spin-orbit-coupled random scattering modes across the
entire momentum range in a completely disordered, inhomogeneous, and
anisotropic optical system. This phenomenon, known as the random optical
Rashba effect[20, 21], is characterized by the presence of disordered spin-
orbit coupling throughout the beam profile or a disordered spatial
distribution of the geometric phase and its gradient [8]. Notably, the impact
of intrinsic orbital angular momentum beams on spin-split modes remains
unexplored, which is the focus of this investigation.

This paper aims to explore the influence of topological charge, spot size, and
spatial structure of the beam on resulting spin-split modes in disordered
anisotropic media. For this purpose, we have separately investigated using
Lagaurre-Gaussian (LG) and Perfect vortex (PV) beams. The main purpose of
using PV beams is to investigate the role of spot size and topological charges
separately, as the size of vortices is independent of topological charges for
PV beams. Here, we have quantified the statistics of spin-split modes in the
momentum domain by the standard Shannon entropy. The study shows the momentum
domain entropy of the spin-split modes are affected by both topological charge
and spatial structure.

The structure of this paper is as follows: Section 2 presents the theoretical
framework of LG beams and perfect vortices, along with the formation of spin-
split modes. Section 3 outlines the experimental procedure, while Section 4
discusses simulations and results of scattered modes in momentum space
entropy. Finally, Section 5 provides concluding remarks and a summary of the
findings.

##  2 Theory

Let $\ket{E_{i}}$ and $\ket{E_{o}}$ be the input and output electric fields
respectively, and the beam (containing RCP or LCP polarized light) is passing
through an inhomogeneous(in the transverse plane with coordinates
$\xi\rightarrow$x/y-plane and z being the propagation direction of light)
anisotropic medium, we get –

$\ket{E_{o}}=e^{i{\phi_{d}(\xi)\pm\phi_{g}(\xi)}}\ket{E_{i}}$ (1)

where $\phi_{d}(\xi)$ is dynamic phases and $\phi_{g}(\xi)$ is PB geometric
phase. When the two phases have equal gradient,
$\frac{d\phi_{d}(\xi)}{d\xi}$=$\frac{d\phi_{g}(\xi)}{d\xi}$=$\Omega_{\xi}$,
then the momentum domain shifts for RCP and LCP polarization states become

$<{k_{\xi}}>_{RCP}=2\Omega_{\xi}\hskip 10.00002ptand\hskip
10.00002pt<{k_{\xi}}>_{LCP}=0$ (2)

where

$<k_{\xi=x/y}>=\frac{<E_{0}|i\frac{\partial}{\partial\xi=x/y}|E_{0}>}{<E_{0}|E_{0}>}$

i.e., for RCP polarization, momentum domain shift becomes twice of spatial
gradient of phase (geometric or dynamic), and for LCP, momentum domain shift
does not occur [22, 23].

Dynamical phase and Pancharatnam-Berry (PB) geometric phase in twisted nematic
liquid crystal layers: Polarized light gains PB geometric phase as well as
dynamical phase while propagating in an anisotropic material. The dynamical
phase for a linear birefringent medium is determined by the extraordinary and
ordinary refractive indices ($n_{e}$ and $n_{o}$) and consequently it also
depends upon the magnitude of linear retardance $\delta$ (defined as $\delta$
= $\frac{2\pi}{\lambda}(n_{e}-n_{0})d$, where $d$ is the path length and
$\lambda$ is the wavelength). The PB geometric phase in such birefringent
medium, on the other hand, is determined by the orientation angle of the
anisotropy axis. Thus, in principle, one can produce equal spatial gradients
of the dynamical phase and PB geometric phase in an inhomogeneous birefringent
medium by controllably and simultaneously changing the magnitude of linear
retardance $\delta$ and the orientation angle of the anisotropy axis in the
transverse plane.

The above can be realized by one of the readily available system, which is a
twisted nematic liquid crystal-based spatial light modulator (SLM). The
evolution of polarization in SLM can also be alternatively modelled using the
effective Jones matrix ($J_{eff}$) as a sequential product of matrices of an
equivalent linear retarder ($J_{reta}$, with effective linear retardance
$\delta_{eff}$ and its orientation angle $\theta_{eff}$) and an effective
optical rotator (with optical rotation $\psi_{eff}$). The rotation does not
truly have a dynamical origin, it is actually related to the twist angle
($\psi$).[22, 24, 25, 26, 27, 28]

$J_{eff}=R(\psi_{eff})J_{reta}(\delta_{eff},\theta_{eff})$ (3)

where $\psi_{eff}=-\psi+2\theta_{eff}$

The total dynamical phase is primarily determined by the total linear
retardance $\delta$, while the PB geometric phase is determined by effective
optical rotation ($\psi_{eff}$). It was shown previously [22], for a certain
range of grey values (n=30-170), the two gradients, $\frac{d\delta(n)}{dn}$
and $\frac{d\psi_{eff(n)}}{dn}$ are equal. This leads to an equal gradient of
geometric and dynamic phases which results in equation 2. The variation of the
$\delta$, $\psi_{eff}$ parameters depends on the grey-level values (n). By
changing the grey level values (n) from 30 to 170, both geometric and
dynamical phases can be simultaneously tailored in an SLM. In our experiment
we have modulated the pixels of SLM by generating random grey values (n) using
a delta correlated uniformly distributed random function
$f^{\epsilon}(\phi_{g})=1/2\pi\epsilon$ for
$-\epsilon\pi\leq\phi_{g}(x)<\epsilon\pi$ otherwise
$f^{\epsilon}(\phi_{g})=0$, $0\leq\epsilon\leq 1,(where\epsilon$ controls the
amount of randomness of $x$ coordinate).(shown in Fig.1)

Figure 1:  The geometric phase distribution projected in the spatial light
modulator (SLM). a) the grid was divided into pixels in x and y directions,
where each pixel has a single random grey value (generated using delta
correlated function as mentioned above). b) grid is divided azimuthally into n
number of uniform divisions, with each azimuthal division having a single
random grey scale value (n).

The variation of phase in transverse plane, manifests as a distribution of
intensities in momentum domain, which can be attributed to either dynamical
phase, or geometrical phase, or a combination of both. Particularly, when this
distribution is generated through the influence of only geometrical phase
($\phi=\pm\phi_{g}$), the SHE, becomes observable for right circularly
polarized (RCP) and left circularly polarized (LCP) light. Here, the SOI of
light results from the light beam’s spatial inversion symmetry being broken by
an inhomogeneous distribution of a geometric and dynamic phase combination[1,
2, 3, 4, 5, 6, 7, 8, 5, 4, 29]. The strength of such an effect depends on the
phase inhomogeneity acquired by the light beam. The corresponding intensity
will be distributed throughout the momentum space as

$I(k_{x},k_{y})=|\iint^{+\infty}_{-\infty}e^{-i(k_{x}x+k_{y}y)}E_{o}(x,y)dxdy|^{2}.$
(4)

[23]. It is also important to note, as long as we can describe the system(the
anisotropic media) by a local phase gradient ($\xi$) , we will observe
momentum domain spin-Hall effect of light as per equation 2. But as the
randomness increases, beyond a certain value of disorderness when the system
cannot be described by a particular local gradient anymore, one would observe
random spin-split scattered modes in the momentum domain. To quantify the
statistics of the momentum domain intensity distribution of the random spin-
split scattered modes, we have defined the well-known Shannon entropy function
in the following way-

$H=-\sum_{i}p_{i}[I(k_{x},k_{y})]log(p_{i}[I(k_{x},k_{y})])$ (5)

where $p_{i}[I(k_{x},k_{y})]$ is the probability density function (PDF) of
intensity distribution of the scattered modes in the momentum domain[8].

For this purpose, in addition to the normal Gaussian beam we have used
Laguerre-Gaussian (LG) and Perfect Vortex (PV) beams with different $l$ values
(topological charges) . Perfect vortex beams have a constant spatial structure
regardless of their topological charge. This characteristic makes it a
valuable tool for decoupling the effects of the spatial beam structure and the
phase vortex.

Laguerre-Gaussian beam at the source plane z=0 has the form -

$E_{pl}(x,y,0)=(\frac{\sqrt{2}\rho}{\omega})^{l}L_{p}^{l}(\frac{2\rho^{2}}{\omega^{2}})exp(-\frac{\rho^{2}}{\omega^{2}}+il\phi)$
(6)

where $\rho=(x^{2}+y^{2})^{\frac{1}{2}}$ and $\phi=tan^{-1}(y_{0}/x_{0})$,
$\omega$ is beam waist. $L_{p}^{l}$ is the associated Laguerre polynomial
where p and $l$ are the radial modes and the angular modes (topological
charges) respectively.

An approximate model of the PV beam [30] to enable experiments is as follows:

$E_{PV}(\rho,\phi_{0})=exp[-\frac{(\rho-\rho_{0})^{2}}{\Delta\rho^{2}}]exp(il\phi_{0})$
(7)

where ($\rho,\phi_{0}$) are the polar coordinates in beam cross section, $l$
is the topological charge $\rho_{0}$ is the radius of annular bright
intensity, $\Delta\rho$ is a small width. The Fourier transformation of an
ideal Bessel beam function [31] may be used to calculate the approximate model
of PV beams, and it can be written as

$E_{BG}(\rho,z)=J_{l}(k_{\rho}\rho)exp(il\phi_{0}+ik_{z}z)$ (8)

where $J_{l}$ is the first kind of $l^{th}$ order bessel function,
$k=\sqrt{k_{r}^{2}+k_{z}^{2}}=2\pi/\lambda,r=(\rho,\phi_{0})$ and
$(k_{r},k_{z})$ are the radial and longitudinal wave vectors respectively.

##  3 Experimental Methods

A schematic of the experimental setup for observing spin-split scattering
modes for LG beam is shown in Figure 2(a). A fundamental Gaussian mode of a
632.8 nm line of a He-Ne laser (HNL050L, Thorlabs, 5mW power) is used in this
set-up. The beam is transmitted through the spatial light modulator (SLM,
LC2012). A computer-generated fork hologram for different topological charges
($\pm$l) is projected to the SLM1 to produce an orbital angular momentum
carrying beam (LG). After passing through the SLM1, the central beam (zeroth
order) remains Gaussian, and its adjacent beams (first order) are LG with
$\pm$l topological charges. An aperture is used after SLM to select and pass
only one LG beam (+$l$ or -$l$) according to necessity. After being reflected
by two mirrors (M1 and M2), the beam passes through PSG (Polarization State
Generator) unit comprised of a Glan-Thomson linear polarizer (P1, GTH10M,
Thorlabs, USA). The PSA (Polarization-State-Analyzer) unit consists of a
linear polarizer and a quarter-waveplate, but positioned in the reverse order
for selecting LCP or RCP polarization. In the middle of PSG and PSA, we have
spatial light modulator (SLM2, same model as SLM1) realized as random
inhomogeneous anisotropic media (shown in Figures 1a and 1b). Finally, there
is a lens to obtain the Fourier image in the momentum domain. A CCD
(1024$\times$768 square pixels, pixel dimension 4.65 $\mu$m, Thorlabs, USA)
has been used at the end to collect the light at the Fourier plane. Spin-
selective random scattering modes are observed for RCP polarization only and
LCP polarization does not show any results. Figure 2(b) and 2(c) are the
momentum domain intensity distribution for input LG beam while Figure 2(d) and
2(f) are for input PV beam .

In the usual case of spin-Hall effect of light, the input beam selected by PSG
is linearly polarized. In SLM2, we have synchronously modulated the geometric
and dynamic phases as shown in [22]. For one circularly polarized light (RCP),
the geometric and dynamic phases add up; for the other, they cancel out giving
no effect for disorder. This is called spin-selective scattering modes. To
observe these modes, the inhomogeneous anisotropic medium was made possible by
modulating the pixels of the SLM2 by randomly distributed grey values
(n=30-170). The spin-selective scattering modes have been observed for all the
beams such as Gaussian, LG, and perfect vortex. Figure 2(b),(c) shows this
effect for LG beam with topological charge $l$=5 and Figure 2(d),(e) is for PV
beam with topological charge $l$=5.

To build a setup for perfect vortex (PV), we need an additional lens and
axicon in the above setup. Axicon lens converts LG to Bessel Gauss
consequently Lens L1 has been used as a Fourier lens to transforms Bessel
Gauss to perfect vortex. The position of lens L1 and axicon is commutative. PV
beam is obtained exactly at the focal plane of lens L1 where the SLM2
[projecting the disordered anisotropic media (shown in Figure 1)] has to be
placed accurately; a slight shift from this plane gives us bessel gauss. To
obtain the Fourier transform lens L2 also has to be placed in the Fourier
plane of L1. This was an experimental challenge to adjust L2 and SLM2 in the
same plane, keeping the distance between them only less than 1 millimeter. We
observed spin-selective scattering modes for perfect vortex as well similar to
LG or Gaussian beam.

In order to calculate the momentum space entropy (Shannon entropy) in the
Fourier plane of LG or PV beam, we calculate the normalized probability
density function (PDF) of intensity distribution collected by CCD camera of
1024$\times$768 pixels. This normalized PDF is used in the expression of
Shannon entropy using Eqn 5 [8], to get the momentum space entropy.

Figure 2:  (a) Schematic of the experimental arrangement for observing spin-
selective random scattering modes using Gaussian and Laguerre-Gaussian (LG)
beams. He-Ne Laser: light source; SLM1, SLM2: spatial light modulators (SLM1
was used to generate the LG beam with different topological charge, SLM2 was
used to realize the disordered anisotropic media); A: aperture; M1, M2:
mirrors; P1, P2: polarizers; QWP: quarter wave-plate; L: lens; CCD: camera
(placed in the Fourier plane to record momentum domain) here. The
polarization-state analyzer (PSA) unit comprises a linear polarizer and
quarter waveplate. This unit analyzes desirable polarization states of light
(LCP or RCP in this case). Spin-selective scattering modes of LG beam
transmitted through disordered anisotropic media of $\epsilon$=1 are shown for
topological charge $l$=5 for RCP (b) and LCP (c) by projecting 50 azimuthal
divisions in SLM2 (each division has different grey-scale values (n) which
have been assigned randomly). Similar result are shown for input PV beam of
topological charge $l$ = 5 passing through the same grey scale pattern
projected on SLM2 as mentioned above; in (c) LCP and (d) RCP. The spin
selective property of the random inhomogeneous media [23] holds not only for
Gaussian beam but also for LG and PV beams. (f) Schematic of the experimental
setup to observe spin-selective scattering modes of input perfect vortex beams
Most components are similar to the setup shown in (a). Additional components
to produce the PV beam are:- L1, Axicon. L1 , L2: lenses which are used as
Fourier lenses; Axicon lens converts LG beam to Bessel Gauss.

##  4 Results and Discussions

At first, we consider the effect of input LG beam with the varying $l$ values
and input disorder parameter $\epsilon$ of the disordered anisotropic optical
media. As shown in Figure 3, which displays the experimentally observed
momentum space intensity distributions of Gaussian and LG beam with $l$=0, 3,
5, and with varying input disorder parameters, $\epsilon$ = 0, 0.5, 1. It is
observed as the input disorder parameter $\epsilon$ increases, the number of
random scattered spin-split modes increases gradually which indicates increase
in momentum space entropy. This variation has been quantified by using
equation 5 and is presented in Fig 4.

Figure 3:  Variation of momentum domain intensity distribution of Gaussian and
Laguerre-Gaussian beams for different topological charges ($l$=0, 3, 5)
transmitting through disordered anisotropic media with input disorder
parameter $\epsilon$ = 0, 0.5, 1. As the $\epsilon$ value increases, the beam
gets more and more scattered. SLM2 was projected with 50 azimuthal divisions
of random grey scale (n). Results shown are only for RCP polarization selected
using PSA.

Figure 4, depicts the variation of momentum space entropy with input
$\epsilon$ parameter along with varying topological charges $l$ for input LG
beam. We observe a sudden surge in momentum space entropy at a critical value
of the disorder parameter $\epsilon$ i.e., there exists a threshold value of
$\epsilon$. This provides strong evidence for existence of a phase transition
within the system under study. This occurs because, as long as a local spatial
phase gradient (both geometrical and dynamical) can be clearly defined within
the system, we anticipate the Spin Hall Shift. However, as the heterogeneity
of the optical media intensifies, a critical point is reached where the local
spatial phase gradient can no longer be distinctly defined. i.e., the media
becomes completely random. We observe the emergence of random spin split
modes. In Fig. 4 (a) and (b) we see a clear threshold value of the disorder
parameter $\epsilon$ around 0.25. After this, as input disorder parameter
$\epsilon$ increases the variation of momentum space entropy increases
rapidly. Also, we see greater $l$ values imply greater total momentum space
entropy for a input disorder parameter $\epsilon$. It is vital to note that as
one varies the $l$ values, one would not only vary the topological charge but
also the effective vortex size or the input spot size of the beam in the
spatial domain. This correspondingly changes the amount of inhomogeneity
probed by the input beam in the spatial domain and consequently their momentum
domain span. It may so happen that $l$ influences the spin split modes but in
addition to that the effective spot size will also affect them as larger spot
size will probe greater area of inhomogenity. This is why we anticipated that
not only $l$ but the spin split modes will also be affected by the spot size
of the input beam. Thus, it is necessary to investigate the role of beam spot
size separately.

Figure 4:  Momentum space entropy ($H$) variation of LG beams with different
topological charges , $l$ = 0, 3, 5, 7, 10 , passing through disordered
anisotropic media while disorder parameter $\epsilon$ varies from 0 to 1 . The
results are for different effective beam waist in the spatial domain, (a)
2.5mm and (b) 1.88mm . The LG beam is passed through the SLM2 which is
projecting with Figure 1(a) random media. Results are shown only for RCP
polarization, was selecting using PSA.

We have investigated the role of spot size of the input beam by keeping the
$l$ value same and varying the effective spot size of the input beam. The
momentum space entropy variation shown in Fig 5. follows the same trend as in
Fig 4. A clear threshold is observed here for disordered parameter around
$\epsilon$ = 0.25, similar to the previous case. However, the variation of
momentum space entropy with the input disorder parameter $\epsilon$, is much
steeper when we have larger spotsize of the input beam. The important concept
here is that if we take the same heterogeneous system and use a smaller spot
size of the beam then it probes a lesser inhomogeneity but the corresponding
momentum span is large. On the other hand, for a larger spot size the scale of
inhomogeneity probed by the input beam is more, for a corresponding small
momentum domain span. So, one would expect a higher momentum space entropy in
the latter case. Thus, the effects of spotsize and topological charge $l$ on
the spin split modes are still combined.

Figure 5:  Momentum space entropy ($H$) variation of LG beams for different
effective beam waist in spatial domain, $\omega_{0}$=2.5mm, 1.88mm, 1.5mm,
1.25mm, 1.07mm; for (a) $l$=5 and for (b) $l$=7 for disorder parameter
$\epsilon$ varying from 0 to 1. The LG beam is passed through the SLM2 which
is projecting with Figure 1(a) random media. Results are shown only for RCP
polarization, was selecting using PSA.

To decouple their effects, we have used the Perfect Vortex (PV) beams. As
mentioned in the theory section the spatial structure of the PV beam remains
constant for all $l$ values. The experimentaly obtained PV beams are shown in
figure Fig. 6 (a). The figure shows for topological charges $l$ = 3, 5, 7 the
spatial structure of the beam remains constant in spatial domain. The
corresponding momentum space intensity distribution for $l$ = 3, 5, 7 for
input disorder parameter $\epsilon$ = 0, 0.5, 1 is shown in Fig. 6 (b).
Similar to the previous case, with increasing input disorder parameter
$\epsilon$, the beam gets more scattered. For the maximum input disorder
parameter, for $\epsilon$ = 1, the maximum intensity of the central beam
decreases as the number of random spin split modes is highest for this case.
Since spatial structure for PV beams is independent of $l$, the inhomogeneity
probed by the input PV beam remains fixed here for all values of $l$; unlike
the case of LG beam, where the inhomogeneity probed by the input beam changes
with the change of $l$ values due to change in the spatial structure of the
beam.

Figure 6: (a) Input Perfect Vortex beam in spatial domain where the spatial
structure of the beam remains same for different topological charges $l$=0, 5,
10. (b) The corresponding momentum space intensity distribution after the
input perfect vortex beam with the $l$ = 3, 5, 7 passes through inhomogeneous
anisotropic media (projected in SLM 2) with disorder parameter $\epsilon$
values = 0, 0.5, 1.

The corresponding momentum space entropies are quantified in figure 7. Even
though in the space domain, the amount of heterogeneity probed by the input PV
beam remains same for all $l$ values, their momentum space (Fourier domain
span) may get changed as PV transforms into Bessel-Gauss beam in the Fourier
domain. Despite this what we observed is with variation in values of $l$ and
with variation in disorder parameter $\epsilon$, the change of momentum space
entropy is rather minimal. Therefore, this indicates, that the influence of
topological charge $l$ alone on the statistics of momentum domain spin split
modes and on their corresponding momentum domain entropy is rather weak. For
the LG beam, the momentum domain entropy of the spin split modes is also
primarily dominated by the change in spatial structure of the beam i.e. the
effective spot size rather than its topological charge. This is further
confirmed by the simulation in Fig 8.

Figure 7: The variation in momentum space entropy as we change the input
disorder parameter ($\epsilon$, ranging from 0.2 to 1) using PV beams with
distinct topological charges ($l$ values of 0, 3, 5, 7) as input beams. The
optical medium consists of 300 azimuthal divisions, each division is assigned
a random grey value generated by the delta correlated function as mentioned
before, and these results are for right-circularly polarized (RCP) light.

Fig. 8 shows the variation of momentum space entropy with the variation input
disorder parameter for LG beam of different input spot sizes. Fig 8 (a)
corresponds to input beam waist 10 mm and Fig. 8 (b) corresponds to input beam
waist of 2.5 mm. The momentum space entropy is much higher for smaller input
spot size: 2.5 mm (fig. 8 (b)) than for the larger input spot size: 10mm (fig.
8 (a)). It is known that effective spot size of LG beam decreases with
lowering in $l$ value. Here, we clearly see a rise in threshold value as the
overall input spot size increases. As discussed before, a smaller input spot
size probes smaller area of inhomogeneity but corresponds to a larger momentum
domain span. Thus, we can conclude from fig. 8 that the role of input spatial
structure of the beam is more dominant than its topological charge ($l$).

Figure 8: Simulation showing momentum space entropy variation of the randomly
scattered spin split modes when LG beams (of different topological charges) of
different beam waists pass through the disordered inhomogeneous anisotropic
optical media. (a) beam waist = 10mm in the spatial domain. (b) beam waist =
2.5 mm (smaller spot size compared to the previous one but corresponding to a
larger Fourier domain span). As the overall input spot size increases, the
threshold value increases gradually.

##  5 Conclusion

We have investigated the role of topological charges ($l$) and the beam waist
($\omega_{0}$) in orbital angular carrying beam on the momentum space entropy
of disordered system. Our investigation shows when we are dealing with both LG
and perfect vortex beam, the spot size plays an important role in momentum
space entropy $H$ of spin-split modes compared to topological charge $l$.

## References

  * [1] K. Y. Bliokh, F. J. Rodríguez-Fortuño, F. Nori and A. V. Zayats  “Erratum: Spin–orbit interactions of light”  In _Nature Photonics_ 11.137, 2015, pp. 796–808
  * [2] Zengkai Shao et al.  “Spin-orbit interaction of light induced by transverse spin angular momentum engineering”  In _Nature communications_ 9.1 Nature Publishing Group UK London, 2018, pp. 926
  * [3] Andrea Aiello, Norbert Lindlein, Christoph Marquardt and Gerd Leuchs  “Transverse angular momentum and geometric spin Hall effect of light”  In _Physical review letters_ 103.10 APS, 2009, pp. 100401
  * [4] Xiaohui Ling et al.  “Giant photonic spin Hall effect in momentum space in a structured metamaterial with spatially varying birefringence”  In _Light: Science & Applications_ 4.e290, 2015, pp. 2047–7538
  * [5] Onur Hosten and Paul Kwiat  “Observation of the Spin Hall Effect of Light via Weak Measurements”  In _Science_ 319.787 Science, 2008
  * [6] L. B. Ma et al.  “Spin–orbit coupling of light in asymmetric microcavities”  In _Nature Communications_ 7 Nature Communications, 2016, pp. 2041–1723
  * [7] K. Yu. Bliokh, D. Yu. Frolov and Yu. A. Kravtsov  “Non-Abelian evolution of electromagnetic waves in a weakly anisotropic inhomogeneous medium”  In _Phys. Rev. A_ 75 American Physical Society, 2007, pp. 053821
  * [8] Elhanan Maguid et al.  “Disorder-induced optical transition from spin Hall to random Rashba effect”  In _Science_ 358.6369, 2017, pp. 1411–1415
  * [9] Konstantin Y Bliokh, Avi Niv, Vladimir Kleiner and Erez Hasman  “Geometrodynamics of spinning light”  In _Nature Photonics_ 2.12 Nature Publishing Group UK London, 2008, pp. 748–753
  * [10] Shiyi Xiao et al.  “Spin-dependent optics with metasurfaces”  In _Nanophotonics_ 6.1, 2017, pp. 215–234
  * [11] Xiaohui Ling et al.  “Realization of tunable photonic spin Hall effect by tailoring the Pancharatnam-Berry phase”  In _Scientific reports_ 4.1 Nature Publishing Group UK London, 2014, pp. 5557
  * [12] Konstantin Y Bliokh, Daria Smirnova and Franco Nori  “Quantum spin Hall effect of light”  In _Science_ 348.6242 American Association for the Advancement of Science, 2015, pp. 1448–1451
  * [13] Xiaobo Yin et al.  “Photonic spin Hall effect at metasurfaces”  In _Science_ 339.6126 American Association for the Advancement of Science, 2013, pp. 1405–1407
  * [14] Nir Shitrit et al.  “Spin-Optical Metamaterial Route to Spin-Controlled Photonics”  In _Science_ 340.6133, 2013, pp. 724–726
  * [15] Xiaohui Ling et al.  “Recent advances in the spin Hall effect of light”  In _Reports on Progress in Physics_ 80.6 IOP Publishing, 2017, pp. 066401
  * [16] Jian Wang and Jing Du  “Plasmonic and Dielectric Metasurfaces: Design, Fabrication and Applications”  In _Applied Sciences_ 6.9, 2016
  * [17] Jan Petersen, Jürgen Volz and Arno Rauschenbeutel  “Chiral nanophotonic waveguide interface based on spin-orbit interaction of light”  In _Science_ 346.6205 American Association for the Advancement of Science, 2014, pp. 67–71
  * [18] Andrea Aiello, Peter Banzer, Martin Neugebauer and Gerd Leuchs  “From transverse angular momentum to photonic wheels”  In _Nature Photonics_ 9.12 Nature Publishing Group UK London, 2015, pp. 789–795
  * [19] Yuri Gorodetski, Sergey Nechayev, Vladimir Kleiner and Erez Hasman  “Plasmonic Aharonov-Bohm effect: Optical spin as the magnetic flux parameter”  In _Physical Review B_ 82.12 APS, 2010, pp. 125433
  * [20] Nir Shitrit et al.  “Rashba-type plasmonic metasurface”  In _Opt. Lett._ 38.21 Optica Publishing Group, 2013, pp. 4358–4361
  * [21] Kobi Frischwasser, Igor Yulevich, Vladimir Kleiner and Erez Hasman  “Rashba-like spin degeneracy breaking in coupled thermal antenna lattices”  In _Opt. Express_ 19.23 Optica Publishing Group, 2011, pp. 23475–23482
  * [22] Mandira Pal et al.  “Tunable Spin dependent beam shift by simultaneously tailoring geometric and dynamical phases of light in inhomogeneous anisotropic medium”  In _Scientific Reports_ 6.39582, 2016, pp. 2045–2322
  * [23] Ankit Kumar Singh, Antariksha Das, Sourin Das and Nirmalya Ghosh  “Spin-selective scattering modes in a disordered anisotropic optical medium”  In _Physical Review A_ 102.3 APS, 2020, pp. 033518
  * [24] Subhasish Dutta Gupta, Nirmalya Ghosh and Ayan Banerjee  “Wave optics: Basic concepts and contemporary trends”  CRC Press, 2015
  * [25] Charles Loussert, Ulysse Delabre and Etienne Brasselet  “Manipulating the Orbital Angular Momentum of Light at the Micron Scale with Nematic Disclinations in a Liquid Crystal Film”  In _Phys. Rev. Lett._ 111 American Physical Society, 2013, pp. 037802
  * [26] Kapil Dev and Anand Asundi  “Mueller–Stokes polarimetric characterization of transmissive liquid crystal spatial light modulator” Computational Optical Measurement  In _Optics and Lasers in Engineering_ 50.4, 2012, pp. 599–607
  * [27] Vicente Durán, Jesús Lancis, Enrique Tajahuerce and Zbigniew Jaroszewicz  “Cell parameter determination of a twisted-nematic liquid crystal display by single-wavelength polarimetry”  In _Journal of Applied Physics_ 97.4, 2005, pp. 043101
  * [28] Kanghua Lu and Bahaa E. A. Saleh  “Theory and design of the liquid crystal TV as an optical spatial phase modulator”  In _Optical Engineering_ 29.3 SPIE, 1990, pp. 240 –246
  * [29] Charles Loussert, Ulysse Delabre and Etienne Brasselet  “Manipulating the orbital angular momentum of light at the micron scale with nematic disclinations in a liquid crystal film”  In _Physical Review Letters_ 111.3 APS, 2013, pp. 037802
  * [30] Yachao Liu et al.  “Generation of perfect vortex and vector beams based on Pancharatnam-Berry phase elements”  In _Scientific Reports_ 7 Scientific Reports, 2017, pp. 2045–2322
  * [31] Pravin Vaity and Leslie Rusch  “Perfect vortex beam: Fourier transformation of a Bessel beam”  In _Opt. Lett._ 40.4 Optica Publishing Group, 2015, pp. 597–600


##  1 Abstract

We examine the circular, self-similar expansion of frictional rupture due to
fluid injected at a constant rate. In a prior problem [_Viesca_ , 2021], the
author examined rupture of a planar frictional fault driven by line source of
fluid injected at constant pressure. Here, fluid injection occurs at constant
volumetric rate at a point on the fault and fluid migration occurs within a
thin, relatively permeable layer containing and parallel to the fault plane.
For the particular case when the Poisson ratio $\nu=0$, self-similarity of the
fluid pressure profile implies that fault slip will also evolve in an
axisymmetric, self-similar manner, reducing the three-dimensional problem for
the spatiotemporal evolution of fault slip to a single self-similar dimension.
The rupture radius grows as $\lambda\sqrt{4\alpha_{hy}t}$ , where $t$ is time
since the start of injection, $\alpha_{hy}$ is the hydraulic diffusivity of
the pore fluid pressure, and $\lambda$ is a prefactor determined by a single
parameter, $T$, which depends on the pre-injection stress state and injection
conditions. The prefactor has the range $0<\lambda<\infty$, where the lower
and upper limits of $\lambda$ correspond to, respectively marginal
pressurization of the fault and critically stressed conditions, in which the
fault-resolved shear stress is close to the pre-injection fault strength. In
both limits of $\lambda$, we derive asymptotic solutions for the slip by
perturbation expansion, to arbitrary order. For the critically stressed limit,
the solution for slip contains an interior boundary layer and an outer
solution. We perform a matched asymptotic expansion to provide a composite
solution with uniform convergence over the entire rupture. The leading order
of these solutions was recently used by _Saez et al. [2022]_ to verify three-
dimensional boundary-element solutions in the limit $\nu=0$.

##  2 Problem formulation

For a circular shear fracture of radius $R(t)$ undergoing relative slip
$\delta(r,t)$ in a single direction, the shear traction on the crack faces
$\tau(r,t)$, opposing the direction of slip, is given by [_Salamon and
Dundurs_ , 1971, 1977; _Viesca and Bhattacharya_ , 2019]

$\tau(r,t)=\tau_{b}+\frac{\mu}{2\pi}\int_{0}^{R(t)}\left(\frac{E[k(r/s)]}{s-r}+\frac{K[k(r/s)]}{r+s}\right)\frac{\partial\delta(s,t)}{\partial
s}ds$ (1)

where $k(u)=2\sqrt{u}/(1+u)$, and $E$ and $K$ are the complete elliptical
integrals of the second and first kind. The background shear stress resolved
on the plane of the fracture in the absence of slip, $\tau_{b}$ is here taken
to be uniform. The relation (1) presumes that the Poisson ratio $\nu=0$. This
expression is equivalent to that relating relative crack opening and crack-
face tractions for a circular mode-I crack. In that case, this expression
holds for any value of the Poisson ratio, provided that $\mu$ above is
replaced by $\mu/(1-\nu)$. The inversion of this equation for the relative
displacement $\delta$ was done in the context of mode-I equation [_Sneddon,
1951_], which we adapt here for our particular mixed-mode II/III case with
$\nu=0$

$\delta(r,t)=\frac{4}{\pi}\frac{R(t)}{\mu}\int_{\rho}^{1}\frac{\eta}{\sqrt{\eta^{2}-\rho^{2}}}\int_{0}^{1}\frac{x\,\Delta\tau(x\eta
R,t)}{\sqrt{1-x^{2}}}dx\,d\eta$ (2)

where

$\Delta\tau(r,t)=\tau_{b}-\tau(r,t)$ (3)

and $\rho=r/R(t)$.

We may rearrange this inversion to isolate the contribution of a crack-tip
singularity to the slip distribution. To do so, we first denote the innermost
integral

$f(\eta)=\int_{0}^{1}\frac{x\,\Delta\tau(x\eta R,t)}{\sqrt{1-x^{2}}}dx$

where for brevity and without loss of generality, we have suppressed explicit
mention of the time dependence of $f$. The inversion (2) reduces to

$\delta(r,t)=\frac{4}{\pi}\frac{R(t)}{\mu}\int_{\rho}^{1}\frac{\eta}{\sqrt{\eta^{2}-\rho^{2}}}f(\eta)\,d\eta$
(4)

Integration of (4) by parts yields

$\frac{\delta(r,t)}{4R(t)/(\pi\mu)}=f(1)\sqrt{1-\rho^{2}}-\int_{\rho}^{1}\sqrt{\eta^{2}-\rho^{2}}\frac{df}{d\eta}d\eta$
(5)

where we identify that

$f(1)=\int_{0}^{1}\frac{x\Delta\tau(xR,t)}{\sqrt{1-x^{2}}}dx$ (6)

is a stress intensity factor as it is a pre-factor for the relative-
displacement field behind the crack tip, $\delta\sim f(1)\sqrt{1-r}$
corresponding to a singular stress field ahead of the crack tip $\tau\sim
f(1)/\sqrt{r-1}$. This stress intensity factor is identically zero for
problems in which the fault has finite strength.

We now look to consider a particular form of $\Delta\tau$ corresponding to the
injection of fluid at a constant volumetric rate $Q$ into a thin fault zone
whose friction coefficient remains at a constant value, $f$. The strength of
the fault is given by

$\tau_{s}=f(\sigma-p(r,t))$ (7)

where $p$ is the local fluid pressure at a point on the fault, which is
presumed to have an initial, uniform value $p_{o}$ plus a perturbation due to
injection, given by the solution to the axisymmetric diffusion equation
governing fluid pressure in the plane of the thin fault zone

$\frac{\partial p}{\partial t}=\alpha_{hy}\frac{1}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial p}{\partial r}\right)$

where $\alpha_{hy}=k/(\eta_{f}\beta)$ is the hydraulic diffusivity $k$ is the
fault-zone permeability, $\eta_{f}$ is the fluid viscosity and $\beta$ is the
composite compressibility of the porous matrix and fluid. The corresponding
solution of interest is that subject to the injection condition

$-2\pi rk\frac{\partial p}{\partial r}=\frac{Q}{h}$

and is given by

$p(r,t)=p_{o}+\Delta pE_{1}[r^{2}/(\alpha t)]$ (8)

where $E_{1}$ is the elliptic integral, $\alpha=4\alpha_{hy}$, and $\Delta
p=Q/(4\pi kh)$. The thinness of the fault zone is relevant when considering
that the fault zone as a permeable, poroelastic conduit confined within an
elastic medium, which may have a stiffness dissimilar to the fault zone. Here
we consider the case where the fault-zone thickness $h$ is much smaller than
the diffusive lengthscale $\sqrt{\alpha t}$. Under this condition, there is no
change in the fault-normal stress over distances $\sim\sqrt{\alpha t}$ due to
localized, constrained swelling of the fault zone around the injection source,
and fluid pressure is decoupled from volumetric straining of the fault zone
and follows the autonomous diffusion equation [_Marck et al._ , 2015 consider
the injection scenario here; see also _Jacquey and Viesca_ , 2023, for similar
results under plane-strain conditions].

Self-similarity of the pore-fluid pressure suggests self-similarity of the
rupture front and slip profile. In this case, the rupture radius may be
written as

$R(t)=\lambda\sqrt{\alpha t}$ (9)

where $\lambda$ is a pre-factor to be determined. The non-singular stress
condition that $f(1)=0$ provides a condition to solve for the pre-factor
$\lambda$ in terms of other problem parameters. This condition is imposed as
the finite-frictional strength of the fault precludes a stress singularity. To
impose this non-singular condition, we note that the frictional strength
condition requires fault shear stress to equal fault shear strength,
$\tau(r,t)=\tau_{s}(r,t)$, when and where sliding occurs. Combining this
equation with (3), and (6)–9), and rearranging,

$T=\int_{0}^{1}\frac{x\,E_{1}\\!\left[(\lambda
x)^{2}\right]}{\sqrt{1-x^{2}}}dx$ (10)

where the dimensionless parameter

$T=\frac{\sigma^{\prime}}{\Delta p}\left(1-\frac{\tau_{b}}{\tau_{p}}\right)$
(11)

is a measure of distance of the fault stress state from critically stressed
conditions, $\tau_{b}\rightarrow\tau_{p}$, in which $T\rightarrow 0$ and the
pre-injection strength of the fault is $\tau_{p}=f\sigma^{\prime}$, and also
where $\sigma^{\prime}=\sigma-p_{o}$ is the pre-injection effective fault
normal stress. When the initial stress $\tau_{b}$ is a finite distance away
from $\tau_{p}$, but the injection-induced pressure coefficient $\Delta p$ is
small compared to $\sigma^{\prime}$, the fault is marginally pressurized and
$T\rightarrow\infty$.

We may anticipate the behavior of $\lambda$ in these two end-member regimes:
i.e., we expect that the rupture will progress quickly relatively to pore-
fluid pressure diffusion ($\lambda\gg 1$) when the fault is critically
stressed such that and relatively slowly ($\lambda\ll 1$) when marginally
pressurized. We can derive the relation between $T$ and $\lambda$ in these
limits by asymptotically expanding the integral (12).

For large values of $\lambda$, we apply the change of variable $u=x\lambda$

$T=\frac{1}{\lambda^{2}}\int_{0}^{\lambda}\frac{u\,E_{1}\\!\left[u^{2}\right]}{\sqrt{1-(u/\lambda)^{2}}}du$
(12)

followed by the Taylor expansion

$\frac{1}{\sqrt{1-(u/\lambda)^{2}}}\approx
1+\frac{(u/\lambda)^{2}}{2}+...+\frac{1}{\sqrt{\pi}}\frac{[(n-1)/2]!}{(n/2)!}(u/\lambda)^{n}$
(13)

for even $n$, and

$T\approx\sum_{n=0,2,...}^{\infty}\frac{1}{\lambda^{2+n}}\frac{1}{\sqrt{\pi}}\frac{[(n-1)/2]!}{(n/2)!}\int_{0}^{\lambda}u^{n+1}\,E_{1}\\!\left(u^{2}\right)du$
(14)

Passing the limit of the integral to infinity without loss of generality of
the power-law expansion, owing to the exponentially fast decay of the
integrand, we may first evaluate the integral by parts, second by subsequently
using another change of variable $w=u^{2}$, and third by using the definition
of the gamma function $\Gamma(z)=\int_{0}^{\infty}s^{z-1}e^{-s}ds$

$\displaystyle\begin{split}\int_{0}^{\infty}u^{n+1}\,E_{1}\\!\left(u^{2}\right)du&=\frac{2}{n+2}\int_{0}^{\infty}u^{n+1}e^{-u^{2}}du\\\\[9.0pt]
&=\frac{1}{n+2}\int_{0}^{\infty}w^{n/2}e^{-w}dw=\frac{\Gamma(n/2+1)}{n+2}=\frac{(n/2)!}{n+2}\end{split}$
(15)

Combining (14) with (15), and changing the index of summation $n=2m$, we find
the expansion of $T$ for large $\lambda$ in the critically stressed limit

$T\approx\sum_{m=0}^{\infty}\frac{1}{\sqrt{\pi}}\frac{[(2m-1)/2]!}{2m+2}\frac{1}{\lambda^{2+2m}}=\frac{1}{2\lambda^{2}}+\frac{1}{8\lambda^{4}}+\frac{1}{8\lambda^{6}}+...$
(16)

For small values of $\lambda$, we may evaluate (12) using the asymptotic
expansion of the exponential integral for small values of its argument [e.g.,
_Bender and Orszag_ , 1999]

$E_{1}\\!\left[(\lambda x)^{2}\right]\approx-\gamma-2\ln(\lambda
x)-\sum_{n=1}^{\infty}(-1)^{n}\frac{(\lambda x)^{2n}}{n\cdot n!}$ (17)

where $\gamma=0.57721...$ is the Euler-Maraschoni constant, such that

$\displaystyle\begin{split}T&\approx-\gamma\int_{0}^{1}\frac{x}{\sqrt{1-x^{2}}}dx-2\int_{0}^{1}\frac{\ln(\lambda
x)}{\sqrt{1-x^{2}}}dx-\sum_{n=1}^{\infty}(-1)^{n}\frac{\lambda^{2n}}{n\cdot
n!}\int_{0}^{1}\frac{x^{2n+1}}{\sqrt{1-x^{2}}}dx\\\\[9.0pt]
&\approx-\gamma+2-2\ln(2\lambda)-\sum_{n=1}^{\infty}(-1)^{n}\frac{\lambda^{2n}}{n\cdot
n!}\left(n!\frac{\sqrt{\pi}/2}{(n+1/2)!}\right)\\\\[9.0pt]
&\approx-\gamma+2-2\ln(2\lambda)-\sum_{n=1}^{\infty}(-1)^{n}\frac{\sqrt{\pi}/2}{n\cdot(n+1/2)!}\lambda^{2n}\\\\[9.0pt]
&\approx-\gamma+2-2\ln(2\lambda)+\frac{2}{3}\lambda^{2}-\frac{2}{15}\lambda^{4}+\frac{8}{315}\lambda^{6}-...\end{split}$
(18)

in which the first integral is readily evaluated, the second likewise but by
parts, and the evaluation of the last integral is found in Appendix C, where
it occurs for another problem and we identify it as $I_{2n+1}$.

We are interested in deriving an asymptotic expansion for slip in the limits
of large and small $\lambda$. In the marginally pressurized limit, $\lambda\ll
1$, the rupture radius lags far behind the fluid diffusive lengthscale
$\sqrt{\alpha t}$. In the critically stressed limit, where, $\lambda\gg 1$,
the rupture radius is much larger than the diffusive lengthscale,
$R(t)\gg\sqrt{\alpha t}$, and with this scale separation, there will be an
outer solution valid over distances $\sim R$ and an inner solution relevant
over distances $\sim\sqrt{\alpha t}$ near $r=0$. This inner solution is in-
effect, an interior boundary layer. In what follows, we present in section 3
the asymptotic expansion of the inner and outer solutions for the critically
stressed case, as well as their asymptotic matching and the construction of a
composite solution; and we show in section 4 the asymptotic expansion of slip
in the marginally pressurized limit in section. Details of the solution
procedure will be relegated to appendices A and B for the inner and outer
solutions, and appendix C for the marginally pressurized solution.

##  3 Critically stressed faults: inner and outer solutions

For small values of the stress parameter $T$ (large $\lambda$), the critically
stressed fault is highly sensitive to the injection of fluids and rupture
outpaces the migration of injected fluids. We perform a perturbation expansion
of the self-similar slip profile in terms of the small parameter $1/\lambda$.
Below, we provide the behavior of slip on an interior boundary layer about the
origin in section 3.1. This solution captures behavior near the origin and
will asymptotically decay with away from it. Subsequently, in section 3.2, we
provide an outer solution that captures behavior at distances greater than the
diffusive length away from the origin and up to the rupture front. The outer
solution’s approximation breaks down at the origin, where that solution is
singular. We will find that the manner of the asymptotic decay of the inner
solution matches the singularity of the outer solution and that the two
overlap over distances intermediate to the diffusive lengthscale and rupture
front position. In section 3.3, we construct a composite solution as the sum
of the inner and outer solutions, less the common overlap.

###  3.1 Inner solution

To examine the variation of slip over distances of the order of the diffusive
lengthscale $\sqrt{\alpha t}$ away from the origin, we use the scaled radial
position $s=r/\sqrt{\alpha t}$, which has the domain $0\leq s\leq\lambda$.
While the details of the solution are deferred to appendix A, the final result
is that the asymptotic expansion for slip in the large $\lambda$ limit,
$\delta_{in}$, is given by

$\displaystyle\frac{\delta_{in}(s,t)}{\sqrt{\alpha t}f\Delta
p/\mu}=\delta_{0}(s,t)+\frac{4}{\pi^{2}}\sum_{j=0}^{\infty}\frac{1}{\lambda^{2j+1}}\frac{1}{2j+1}\left[\sum_{k=0}^{j}s^{2k}\frac{\Gamma(k-1/2)}{\Gamma(k+1)}\frac{\Gamma(3/2+j-k)}{2(j-k)+1}\right]$

where the leading, zeroth-order term is

$\delta_{0}(s,t)=2\sqrt{\pi}e^{-s^{2}/2}\left[(1+s^{2})I_{0}(s^{2}/2)+s^{2}I_{1}(s^{2}/2)\right]-4s$

and $I_{0}$ and $I_{1}$ are modified Bessel functions of the first kind. The
subsequent, higher-order terms are polynomial and, to fifth-order in the small
parameter $1/\lambda$, the expansion is

$\frac{\delta_{in}(s,t)}{\sqrt{\alpha t}f\Delta
p/\mu}=\delta_{0}(s,t)-\frac{1}{\lambda}\frac{4}{\pi}-\frac{1}{\lambda^{3}}\frac{2}{3\pi}(1-s^{2})-\frac{1}{\lambda^{5}}\frac{1}{5\pi}\left(3-s^{2}-s^{4}/2\right)+O(\lambda^{-7})$
(19)

###  3.2 Outer solution

We now examine behavior over radial distances comparable to the current
rupture radius $R(t)$, relatively much larger than the diffusive lengthscale
$\sqrt{\alpha t}$. We rescale radial distance as $\rho=r/R(t)$, and find the
asymptotic expansion of the outer solution for slip to be

$\displaystyle\begin{split}\frac{\delta_{out}(\rho,t)}{Rf\Delta
p/\mu}=\frac{4}{\pi^{2}}\sum_{n=0}^{\infty}&\frac{1}{\lambda^{2n+2}}\frac{1}{\rho^{2n+1}}\frac{\Gamma(n+3/2)}{2n+1}\frac{\Gamma(n+1/2)}{\Gamma(n+2)}\\\\[9.0pt]
&\cdot\left[\left(\text{acos}(\rho)-\rho\sqrt{1-\rho^{2}}\right)+2\sqrt{\pi}\frac{(1-\rho^{2})^{3/2}}{\rho}\sum_{k=1}^{n}\frac{\rho^{2k}}{2k+2}\frac{\Gamma(k+2)}{\Gamma(k+1/2)}\right]\end{split}$
(20)

where we have deferred discussion of the solution details to appendix B. This
expansion, to sixth-order in the small parameter $1/\lambda$ when scaled as
such, reduces to

$\displaystyle\frac{\delta_{out}(\rho,t)}{Rf\Delta p/\mu}=$
$\displaystyle\frac{1}{\lambda^{2}}\frac{2}{\pi}\left(\frac{\text{acos}(\rho)}{\rho}-\sqrt{1-\rho^{2}}\right)+\frac{1}{\lambda^{4}}\frac{1}{4\pi}\left[\frac{\text{acos}(\rho)}{\rho^{3}}-\sqrt{1-\rho^{2}}\left(2-\frac{1}{\rho^{2}}\right)\right]$
$\displaystyle+\frac{1}{\lambda^{6}}\frac{3}{16\pi}\left[\frac{\text{acos}(\rho)}{\rho^{5}}-\sqrt{1-\rho^{2}}\left(\frac{8}{3}-\frac{2}{3}\frac{1}{\rho^{2}}-\frac{1}{\rho^{4}}\right)\right]+O(\lambda^{-8})$

###  3.3 Matched asymptotic expansion and composite solution

The inner and outer expansions can be used to construct a composite asymptotic
solution, which will provide uniform convergence over the entire domain [e.g.,
_Hinch_ , 1991]. The composite solution is given by

$\frac{\delta_{comp}(\rho,t)}{\sqrt{\alpha t}f\Delta
p/\mu}=\frac{\delta_{in}(\lambda\rho,t)+\delta_{out}(\rho,t)-\delta_{overlap}(\rho,t)}{\sqrt{\alpha
t}f\Delta p/\mu}$ (21)

where $\delta_{overlap}$ is the common asymptotic behavior of the inner and
outer solutions in an intermediate region corresponding to the outer region
($s\rightarrow\infty$) of the inner solution, and the inner region
($\rho\rightarrow 0$) of the outer solution. This overlap is subtracted to
avoid double-counting contributions in this intermediate region when summing
the inner and outer solutions.

The expression for the inner solution is given as before,

$\displaystyle\frac{\delta_{in}(s,t)}{\sqrt{\alpha t}f\Delta p/\mu}=\,\,$
$\displaystyle
2\sqrt{\pi}e^{-s^{2}/2}\left[(1+s^{2})I_{0}(s^{2}/2)+s^{2}I_{1}(s^{2}/2)\right]-4s$
$\displaystyle-\frac{4}{\pi\lambda}-\frac{2}{3\pi\lambda^{3}}\left(1-s^{2}\right)-\frac{1}{5\pi\lambda^{5}}\left(3-s^{2}-\frac{s^{4}}{2}\right)+O(\lambda^{-7})$

Regarding the outer solution, we must carefully note that the above scaling
for the composite and inner solutions do not match that of the prior
expression of the outer solution (20): their scalings differ by a factor of
$R(t)/\sqrt{\alpha t}=\lambda$. Rescaling the outer solution in the manner of
(21),

$\displaystyle\frac{\delta_{out}(\rho,t)}{\sqrt{\alpha t}f\Delta p/\mu}=$
$\displaystyle\frac{1}{\lambda}\frac{2}{\pi}\left(\frac{\text{acos}(\rho)}{\rho}-\sqrt{1-\rho^{2}}\right)+\frac{1}{\lambda^{3}}\frac{1}{4\pi}\left[\frac{\text{acos}(\rho)}{\rho^{3}}-\sqrt{1-\rho^{2}}\left(2-\frac{1}{\rho^{2}}\right)\right]$
$\displaystyle+\frac{1}{\lambda^{5}}\frac{3}{16\pi}\left[\frac{\text{acos}(\rho)}{\rho^{5}}-\sqrt{1-\rho^{2}}\left(\frac{8}{3}-\frac{2}{3}\frac{1}{\rho^{2}}-\frac{1}{\rho^{4}}\right)\right]+O(\lambda^{-8})$

we may now evaluate the so-called overlap of the inner and outer solutions,
here both expanded to fifth-order in $1/\lambda$ when scaled as such.

The overlap of these two solutions, $\delta_{overlap}$, is determined by
comparing the outer behavior of inner solution, i.e., the asymptotic behavior
of $\delta_{in}(s\rightarrow\infty,t)$ with the inner behavior of the outer
solution, i.e., the asymptotic behavior of $\delta_{out}(\rho\rightarrow
0$,t). The latter can be determined to be

$\displaystyle\frac{\delta_{out}(\rho\rightarrow 0,t)}{\sqrt{\alpha t}f\Delta
p/\mu}$
$\displaystyle=\frac{1}{\lambda}\left({\color[rgb]{1,0,0}\frac{1}{\rho}-\frac{4}{\pi}}+{\color[rgb]{0,0,1}\frac{2}{3}\frac{\rho^{2}}{\pi}}+{\color[rgb]{0,1,0}\frac{1}{10\pi}\rho^{4}}+O(\rho^{6})\right)$
$\displaystyle+\frac{1}{\lambda^{3}}\left({\color[rgb]{0,0,1}\frac{1}{8\rho^{3}}-\frac{2}{3\pi}}+{\color[rgb]{0,1,0}\frac{\rho^{2}}{5\pi}}+\frac{\rho^{4}}{28\pi}+O(\rho^{6})\right)$
$\displaystyle+\frac{1}{\lambda^{5}}\left({\color[rgb]{0,1,0}\frac{3}{32\rho^{5}}-\frac{3}{5\pi}}+\frac{3\rho^{2}}{14\pi}+\frac{\rho^{4}}{24\pi}+O(\rho^{6})\right)+O(\lambda^{-7})$

where we again note here that the scaling of $\delta_{out}$ on the left hand
side is different from the scaling of slip in (20), in which
$R(t)=\lambda\sqrt{\alpha t}$ was used in place of $\sqrt{\alpha t}$. In the
equation above, we scale slip with $\sqrt{\alpha t}$ so that we may correctly
compare this inner behavior with the outer behavior of the inner solution,
which is

$\frac{\delta_{in}(s\rightarrow\infty,t)}{\sqrt{\alpha t}f\Delta
p/\mu}=\left({\color[rgb]{1,0,0}\frac{1}{s}}+{\color[rgb]{0,0,1}\frac{1}{8s^{3}}}+{\color[rgb]{0,1,0}\frac{3}{32s^{5}}}+O(s^{-7})\right){\color[rgb]{1,0,0}-\frac{1}{\lambda}\left(\frac{4}{\pi}\right)}{\color[rgb]{0,0,1}-\frac{1}{\lambda^{3}}\left(\frac{2}{3\pi}(1-s^{2})\right)}{\color[rgb]{0,1,0}-\frac{1}{\lambda^{5}}\left(\frac{1}{5\pi}(3-s^{2}-s^{4}/2)\right)}$

Recognizing that $s=\rho\lambda$, we can identify the common, overlapping
terms for the two expansions, which we color in red, blue, and green to
facilitate comparisons. This overlap is clearly

$\displaystyle\begin{split}\frac{\delta_{overlap}(\rho,t)}{\sqrt{\alpha
t}f\Delta
p/\mu}&={\color[rgb]{1,0,0}\frac{1}{\lambda}\left(\frac{1}{\rho}-\frac{4}{\pi}\right)}+{\color[rgb]{0,0,1}\frac{1}{\lambda^{3}}\left[\frac{1}{8\rho^{3}}-\frac{2}{3\pi}\left(1-(\rho\lambda^{2})\right)\right]}\\\\[9.0pt]
&+{\color[rgb]{0,1,0}\frac{1}{\lambda^{5}}\left[\frac{3}{32\rho^{5}}-\frac{1}{5\pi}\left(3-(\rho\lambda)^{2}-\frac{(\rho\lambda)^{4}}{2}\right)\right]}+O(\lambda^{-7})\end{split}$
(22)

For example, the composite solution to fifth order in the small parameter
$1/\lambda$ is

$\displaystyle\begin{split}\frac{\delta_{comp}(\rho,t)}{\sqrt{\alpha t}f\Delta
p/\mu}&=2\sqrt{\pi}e^{-(\lambda\rho)^{2}/2}\left\\{[1+(\lambda\rho)^{2}]I_{0}\left[\frac{(\lambda\rho)^{2}}{2}\right]+(\lambda\rho)^{2}I_{1}\left[\frac{(\lambda\rho)^{2}}{2}\right]\right\\}-4\lambda\rho\\\\[9.0pt]
&+\frac{1}{\lambda}\frac{2}{\pi}\left(\frac{\text{acos}(\rho)}{\rho}-\sqrt{1-\rho^{2}}\right)-\frac{1}{\lambda\rho}\\\\[9.0pt]
&+\frac{1}{\lambda^{3}}\frac{1}{4\pi}\left[\frac{\text{acos}(\rho)}{\rho^{3}}-\sqrt{1-\rho^{2}}\left(2-\frac{1}{\rho^{2}}\right)\right]-\frac{1}{(2\lambda\rho)^{3}}\\\\[9.0pt]
&+\frac{1}{\lambda^{5}}\frac{3}{16\pi}\left[\frac{\text{acos}(\rho)}{\rho^{5}}-\sqrt{1-\rho^{2}}\left(\frac{8}{3}-\frac{2}{3}\frac{1}{\rho^{2}}-\frac{1}{\rho^{4}}\right)\right]-\frac{3}{(2\lambda\rho)^{5}}\end{split}$
(23)

where the first two rows are the composite to first order and the first three
are that to third order.

##  4 Fault slip under marginally pressurized conditions

In the limit of large $T$ (small $\lambda$), the self-similar fluid pressure
profile stretches faster than crack growth: rupture lags the injected fluid.
In this case, we may perform a perturbation expansion for the slip profile in
powers of the small parameter $\lambda$. The details of the solution may be
found in appendix C, but the perturbation to arbitrary order is given by

$\frac{\delta_{mp}(\rho,t)}{R(t)f\Delta
p/\mu}=\frac{8}{\pi}\left(\sqrt{1-\rho^{2}}-\rho\text{
acos}(\rho)\right)+2\left(1-\rho^{2}\right)^{3/2}\sum_{n=1}^{\infty}\frac{\lambda^{2n}}{2n+1}\frac{1}{\Gamma(n+3/2)}\sum_{k=1}^{n}\frac{\Gamma(k)}{\Gamma(k+1/2)}\rho^{2k-2}$

To fourth order in $\lambda$, the asymptotic solution is

$\frac{\delta_{mp}(\rho,t)}{R(t)f\Delta
p/\mu}=\frac{8}{\pi}\left(\sqrt{1-\rho^{2}}-\rho\text{
acos}(\rho)\right)+\lambda^{2}\frac{16}{9\pi}(1-\rho^{2})^{3/2}+\lambda^{4}\frac{32}{225\pi}(1-\rho^{2})^{3/2}(3+2\rho^{2})+O(\lambda^{6})$

References

Bhattacharya, P., and R. C. Viesca (2019) Fluid-induced aseismic fault slip
outpaces pore-fluid migration, _Science_ , 364, 464–468,
doi:10.1126/science.aaw7354

Bender, C. M. and S. A. Orszag (1999) Advanced Mathematical Methods for
Scientists and Engineers, _Springer-Verlag_ , doi:10.1007/978-1-4757-3069-2

Hinch, E. J. (1991), Perturbation Methods, _Cambridge University Press_ ,
doi:10.1017/CBO9781139172189

Jacquey, A. B. and R. C. Viesca (2023) Nucleation and arrest of fluid-induced
aseismic slip, _Geophys. Res. Lett._ , 50, e2022GL101228,
doi:10.1029/2022GL101228

Marck, J., A. A. Savitski, and E. Detournay (2015), Line source in a
poroelasic layer bounded by an elastic space, _Int. J. Numer. Anal. Meth.
Geomech_ , 39, 1484–1505, doi:10.1002/nag.2405

Salamon and Dundurs (1971) Elastic fields of a dislocation loop in a two-phase
material, _J. Elasticity_ , 1(2), 153–164, doi:10.1007/BF00046466

Salamon and Dundurs (1977) A circular glide dislocation loop in a two-phase
material, _J. Phys. C Solid State_ , 10, 497–507,
doi:10.1088/0022-3719/10/4/007

Sáez, A., B. Lecampion, P. Bhattacharya, R. C. Viesca (2022) Three-dimensional
fluid-drive stable frictional ruptures, _J. Mech. Phys. Solids._ , 160,
104754, doi:10.1016/j.jmps.2021.104754

Sneddon, I. N. (1951) Fourier Transforms, _McGraw-Hill_.

Appendices

##  Appendix A Asymptotic solution details of inner solution for critically
stressed conditions

Here we derive the boundary-layer solution for slip near $r=0$ valid over
distances comparable to the diffusion lengthscale $\sqrt{\alpha t}$. However,
we note that in our original inversion for slip (2) and (4), radial distances
are scaled by the current rupture radius $R$. Here we revisit the problem of
solving for slip, but now scale radial distances by $\sqrt{\alpha t}$,
denoting

$s=r/\sqrt{\alpha t}$

Substituting the variable of integration in (4) with $\eta=\xi/\lambda$, we
may thus rewrite (4) as

$\delta(s,t)=-\frac{4}{\pi}\frac{\sqrt{\alpha
t}}{\mu}\int_{s}^{\lambda}\sqrt{\xi^{2}-s^{2}}\frac{df}{d\xi}d\xi$ (24)

where

$\frac{df}{d\xi}=\frac{d}{d\xi}\int_{0}^{1}\frac{x\Delta\tau(x\xi\sqrt{\alpha
t})}{\sqrt{1-x^{2}}}dx$ (25)

in which we have suppressed explicit time dependence of $\Delta\tau$, except
through the term $\sqrt{\alpha t}$. That we may do so is apparent considering
that the particular expression for $\Delta\tau$ of interest follows the
combination of (3) and (6)–9) and the requirement that $\tau=\tau_{s}$,

$\Delta\tau(r,t)=\tau_{b}-f\left[\sigma^{\prime}-\Delta pE_{1}(r^{2}/\alpha
t)\right]$

evaluating

$\Delta\tau(x\xi\sqrt{\alpha t},t)=\tau_{b}-f\left(\sigma^{\prime}-\Delta
pE_{1}[(x\xi)^{2}]\right)$

In this case, we can simplify $\Delta\tau(x\xi\sqrt{\alpha
t},t)\Rightarrow\Delta\tau(x\xi)$ and do so in our evaluation of (25),

$\displaystyle\frac{df}{d\xi}$
$\displaystyle=\frac{d}{d\xi}\int_{0}^{1}\frac{x\Delta\tau(x\xi)}{\sqrt{1-x^{2}}}dx$
$\displaystyle=\int_{0}^{1}\frac{x}{\sqrt{1-x^{2}}}\displaystyle\frac{d\Delta\tau}{dw}{\displaystyle\frac{dw}{d\xi}}dx$
$\displaystyle=-2f\Delta
p\,\frac{1}{\xi}\int_{0}^{1}\frac{xe^{-(x\xi)^{2}}}{\sqrt{1-x^{2}}}dx$

where in the above $w=x\xi$ and $d\Delta\tau(w)/dw=2\exp[-(x\xi)^{2}]/(x\xi)$.
We now label the integral

$F(\xi)=\int_{0}^{1}\frac{xe^{-(x\xi)^{2}}}{\sqrt{1-x^{2}}}dx$

such that (24) becomes

$\delta(s,t)=\frac{\sqrt{\alpha t}}{\mu}f\Delta
p\frac{8}{\pi}\int_{s}^{\lambda}\sqrt{\xi^{2}-s^{2}}\frac{F(\xi)}{\xi}d\xi$

or more simply

$\frac{\delta(s,t)}{\sqrt{\alpha t}f\Delta
p/\mu}=\frac{8}{\pi}\int_{s}^{\lambda}\sqrt{1-(s/\xi)^{2}}\,F(\xi)\,d\xi$ (26)

We now look to asymptotically evaluate the expression (26), which is exact, in
the limit of large $\lambda$ to arrive to our expansion approximation of the
inner solution for critically stressed conditions. To do so, we will first
rewrite (26) as:

$\frac{\delta(s,t)}{\sqrt{\alpha t}f\Delta
p/\mu}=\frac{8}{\pi}\left[\int_{s}^{\infty}\sqrt{1-(s/\xi)^{2}}F(\xi)-\int_{\lambda}^{\infty}\sqrt{1-(s/\xi)^{2}}F(\xi)\right]$
(27)

where the leading-order term, $O(1)$, is

$\frac{8}{\pi}\int_{s}^{\infty}\sqrt{1-(s/\xi)^{2}}F(\xi)d\xi=2\sqrt{\pi}e^{-s^{2}/2}\left[(1+s^{2})I_{0}(s^{2}/2)+s^{2}I_{1}(s^{2}/2\right]-4s$
(28)

where $I_{0}$ and $I_{1}$ modified Bessel functions of the first kind. We turn
our attention to the second integral in (27) and note that, for large
$\lambda$, the integration variable $\xi$ takes on only large values
warranting a Taylor expansion of the integrand’s square-root term

$-\int_{s}^{\infty}\sqrt{1-(s/\xi)^{2}}F(\xi)d\xi\approx-\int_{s}^{\infty}\left(1-\frac{1}{2}(s/\xi)^{2}-\frac{1}{8}(s/\xi)^{4}+...\right)F(\xi)d\xi$
(29)

We look to evaluate the series of integrals of (29), which are the remaining
$O(1/\lambda)$ terms of (27). We start with the first integral

$I(\lambda)=-\int_{\lambda}^{\infty}F(\xi)d\xi$

To expand $I$ for large $\lambda$, we note that

$\frac{dI}{d\lambda}=F(\lambda)$

and seek the expansion of

$F(\lambda)=\int_{0}^{1}\frac{xe^{-(\lambda x)^{2}}}{\sqrt{1-x^{2}}}dx$ (30)

Since $(\lambda x)^{2}$ has a minimum at $x=0$, we use Laplace’s method for
asymptotic expansion, denoting $g(x)=x/\sqrt{1-x^{2}}$:

$F(\lambda)\approx\int_{0}^{\infty}\left[g(0)+g^{\prime}(0)x+\frac{g^{\prime\prime}(0)}{2}x^{2}+\frac{g^{\prime\prime\prime}(0)}{6}x^{3}+...\right]e^{-(\lambda
x)^{2}}dx$

with $g(0)=0$, $g^{\prime}(0)=1$, $g^{\prime\prime}(0)=0$,
$g^{\prime\prime\prime}(0)=1/2$, etc., or more compactly written as

$F(\lambda)\approx\int_{0}^{\infty}\left[\sum_{m=0}^{\infty}\begin{cases}\frac{2^{m}}{\pi}\frac{\Gamma(1+m/2)\Gamma(m/2)}{m!}x^{m}&m\text{
odd}\\\ 0&m\text{ even}\end{cases}\right]e^{-(\lambda x)^{2}}dx$ (31)

This requires the evaluation of

$\int_{0}^{\infty}x^{m}e^{-(\lambda x)^{2}}dx$

which can be done with the substitution $u=(\lambda x)^{2}$

$\frac{1}{2\lambda^{m+1}}\int_{0}^{\infty}u^{\frac{m}{2}-1}e^{-u}du=\frac{1}{2\lambda^{m+1}}\Gamma\left(\frac{m+1}{2}\right)$
(32)

where we recall the definition of the gamma function
$\Gamma(z)=\int_{0}^{\infty}u^{z-1}e^{-z}dz=(z-1)!$. Combining (31) and (32)

$F(\lambda)=\sum_{m=1,3,...}^{\infty}\frac{1}{2\lambda^{m+1}}\frac{2^{m}}{\pi}\frac{\Gamma(1+m/2)\Gamma(m/2)\Gamma(m/2+1/2)}{m!}$

which can be further simplified using the Legendre duplication formula
$\Gamma(m)=2^{m-1}\Gamma(m/2)\Gamma(m/2+1/2)/\sqrt{\pi}$

$F(\lambda)=\sum_{m=1,3,...}^{\infty}\frac{1}{\sqrt{\pi}}\frac{1}{\lambda^{m+1}}\frac{\Gamma(1+m/2)\Gamma(m)}{m!}$

Recognizing that $\Gamma(m)=(m-1)!$

$F(\lambda)=\sum_{m=1,3,...}^{\infty}\frac{1}{\sqrt{\pi}}\frac{1}{\lambda^{m+1}}\frac{\Gamma(1+m/2)}{m}$

and changing the summation index $m=2n+1$ with $n=0,1,2,...$

$F(\lambda)=\sum_{n=0}^{\infty}\frac{1}{\sqrt{\pi}}\frac{1}{\lambda^{2n+2}}\frac{\Gamma(n+3/2)}{2n+1}$
(33)

we may finally calculate the integral

$I(\lambda)=-\int_{\lambda}^{\infty}F(\lambda)d\lambda=-\frac{1}{\sqrt{\pi}}\sum_{n=0}^{\infty}\frac{1}{\lambda^{2n+1}}\frac{\Gamma(n+3/2)}{(2n+1)^{2}}$

We now look to evaluate all terms in the expansion (29), which can be written
in summation form as

$-\int_{\lambda}^{\infty}\sqrt{1-(s/\xi)^{2}}F(\xi)d\xi\approx-\int_{\lambda}^{\infty}\left[-\sum_{k=0}^{\infty}\frac{(k-3/2)!}{2\sqrt{\pi}k!}\left(\frac{s}{\xi}\right)^{2k}\right]F(\xi)d\xi$

where the summation in brackets is the Taylor series of $\sqrt{1-(s/\xi)^{2}}$
about $s/\xi=0$. Passing factorials to Gamma functions, via the identity
$\Gamma(n+1)=n!$, and rearranging we have

$-\int_{\lambda}^{\infty}\sqrt{1-(s/\xi)^{2}}F(\xi)d\xi=\sum_{k=0}^{\infty}\frac{1}{2\sqrt{\pi}}\frac{\Gamma(k-1/2)}{\Gamma(k+1)}s^{2k}\int_{\lambda}^{\infty}\frac{F(\xi)}{\xi^{2k}}d\xi$
(34)

we denote the integral

$J(\lambda)=\int_{\lambda}^{\infty}\frac{F(\xi)}{\xi^{2k}}d\xi$

such that

$\frac{dJ}{d\lambda}=-\frac{F(\lambda)}{\lambda^{2k}}=-\frac{1}{\sqrt{\pi}}\frac{1}{\lambda^{2k}}\left(\sum_{m=1,3,...}^{\infty}\frac{1}{\lambda^{2n+2}}\frac{\Gamma(n+3/2)}{2n+1}\right)$

in which we have substituted (33) and may subsequently integrate to find

$J(\lambda)=-\int_{\lambda}^{\infty}\frac{dJ}{d\lambda}d\lambda=\frac{1}{\sqrt{\pi}}\sum_{n=0}^{\infty}\frac{1}{\lambda^{2(n+k)+1}}\frac{1}{2(n+k)+1}\frac{\Gamma(n+3/2)}{2n+1}$

Substituting this expression for $J$ into (34)

$\displaystyle-\int_{\lambda}^{\infty}\sqrt{1-(s/\xi)^{2}}F(\xi)d\xi$
$\displaystyle=\frac{1}{2\pi}\sum_{k=0}^{\infty}\sum_{n=0}^{\infty}s^{2k}\frac{\Gamma(k-1/2)}{\Gamma(k+1)}\frac{1}{\lambda^{2(n+k)+1}}\frac{1}{2(n+k)+1}\frac{\Gamma(n+3/2)}{2n+1}$
(35)
$\displaystyle=\frac{1}{2\pi}\sum_{k=0}^{\infty}\sum_{j=k}^{\infty}s^{2k}\frac{\Gamma(k-1/2)}{\Gamma(k+1)}\frac{1}{\lambda^{2j+1}}\frac{1}{2j+1}\frac{\Gamma(3/2+j-k)}{2(j-k)+1}$
(36)
$\displaystyle=\frac{1}{2\pi}\sum_{j=0}^{\infty}\frac{1}{\lambda^{2j+1}}\frac{1}{2j+1}\left[\sum_{k=0}^{j}s^{2k}\frac{\Gamma(k-1/2)}{\Gamma(k+1)}\frac{\Gamma(3/2+j-k)}{2(j-k)+1}\right]$
(37)

where, in passing from (35) to (36) we replaced the summation index $n$ with
another $j$ chosen such that $2j+1=2(n+k)+1$, or $n=j-k$, and in passing from
(36) to (37) we swapped the order of summation. We thus find that remaining
term in (27) is asymptotically expanded as

$\frac{8}{\pi}\left(-\int_{\lambda}^{\infty}\sqrt{1-(s/\xi)^{2}}F(\xi)d\xi\right)\approx-\frac{1}{\lambda}\frac{4}{\pi}-\frac{1}{\lambda^{3}}\frac{2}{3\pi}(1-s^{2})-\frac{1}{\lambda^{5}}\frac{1}{5\pi}\left(3-s^{2}-s^{4}/2\right)+O(\lambda^{-7})$
(38)

We therefore arrive to the asymptotic expansion of slip of the inner solution
for slip in the critically stressed limit:

$\frac{\delta_{in}(s,t)}{\sqrt{\alpha t}f\Delta
p/\mu}=\delta_{0}(s,t)-\frac{1}{\lambda}\frac{4}{\pi}-\frac{1}{\lambda^{3}}\frac{2}{3\pi}(1-s^{2})-\frac{1}{\lambda^{5}}\frac{1}{5\pi}\left(3-s^{2}-s^{4}/2\right)+O(\lambda^{-7})$
(39)

where

$\delta_{0}(s,t)=2\sqrt{\pi}e^{-s^{2}/2}\left[(1+s)I_{0}(s^{2}/2)+s^{2}I_{1}(s^{2}/2\right]-4s$

the outer behavior of which as $s\rightarrow\infty$ is

$\frac{\delta_{in}(s\rightarrow\infty,t)}{\sqrt{\alpha t}f\Delta
p/\mu}=\left(\frac{1}{s}+\frac{1}{8s^{3}}+\frac{3}{32s^{5}}+O(s^{-7})\right)-\frac{1}{\lambda}\frac{4}{\pi}-\frac{1}{\lambda^{3}}\frac{2}{3\pi}(1-s^{2})-\frac{1}{\lambda^{5}}\frac{1}{5\pi}\left(3-s^{2}-s^{4}/2\right)+O(\lambda^{-7})$
(40)

##  Appendix B Asymptotic solution details of outer solution for critically
stressed conditions

Here we turn our attention to the solution for slip for distances $r\sim R(t)$
when $R(t)\gg\sqrt{\alpha t}$, i.e., when $\lambda$ is large. We return to the
inversion for slip (5), applying the non-singular condition $f(1)=0$:

$\delta(r,t)=-\frac{4}{\pi}\frac{R}{\mu}\int_{\rho}^{1}\sqrt{\mu^{2}-\rho^{2}}\frac{d}{d\eta}\int_{0}^{1}\frac{x\Delta\tau(x\eta
R,t)}{\sqrt{1-x^{2}}}dx$ (41)

where we recall that

$\Delta\tau(r,t)=\tau_{b}-f[\sigma^{\prime}-\Delta pE_{1}(r^{2}/\alpha t)]$

When evaluating

$\Delta\tau(x\eta R(t),t)=\tau_{b}-f\left(\sigma^{\prime}-\Delta
pE_{1}[(x\eta\lambda)^{2}]\right)$

we note that we can simplify the left hand side $\Delta\tau(x\xi\sqrt{\alpha
t},t)\Rightarrow\Delta\tau(x\eta\lambda)$ and do so in our evaluation of (41),
as was done similarly for the inner solution in Appendix A:

$\displaystyle\frac{d}{d\eta}\int_{0}^{1}\frac{x\Delta\tau(x\eta\lambda)}{\sqrt{1-x^{2}}}dx$
$\displaystyle=\int_{0}^{1}\frac{x}{\sqrt{1-x^{2}}}\frac{d\Delta\tau}{dw}\frac{dw}{d\zeta}dx$
$\displaystyle=-2f\Delta
p\frac{1}{\eta}\int_{0}^{1}\frac{xe^{-(x\eta\lambda)^{2}}}{\sqrt{1-x^{2}}}dx$

where now $w=x\eta\lambda$, and (41) becomes

$\delta(r,t)=\frac{Rf\Delta
p}{\mu}\frac{8}{\pi}\int_{\rho}^{1}\frac{\sqrt{\eta^{2}-\rho^{2}}}{\eta}\int_{0}^{1}\frac{xe^{-(x\eta\lambda)^{2}}}{\sqrt{1-x^{2}}}dxd\eta$
(42)

We identify the innermost integral, following (30), as $F(\eta\lambda)$ which
has the known asymptotic expansion for large values of its argument (33). We
substitute this expansion in (42) to find the asymptotic expansion of the
outer solution, $\delta_{out}$

$\delta_{out}(\rho,t)=\frac{Rf\Delta
p}{\mu}\frac{8}{\pi}\sum_{n=0}^{\infty}\frac{1}{\lambda^{2n+2}}\frac{1}{\sqrt{\pi}}\frac{\Gamma(n+3/2)}{2n+1}\int_{\rho}^{1}\frac{\sqrt{\eta^{2}-\rho^{2}}}{\eta^{2n+3}}d\eta$

To evaluate the integral, we use the trigonemetric variable substitution
$\cos\theta=\rho/\eta$ on

$\int_{\rho}^{1}\frac{\sqrt{\eta^{2}-\rho^{2}}}{\eta^{2n+3}}d\eta=\frac{1}{\rho^{2n+1}}\int_{0}^{\text{acos}(\rho)}\cos^{2n}\theta\sin^{2}\theta
d\theta$

and to evaluate this last integral in turn, we’ll consider the general
integral

$I(m,p,\rho)=\int_{0}^{\text{acos}(\rho)}\cos^{m}\theta\sin^{p}\theta d\theta$

whose solution is as follows. We begin by rearranging the integrand

$\displaystyle I(m,p,\rho)$
$\displaystyle=\int_{0}^{\text{acos}(\rho)}\left(\cos^{m-1}\theta\right)\left(\sin^{p}\theta\cos\theta\right)d\theta$
$\displaystyle=\frac{\rho^{m-1}(1-\rho^{2})^{\frac{p+1}{2}}}{p+1}+\frac{m-1}{p-1}\int_{0}^{\text{acos}(\rho)}\cos^{m-2}\theta\sin^{p+2}d\theta$
$\displaystyle=\frac{\rho^{m-1}(1-\rho^{2})^{\frac{p+1}{2}}}{p+1}-\frac{m-1}{p+1}\int_{0}^{\text{acos}(\rho)}\cos^{m}\theta\sin^{p}\theta
d\theta+\frac{m-1}{p+1}\int_{0}^{\text{acos}(\rho)}\cos^{m-2}\theta\sin^{p}\theta
d\theta$

and, in the above, we pass from the first to second line using integration by
parts, and from the second to third via $\sin^{2}\theta=1-\cos^{2}\theta$. We
identify the first integral as $I(m,p,\rho)$ and the second as
$I(m-2,p,\rho)$. Solving for $I(m,p,\rho)$, we find the recursion relation:

$I(m,p,\rho)=\frac{\rho^{m-1}(1-\rho^{2})^{\frac{p+1}{2}}}{m+p}+\frac{m-1}{m+p}I(m-2,p,\rho)$
(43)

Returning to the specific case $m=2n$ and $p=2$, we find the recursion
relation for which we need only calculate the case $n=0$ to find all others
$n>0$:

$\displaystyle I(0,2,\rho)$
$\displaystyle=\frac{1}{2}\left(\text{acos}(\rho)-\rho\sqrt{1-\rho^{2}}\right)$
$\displaystyle I(2,2,\rho)$
$\displaystyle=\frac{1}{8}\left(\text{acos}(\rho)-\rho\sqrt{1-\rho^{2}}\right)+\frac{\rho(1-\rho^{2})^{3/2}}{4}$
$\displaystyle=\frac{1}{8}\left[\left(\text{acos}(\rho)-\rho\sqrt{1-\rho^{2}}\right)+\frac{(1-\rho^{2})^{3/2}}{\rho}2\rho^{2}\right]$
$\displaystyle I(4,2,\rho)$
$\displaystyle=\frac{1}{16}\left(\text{acos}(\rho)-\rho\sqrt{1-\rho^{2}}\right)+\frac{\rho(1-\rho^{2})^{3/2}}{8}+\frac{\rho^{3}(1-\rho^{2})^{3/2}}{6}$
$\displaystyle=\frac{1}{16}\left[\left(\text{acos}(\rho)-\rho\sqrt{1-\rho^{2}}\right)+\frac{(1-\rho^{2})^{3/2}}{\rho}\left(2\rho^{2}+\frac{8}{3}\rho^{4}\right)\right]$

From this recursion, we may deduce an explicit expression for $I(m,2,\rho)$,
recognizing from the sequence above that such an expression will have the form

$I(m,2,\rho)=a_{m}\left[\left(\text{acos}(\rho)-\rho\sqrt{1-\rho^{2}}\right)+\frac{(1-\rho^{2})^{3/2}}{\rho}\sum_{k=2,4,...}^{m}b_{k}\rho^{k}\right]$

where $a_{m}$ and $b_{k}$ are to be determined. From the factor $(m-1)/(m+2)$
in the recursion, we recognize that

$\displaystyle
a_{m}=\frac{(m-1)}{(m+2)}\cdot\frac{(m-3)}{m}\cdot\frac{(m-5)}{(m-2)}\cdot...\cdot\frac{7}{10}\cdot\frac{5}{8}\cdot\frac{3}{6}\cdot\frac{1}{4}\cdot\frac{1}{2}$
(44)

Substituting $m=2n$

$\displaystyle a_{2n}$
$\displaystyle=\frac{(2n-1)(2n-3)(2n-5)...}{(2n+2)(2n)(2n-2)...}$ (45)
$\displaystyle=\frac{(2n-1)(2n-3)(2n-5)...}{2^{n+1}(n+1)(n)(n-1)...}$ (46)
$\displaystyle=\frac{(2n-1)!/[2^{n-1}(n-1)!]}{2^{n+1}(n+1)!}$ (47)
$\displaystyle=\frac{(2n-1)!}{2^{2n}(n+1)!(n-1)!}$ (48)

Passing to the gamma function $n!=\Gamma(n+1)$

$a_{2n}=\frac{\Gamma(2n)}{2^{2n}\Gamma(n)\Gamma(n+2)}$

Using the Legendre duplication identity
$\Gamma(2z)/[2^{2z}\Gamma(z)]=\Gamma(z+1/2)/(2\sqrt{\pi})$

$a_{2n}=\frac{1}{2\sqrt{\pi}}\frac{\Gamma(n+1/2)}{\Gamma(n+2)}$
$a_{m}=\frac{1}{2\sqrt{\pi}}\frac{\Gamma\left(\frac{m}{2}+\frac{1}{2}\right)}{\Gamma\left(\frac{m}{2}+2\right)}$

Regarding the coefficients of the series $b_{k}$, we can immediately deduce
that $b_{0}=0$, and may find from the recursion (43) and the factoring of
$a_{m}$ that, for $k\geq 2$

$b_{k}=\frac{1}{a_{k}}\frac{1}{k+2}$

We may now write the asymptotic expansion for the outer solution under
critically stressed conditions:

$\frac{\delta_{out}(\rho,t)}{Rf\Delta
p/\mu}=\frac{8}{\pi}\sum_{n=0}^{\infty}\frac{1}{\lambda^{2n+2}}\frac{1}{\sqrt{\pi}}\frac{\Gamma(n+3/2)}{2n+1}\frac{1}{\rho^{2n+1}}a_{2n}\left[\left(\text{acos}(\rho)-\rho\sqrt{1-\rho^{2}}\right)+\frac{(1-\rho^{2})^{3/2}}{\rho}\sum_{k=1}^{n}b_{2k}\rho^{2k}\right]$
$\displaystyle\begin{split}\frac{\delta_{out}(\rho,t)}{Rf\Delta
p/\mu}=\frac{4}{\pi^{2}}\sum_{n=0}^{\infty}&\frac{1}{\lambda^{2n+2}}\frac{1}{\rho^{2n+1}}\frac{\Gamma(n+3/2)}{2n+1}\frac{\Gamma(n+1/2)}{\Gamma(n+2)}\\\\[9.0pt]
&\cdot\left[\left(\text{acos}(\rho)-\rho\sqrt{1-\rho^{2}}\right)+2\sqrt{\pi}\frac{(1-\rho^{2})^{3/2}}{\rho}\sum_{k=1}^{n}\frac{\rho^{2k}}{2k+2}\frac{\Gamma(k+2)}{\Gamma(k+1/2)}\right]\end{split}$
(49)

##  Appendix C Asymptotic solution details for marginally pressurized
conditions

We pursue the asymptotic expansion of the solution in the marginally
pressurized limit, $\lambda\ll 1$. We return to the integral expression for
slip due to the constant-injection-rate source

$\delta(r,t)=\frac{8}{\pi}\frac{Rf\Delta
p}{\mu}\int_{\rho}^{1}\frac{\sqrt{\eta^{2}-\rho^{2}}}{\eta}\int_{0}^{1}\frac{xe^{-(x\eta\lambda)^{2}}}{\sqrt{1-x^{2}}}dx$

where we again identify the innermost inegral as $F(\zeta\lambda)$ and look to
expand it now for $\eta\lambda\ll 1$. In the above, we can substitute the
Taylor expansion

$e^{-(x\eta\lambda)^{2}}=\sum_{n=0}^{\infty}\frac{(-1)^{m}}{n!}(x\eta\lambda)^{n}$

and subsequently evaluate integrals of the form

$I_{2n+1}=\int_{0}^{1}\frac{x^{2n+1}}{\sqrt{1-x^{2}}}dx$

for $n=0,1,2...$ Letting $m=2n+1$ and $x=\sin\theta$ we transform the integral
to

$\displaystyle I_{m}$
$\displaystyle=\int_{0}^{1}\frac{x^{m}}{\sqrt{1-x^{2}}}dx=\int_{0}^{\pi/2}\sin^{m}\theta
d\theta$ $\displaystyle=\int_{0}^{\pi/2}\sin^{m-1}\theta\sin\theta
d\theta=\left(-\sin^{m-1}\theta\cos\theta\right|_{0}^{\pi/2}+(m-1)\int_{0}^{\pi/2}\sin^{m-2}\theta\cos^{2}\theta
d\theta$ $\displaystyle=(m-1)\left(\int_{0}^{\pi/2}\sin^{m-2}\theta
d\theta-\int_{0}^{\pi/2}\sin^{m}\theta d\theta\right)$
$\displaystyle=(m-1)(I_{m-2}-I_{m})$

and arrive to the recursion relation

$I_{m}=\frac{m-1}{m}I_{m-2}$

Given that $I_{1}=1$,

$I_{2n+1}=\frac{(2n)\cdot(2n-2)\cdot...\cdot 6\cdot 4\cdot
2}{(2n+1)\cdot(2n-1)\cdot...\cdot 7\cdot 5\cdot
3}=\frac{2^{n}n!}{(2n+1)!/(2^{n}n!)}$ (50)

which we may rewrite using the Gamma function

$I_{2n+1}=4^{n}\frac{\left[\Gamma(n+1)\right]^{2}}{\Gamma(2n+2)}$

and using the Legendre duplication formula we can further simplify

$I_{2n+1}=\frac{\sqrt{\pi}}{2}\frac{\Gamma(n+1)}{\Gamma(n+3/2)}$

and arrive to an expression for the expansion of $F$ for small values of its
argument

$F(\eta\lambda)\approx\frac{\sqrt{\pi}}{2}\sum_{n=0}^{\infty}(\eta\lambda)^{2n}\frac{(-1)^{n}}{\Gamma(n+3/2)}$

$\frac{\delta(\rho,t)}{Rf\Delta
p/\mu}=\frac{8}{\pi}\int_{\rho}^{1}\frac{\sqrt{\eta^{2}-\rho^{2}}}{\eta}F(\eta\lambda)d\eta=\frac{8}{\pi}\sum_{n=0}^{\infty}\frac{\sqrt{\pi}/2}{\Gamma(n+3/2)}\lambda^{2n}\int_{\rho}^{1}\sqrt{\eta^{2}-\rho^{2}}\eta^{2n-1}d\eta$

We denote the integrals within the sum as

$I_{2n-1}(\rho)=\int_{\rho}^{1}\sqrt{\eta^{2}-\rho^{2}}\eta^{2n-1}d\eta$

for $n=0,1,2...$

For $n=0$, we may evaluate $I_{-1}$ using the substitution
$\rho/\eta=\cos\theta$ and integration by parts

$\displaystyle
I_{-1}(\rho)=\int_{\rho}^{1}\frac{\sqrt{\eta^{2}-\rho^{2}}}{\eta}d\eta$
$\displaystyle=\rho\int_{0}^{\text{acos}(\rho)}\frac{\sin^{2}\theta}{\cos^{2}\theta}d\theta$
$\displaystyle=\rho\left[\left(\frac{\sin\theta}{\cos\theta}\right|_{0}^{\text{acos}(\rho)}-\int_{0}^{\text{acos}(\rho)}d\theta\right]$
$\displaystyle=\sqrt{1-\rho^{2}}-\rho\,\text{acos}(\rho)$

The remaining integrals in the sum ($n>0$) can be found by recursion. Letting
$m=2n-1$ and performing integration by parts

$\displaystyle I_{m}(\rho)$
$\displaystyle=\int_{\rho}^{1}\sqrt{\eta^{2}-\rho^{2}}\eta^{m}d\eta$
$\displaystyle=\frac{1}{3}\left.(\eta^{2}-\rho^{2})^{3/2}\eta^{m-1}\right|_{\rho}^{1}-\frac{m-1}{3}\int_{\rho}^{1}(\eta^{2}-\rho^{2})^{3/2}\eta^{m-2}d\eta$
$\displaystyle=\frac{1}{3}(1-\rho^{2})^{3/2}-\frac{m-1}{3}\int_{\rho}^{1}(\eta^{2}-\rho^{2})^{1/2}\eta^{m}(\eta^{2}-\rho^{2})\eta^{-2}d\eta$
$\displaystyle=\frac{1}{3}(1-\rho^{2})^{3/2}-\frac{m-1}{3}\left(I_{m}(\rho)-\rho^{2}I_{m-2}(\rho)\right)$

which leads to the recursion relation

$I_{m}(\rho)=\frac{(1-\rho^{2})^{3/2}}{m+2}+\frac{m-1}{m+2}\rho^{2}I_{m-2}(\rho)$

We may solve this recurrence relation to have an explicit expression for
$I_{m}(\rho)$. Given the form of the recursion, we posit a solution of the
form

$I_{m}(\rho)=\frac{(1-\rho^{2})^{3/2}}{m+2}\sum_{k=1,3,5...}^{m}a_{k}\rho^{k-1}$

and substitute this into the recursion relation to find that

$\sum_{k=1,3,...}^{m}a_{k}\rho^{k-1}=1+\frac{m-1}{m}\sum_{k=1,3,...}^{m}a_{k}\rho^{k+1}$

Expanding

$a_{1}+a_{3}\rho^{2}+a_{5}\rho^{4}+...=1+\frac{m-1}{m}\left(a_{1}\rho^{2}+a_{3}\rho^{4}+...\right)$

and comparing terms, we find that $a_{1}=1$ and subsequent terms follow the
recursion

$a_{m}=\frac{m-1}{m}a_{m-2}$

or

$a_{m}=\frac{(m-1)}{m}\cdot\frac{(m-3)}{(m-2)}\cdot\frac{(m-5)}{(m-4)}\cdot...\cdot\frac{4}{3}\cdot\frac{2}{1}\cdot
a_{1}$

Substituting $m=2q+1$ with $q=1,2,3...$

$a_{2q+1}=\frac{(2q)\cdot(2q-2)\cdot(2q-4)\cdot...\cdot 6\cdot 4\cdot
2}{(2q+1)\cdot(2q-1)\cdot(2q-3)\cdot...\cdot 5\cdot 3\cdot 1}\cdot a_{1}$

which we recognize as the series (50), which reduces to

$a_{2q+1}=\frac{\sqrt{\pi}}{2}\frac{\Gamma(q+1)}{\Gamma(q+3/2)}$

or

$a_{m}=\frac{\sqrt{\pi}}{2}\frac{\Gamma\left(\frac{m}{2}+\frac{1}{2}\right)}{\Gamma\left(\frac{m}{2}+1\right)}$

Therefore,

$I_{m}(\rho)=\frac{(1-\rho^{2})^{3/2}}{m+2}\sum_{k=1,3,...}^{m}\frac{\sqrt{\pi}}{2}\frac{\Gamma\left(\frac{k}{2}+\frac{1}{2}\right)}{\Gamma\left(\frac{k}{2}+1\right)}\rho^{k-1}$

or reverting to $m=2n-1$, with $n=1,2,3,...$

$I_{2n-1}(\rho)=\frac{(1-\rho^{2})^{3/2}}{2n+1}\sum_{k=1,2,...}^{n}\frac{\sqrt{\pi}}{2}\frac{\Gamma\left(k\right)}{\Gamma\left(k+1/2\right)}\rho^{2k-2}$

Such that the asymptotic expression for slip in the marginally pressurized
limit is

$\displaystyle\frac{\delta(\rho,t)}{R(t)f\Delta p/\mu}$
$\displaystyle=\frac{8}{\pi}\int_{\rho}^{1}\frac{\sqrt{\eta^{2}-\rho^{2}}}{\eta}d\eta+\frac{8}{\pi}\sum_{n=1}^{\infty}\lambda^{2n}\frac{\sqrt{\pi}/2}{\Gamma(n+3/2)}\int_{\rho}^{1}\sqrt{\eta^{2}-\rho^{2}}\eta^{2n-1}d\eta$
$\displaystyle=\frac{8}{\pi}\left(\sqrt{1-\rho^{2}}-\rho\text{
acos}(\rho)\right)+\frac{8}{\pi}\left(1-\rho^{2}\right)^{3/2}\sum_{n=1}^{\infty}\frac{\lambda^{2n}}{2n+1}\frac{\sqrt{\pi}/2}{\Gamma(n+3/2)}\sum_{k=1}^{n}\frac{\sqrt{\pi}}{2}\frac{\Gamma(k)}{\Gamma(k+1/2)}\rho^{2k-2}$
$\displaystyle=\frac{8}{\pi}\left(\sqrt{1-\rho^{2}}-\rho\text{
acos}(\rho)\right)+2\left(1-\rho^{2}\right)^{3/2}\sum_{n=1}^{\infty}\frac{\lambda^{2n}}{2n+1}\frac{1}{\Gamma(n+3/2)}\sum_{k=1}^{n}\frac{\Gamma(k)}{\Gamma(k+1/2)}\rho^{2k-2}$