# Pseudopotentials for Two-dimentional Ultracold Scattering in the Presence of
Synthetic Spin-orbit-coupling

Christiaan R. Hougaard, Brendan C. Mulkerin, Xia-Ji Liu, Hui Hu, Jia Wang
Centre for Quantum and Optical Science, Swinburne University of Technology,
Melbourne 3122, Australia.

###### Abstract

We derive a pseudopotential in two dimensions (2D) with the presence of a 2D
Rashba spin-orbit-coupling (SOC), following the same spirit of frame
transformation in [Phys. Rev. A 95, 020702(R) (2017)]. The frame
transformation correctly describes the non-trivial phase accumulation and
partial wave couplings due to the presence of SOC and gives rise to a
different pseudopotential than the free-space one, even when the length scale
of SOC is significantly larger than the two-body potential range. As an
application, we apply our pseudopotential with the Lippmann-Schwinger equation
to obtain an analytical scattering matrix. To demonstrate the validity, we
compare our results with a numerical scattering calculation of finite-range
potential and shows perfect agreement over a wide range of scattering energy
and SOC strength. Our results also indicate that the differences between our
pseudopotential and the original free-space pseudopotential are essential to
reproduce scattering observables correctly.

Modeling the fundamental two-body interactions is one of the critical steps in
investigating the complex quantum physics of a many-body system. In
particular, for systems with short-range interactions at low energies such as
ultracold quantum gases, the two-body interaction can be replaced by a zero-
range pseudopotential, giving the same wavefunction outside the original
potential. One only needs energy-dependent scattering lengths obtained via
partial expansion of two-body scattering to characterize the strength of such
pseudopotential. For example, in many cases, s-wave scattering dominates at
near-zero temperature, and the Fermi pseudopotential (Fermi, 1934; Huang and
Yang, 1957) gives a highly accurate description of the behavior of degenerated
quantum gases. For higher temperature beyond the Wigner-threshold regime,
generalizing the Fermi pseudopotential with an energy-dependent s-wave
scattering length can give quantitative descriptions (Bolda  _et al._ , 2002),
when the contribution from the higher-partial wave is negligible. However,
further generalization is necessary near resonances of higher-partial waves,
vanishes of s-wave scattering or when spin-orbit effects couples higher
partial waves. These regimes become particularly interesting in the field of
quantum gases, where interactions can be engineered at will via Feshbach
resonances (Chin  _et al._ , 2010) and synthetic spin-orbit-couplings (SOC)
(Lin  _et al._ , 2011). Even back at the initial development of
pseudopotential in the 1950s, Yang and Huang have already made an early
attempt to tackle the generalization for higher partial waves (Huang and Yang,
1957). However, they made an algebraic mistake in their original work leading
to an incorrect prefactor that is discovered and corrected much later (Stock
_et al._ , 2005; Derevianko, 2005; Idziaszek and Calarco, 2006). With the
corrected prefactor, the mean-field energy shift of interacting fermions in a
trap accurately matches experimental measurements (Roth and Feldmeier, 2001).

These pioneer works mentioned so far all focused on pseudopotentials in three-
dimensional (3D) space. Extensions to lower dimensions have been only of
academic interest until an exciting development in the area of ultracold
quantum gases recently: the creation of low-dimensional quantum gases.
Confining quantum gases into lower dimensions can strongly enhance the quantum
correlation of the system, leading to qualitative changes (Bloch  _et al._ ,
2008). Experimental realization of quasi-one-dimensional (quasi-1D) bosons
allowed the verification of the fermionization of 1D Bose gases in the
Tonks–Girardeau regime (Paredes  _et al._ , 2004; Kinoshita  _et al._ , 2004).
More recently, tightly confining the motion of atoms in one direction creates
quasi-two-dimensional (quasi-2D) quantum gases that have several applications
in investigating intrigued physics, such as the observation of the BKT phase
(Hadzibabic  _et al._ , 2004; Cladé  _et al._ , 2009), measurement of the
equation of state (Rath  _et al._ , 2010) and super-Efimov physics (Nishida
_et al._ , 2013; Gao  _et al._ , 2015). The experimental realization of
quasi-2D quantum gas motivates an elegant derivation of 2D pseudopotential for
arbitrary partial waves (Kanjilal and Blume, 2006).

Another invaluable development in ultracold quantum gases in recent years is
the realization of synthetic gauge fields that can be used to simulate
electromagnetic interactions in systems of neutral particles. Artificial gauge
fields can also couple a particle’s canonical momentum with its (pseudo)spin
degrees of freedom (Lin  _et al._ , 2011; Zhang  _et al._ , 2014; Zhai, 2015),
providing an essential ingredient, namely SOC, for the study of topological
insulators (Dalibard  _et al._ , 2011; Goldman  _et al._ , 2014). The
interplay between SOC and short-range interactions might lead to new quantum
behaviors and phases, and soon attracts a lot of interest. In cold-atom
systems, several experimental techniques have been developed to realize SOC
such as lattice shaking (Struck  _et al._ , 2012) and Raman coupling (Lin  _et
al._ , 2011). While the Raman laser scheme has already achieved one-
dimensional SOC (an equal mixture of Rashba and Dresselhaus SOC) (Lin  _et
al._ , 2011; Cheuk  _et al._ , 2012; Wang  _et al._ , 2012; Zhang  _et al._ ,
2012a; Qu  _et al._ , 2013; Khamehchi  _et al._ , 2017), SOC with higher
symmetry such as 2D and 3D isotropic ones are more closely related to the
cases in condensed-matter physics. On the other hand, 2D Rashba (isotropic)
SOC, which is our main focus in this work, is more experimentally accessible
than 3D with less problem from heating (Huang  _et al._ , 2016; Meng  _et al._
, 2016).

In previous theoretical studies, people usually directly apply the original
free-space pseudopotentials (obtained from two-body scattering without SOC) in
the presence of SOC (Zhang  _et al._ , 2014; Zhai, 2015). The justification
base on the argument that the characteristic wave-length of synthetic SOC, in
reality, is much larger than the inverse of the range of short-range
interaction. Therefore the effects of SOC were assumed to have no impact on
behaviors of wavefunctions at short inter-particle distances, and hence the
original pseudopotential remains valid from a perturbation point of view.
Nevertheless, in a foresighted study, Cui points out that the presence of SOC
at short distances intrinsically mixes different partial waves via the
couplings of spin, which might lead to non-trivial influence on the short-
range wavefunction (Cui, 2012). Via some numerical investigations, Cui
concludes that free-space pseudopotentials, especially for higher partial
waves such as $p$-wave, are not satisfactory. Ever since, several studies have
carefully calculated two-body scattering with the presence of 3D (Zhang  _et
al._ , 2012b; Yu, 2012; Zhang  _et al._ , 2013a; Duan  _et al._ , 2013; Wang
and Greene, 2015; Guan and Blume, 2016; Wang  _et al._ , 2018) or 2D (Zhang
_et al._ , 2012c, 2013b) SOC, paving the way for designing a pseudopotential
model. One particularly enlightening study carried out by Guan and Blume
reveals that a frame transformation approach (that we will detail later) can
correctly calculate the scattering phase accumulated at short distances
modified by SOC (Guan and Blume, 2017). However, a proper pseudopotential that
includes the nonperturbative effects of SOC at short-range and correctly
reproduces scattering observables is still missing. In this Rapid
Communication, we derive an analytical form of the pseudopotential in 2D with
the presence of 2D Rashba SOC, following the same spirit of frame
transformation in Ref. (Guan and Blume, 2017). To verify the validity, we
apply the Lippmann-Schwinger equation to obtain the analytical scattering
matrix and compare it with a numerical scattering calculation with finite-
range potential.

We first give a brief review of 2D pseudopotential in the free-space without
the presence of SOC. We consider two identical particles ($n=1,2$) of mass $m$
confined in a 2D $x$-$y$ plane with position vectors $\mathbf{r}_{n}$.
Seperating out the center-of-mass (COM) motion, the Hamiltonian of the
relative motion is given by $H^{{\rm fs}}=\mathbf{p}^{2}/2\mu_{2b}+U(\rho)$,
where $\mu_{2b}=m/2$ is the two-body reduced mass,
$\mathbf{r}=\\{\rho,\phi\\}$ is the relative position in polar coordinates,
and $\mathbf{p}=-i\hbar\\{\partial_{\rho},\rho^{-1}\partial_{\phi}\\}$ is the
relative momentum in 2D. We also assume the potential $U(\rho)$ is isotropic
and short-range, i.e., vanishes beyond a small radius $\rho_{0}$. The
isotropic symmetry allows the wavefunction to be expanded as $\Psi^{{\rm
fs}}(\mathbf{r})=\sum_{m_{\ell}}R_{m_{\ell}}^{{\rm
fs}}(\rho)\Phi_{m_{\ell}}(\phi)$, where
$\Phi_{m_{\ell}}(\phi)=e^{im_{\ell}\phi}/\sqrt{2\pi}$, $R_{m_{\ell}}^{{\rm
fs}}(\rho)$ satisfies the radial Schrödinger equation

$\left[-\frac{\hbar^{2}}{2\mu_{2b}}\left(\frac{\partial^{2}}{\partial\rho^{2}}+\frac{1}{\rho}\frac{\partial}{\partial\rho}-\frac{m_{\ell}^{2}}{\rho^{2}}\right)+U(\rho)-E\right]R_{m_{\ell}}^{{\rm
fs}}(\rho)=0$ (1)

and adopts an asymptotic form $R_{m_{\ell}}^{{\rm fs}}(\rho)\propto
J_{m_{\ell}}(k\rho)-\tan[\delta_{m_{\ell}}(k)]N_{m_{\ell}}(k\rho)$ for
$\rho>\rho_{0}$. Here $E=\hbar^{2}k^{2}/2\mu_{2b}$ and
$J_{m_{\ell}}(k_{\tau}r)$ and $Y_{m_{\ell}}(k_{\tau}r)$ are the Bessel
functions of the first and second kind respectively. $\delta_{m_{\ell}}(k)$
are the energy-dependent phase shifts, satisying threshold law
$\tan[\delta_{0}(k)]\propto 1/\log k$ and $\tan[\delta_{m_{\ell}}(k)]\propto
1/k^{2|m_{\ell}|}$ for $|m_{\ell}|\geq 1$. Reference (Kanjilal and Blume,
2006) shows that replacing $U(\rho)$ by a pseudopotential $V_{m_{\ell}}^{{\rm
fs}}(\rho)$ can give the same asymptotic wavefunction, and hence reproduce the
low-energy observables of the original finite-range potential. The explicit
form of $V_{m_{\ell}}^{{\rm fs}}(\rho)$ in free space is given by

$V_{m_{\ell}}^{{\rm
fs}}(\rho,k)=-\frac{\hbar^{2}}{\mu_{2b}}\frac{\tan\left[\delta_{m_{\ell}}(k)\right]}{c_{m_{\ell}}k^{2m_{\ell}}\rho^{m_{\ell}}}\left[\frac{\delta(\rho-s)}{2\pi\rho}\hat{O}_{m_{\ell}}(\rho,k)\right]_{s\rightarrow
0^{+}},$ (2)

where $c_{m_{\ell}}=(2m_{\ell})!/[\Gamma(m_{\ell}+1)]^{2}2^{2m_{\ell}}$ with
$\Gamma(\cdot)$ being the gamma function. The form of delta shell
$\delta(\rho-s)$ of radius $s$ approaches to a contact potential
$\delta(\rho)$ in the limit $s\rightarrow 0$, and allows us to deal with the
divergence of the regularized operator rigorously. The regularized operator
reads as

$\hat{O}_{m_{\ell}}=\left\\{\begin{array}[]{ll}{\frac{2}{1-\tan\left[\delta_{0}(k)\right]f_{0}(k,\rho)}\frac{\partial}{\partial\rho}\rho}&{;m_{\ell}=0}\\\
{\frac{2}{1-\tan\left[\delta_{m_{\ell}}(k)\right]f_{m_{\ell}}(k,\rho)}\frac{\partial^{2m_{\ell}}}{\partial\rho^{2m_{\ell}}}\rho^{m_{\ell}}}&{;m_{\ell}>0}\end{array}\right.$
(3)

where

$f_{m_{\ell}}(k,\rho)=\left\\{\begin{array}[]{ll}{\frac{2}{\pi}\left[1+\gamma+\log\left(\frac{1}{2}k\rho\right)\right]}&{;m_{\ell}=0}\\\
{\frac{2}{\pi}\left[\sum_{n=0}^{2m_{\ell}-1}\frac{1}{2m_{\ell}-n}-\bar{\psi}+\log\left(\frac{1}{2}k\rho\right)\right]}&{;m_{\ell}>0}\end{array}\right..$
(4)

Here $2\bar{\psi}(m_{\ell})=\psi(1)+\psi(m_{\ell}+1)$, where $\psi$ denotes
the digamma function. For $m_{\ell}<0$, the pseudopotential in Eq. (2) takes
the same form but with $m_{\ell}$ replaced by $|m_{\ell}|$. In contrast to the
3D pseudopotential, the $\tan\left[\delta_{m_{\ell}}(k)\right]$ dependence in
the denominator of the regularized operators originates from the fact that
$\frac{\partial}{\partial\rho}\rho N_{0}(k\rho)$ and
$\frac{\partial^{2m_{\ell}}}{\partial\rho^{2m_{\ell}}}\rho^{m_{\ell}}N_{m_{\ell}}(k\rho)$
does not vanishes at $\rho\rightarrow 0$.

Now we consier the effects of SOC, where each particle feels a 2D Rashba SOC
described by $H_{{\rm SO}}^{(n)}=k_{{\rm
SO}}\mathbf{p}_{n}\cdot\mathbf{s}_{n}/m$, with $\mathbf{p}_{n}$ and
$\mathbf{s}_{n}$ being the 2D momentum and spin operator of particle $n$
respectively. Following the spirit of Refs. (Duan  _et al._ , 2013; Wang and
Greene, 2015; Guan and Blume, 2016, 2017; Wang  _et al._ , 2018), we focus on
the scattering in the COM frame, where the relative Hamiltonian can be written
as $H_{{\rm rel}}=H^{{\rm fs}}+V^{{\rm SO}}$ with $V^{{\rm SO}}=k_{{\rm
SO}}\mathbf{\Sigma}\cdot\mathbf{p}/2\mu_{2b}$ describing the SOC effect. and
$\mathbf{\Sigma}=\mathbf{s_{1}-\mathbf{s_{2}}}$ is the relative spin operator.
$k_{{\rm SO}}$ defines the strength of SOC coupling, and gives an energy scale
$E_{{\rm SO}}=\hbar^{2}k_{{\rm SO}}^{2}/2m$.

A formal way to solve the corresponding relative Schrödinger equation is to
formulate it as a multichannel problem by expanding the $\tau$’th independent
solution as

$\Psi_{\tau}^{{\rm SO}}(\mathbf{r})=\sum_{\nu}R_{\nu\tau}^{{\rm
SO}}(\rho)A_{\nu}(\Omega),$ (5)

where the channel functions $A_{\nu}(\Omega)\equiv\langle\Omega|\nu\rangle$
are functions of $\Omega$ that includes all degrees of freedom except for
$\rho$. Due to the azimutual symmetry, total angular momentum (along $z$-axis)
$m_{j}$ is a good quantum number that equals to $m_{\ell}+m_{S}$. Here
$m_{1}$, $m_{2}$ and $m_{S}=m_{1}+m_{2}$ are quantum number of the projection
of the operator $\mathbf{s}_{1}$, $\mathbf{s}_{2}$ and
$\mathbf{S}=\mathbf{s}_{1}+\mathbf{\mathbf{s}}_{2}$ to the quantization
$z$-axis respectively. Defining the total spin basis
$\left|\chi\right\rangle\equiv\left|(s_{1,}s_{2}),S,m_{S}\right\rangle$ as
usual, we choose the channel functions being
$A_{\nu}(\Omega)=i^{m_{\ell}}\Phi_{m_{\ell}}(\phi)\left|\chi\right\rangle$
with $m_{\ell}+m_{S}=m_{j}$ and $S+m_{\ell}+s_{1}+s_{2}$ being even/odd for
bosons/fermions respectively. The subindex $\nu$ collectively represents
quantum numbers $\\{m_{\ell},S,m_{S};m_{j}\\}$. Here we omit quantum numbers
$s_{1}$ and $s_{2}$ in the channel index notations since they are the same for
all channels. At $\rho>\rho_{{\rm 0}},$ wavefuctions can be expressed as a
linear combination of non-interacting (but with SOC) regular and irregular
solution $\underline{R}^{{\rm
SO}}=\underline{F}-\underline{G}\underline{\mathcal{K}}$, where
$\underline{R}^{{\rm SO}}$ is the matrix form of raidal solutions
$R_{\nu\tau}^{{\rm SO}}(\rho)$. (Through out this paper, underline implies
matrix form.) The matrix elements of regular solution $\underline{F}$ can be
written as
$F_{\nu\tau}(\rho)=N_{\tau}C_{\nu\tau}\sqrt{k_{\tau}}J_{m_{\ell}}(k_{\tau}\rho)$,
where $C_{\nu\tau}$, $k_{\tau}$ and $N_{\tau}$ can be obtained by
diagonalizing the non-interacting Hamiltonian using the same procedure as Ref.
(Wang and Greene, 2015). The corresponding irregular solution can be obtained
as
$G_{\nu\tau}(\rho)=N_{\tau}C_{\nu\tau}\sqrt{k_{\tau}}Y_{m_{\ell}}(k_{\tau}\rho)$.
The scattering $K$ matrix $\underline{\mathcal{K}}$ determines scattering
observables and is related to the more familiar $S$ matrix by
$\underline{\mathcal{S}}=(\underline{I}+i\underline{\mathcal{K}})(\underline{I}-i\underline{\mathcal{K}})^{-1}$.
Our goal is to replace the potential $U(\rho)$ by a potential
$\underline{V}^{(m_{j})}(\rho)$ that acts only at $\rho=0$, and gives the same
asymptotic wavefunction and hence the same $K$ matrix. Here the underline
indicates $\underline{V}^{(m_{j})}(\rho)$ is a matrix and not nessesary
diagonal due to the presene of SOC.

To derive this pseudopotential, we follow the spirit of Ref. (Guan and Blume,
2017), and apply a frame transformation approach. Defining a unitary
transformation $\mathcal{U}_{1}=\exp(-ik_{{\rm
SO}}\mathbf{\Sigma}\cdot\mathbf{r}/2\hbar)$, the “rotated” Hamiltonian
$H^{{\rm temp}}\equiv\mathcal{U}_{1}^{-1}H_{{\rm rel}}\mathcal{U}_{1}=H^{{\rm
fs}}+\mathcal{\epsilon^{{\rm temp}}}+\mathcal{O}(\rho)$ is introduced as an
intermediate step. Here we neglect terms of order $\rho$ and higher denoted by
$\mathcal{O}(\rho)$, since the pseudopotential will only act at $\rho=0$. The
constant term $\mathcal{\epsilon^{{\rm temp}}}$ is given by
$\mathcal{\epsilon^{{\rm temp}}}=-E_{{\rm
SO}}(\Sigma_{x}^{2}+\Sigma_{y}^{2}+S_{z}L_{z})/2\hbar^{2}$, where $\Sigma_{x}$
($\Sigma_{y}$) are the $x$($y$)-component of $\mathbf{\Sigma}$, $S_{z}$ is the
$z$-component of total spin operator $\mathbf{S}$, and
$L_{z}=-i\hbar\partial_{\phi}$ is the 2D angular momentum operator. For two
spin-1/2 particles, this operator expanded by channel functions
$A_{\nu}(\Omega)$ gives a diagonal matrix $\mathcal{\underline{\epsilon}^{{\rm
temp}}}$ . In contrary, for higher spins, $\mathcal{\underline{\epsilon}^{{\rm
temp}}}$ is in general not diagonal, where the only non-zero matrix elements
are the ones couples channels with the same $m_{S}$ and hence the same
$m_{\ell}$. Therefore, one can introduc another $\rho$-indepdent unitary
transformation $\mathcal{U}_{2}$ that is block-diagonal in $m_{\ell}$
subspaces and satisfies
$\underline{\epsilon}=\mathcal{U}_{2}^{-1}\underline{\epsilon}^{{\rm
temp}}\mathcal{U}_{2}$ is diagonal. Automatically,
$\mathcal{U}_{2}^{-1}H^{{\rm fs}}\mathcal{U}_{2}=H^{{\rm fs}}$ is also
diagonal. Therefore, we find a unitary transformation
$\mathcal{U}=\mathcal{U}_{2}\mathcal{U}_{1}$ that leads to a set of uncoupled
radial Schrödinger equation that is at least valid near the origin $\rho=0$,
which is given by,

$\left[-\frac{\hbar^{2}}{2\mu_{2b}}\left(\frac{\partial^{2}}{\partial\rho^{2}}+\frac{1}{\rho}\frac{\partial}{\partial\rho}-\frac{m_{\ell}^{2}}{\rho^{2}}\right)+U(\rho)-E_{\nu}\right]\tilde{R}_{\nu\nu}^{{\rm
SO}}(\rho)=0,$ (6)

where $E_{\nu}\equiv E+\epsilon_{\nu}$ with $\epsilon_{\nu}$ being the
diagonal matrix elements of SOC-induced energy shift $\underline{\epsilon}$.
Comparing with the free-space Schrödinger equation in Eq. (1), a
pseudopotential $\underline{\tilde{V}}^{(m_{j})}(\rho,k)$ with diagonal matrix
elements $\tilde{V}_{\nu\nu}^{(m_{j})}(\rho,k)=V_{m_{\ell}}^{{\rm
fs}}(\rho,k_{\nu})$ where $k_{\nu}^{2}/2\mu_{2b}=E_{\nu}$ can reproduce the
wavefunction $\tilde{R}_{\nu\nu}^{{\rm SO}}(\rho)$ in the rotated frame. The
pseudopotential in the original frame can therefore be obtained by an inverse
rotation:

$\underline{V}^{(m_{j})}(\rho,k)=\mathcal{U}\tilde{\underline{V}}^{(m_{j})}(\rho,k)\mathcal{U}^{-1}.$
(7)

Once we obtained $\underline{V}^{(m_{j})}(\rho,k)$ for all $m_{j}$ (or up to a
cut-off $m_{j}$ in practice), the total pseudopotential can be expressed as

$\begin{gathered}V(\rho,k)\Psi^{{\rm
SO}}(\mathbf{r})=\sum_{m_{j}}\sum_{\nu\nu^{\prime}}A_{\nu}(\Omega)\int\rho^{\prime}d\rho^{\prime}\underline{V}_{\nu\nu^{\prime}}^{(m_{j})}(\rho^{\prime})\\\
\times\int d\Omega^{\prime}A_{\nu^{\prime}}^{*}(\Omega^{\prime})\Psi^{{\rm
SO}}(\mathbf{r}^{\prime})\end{gathered},$ (8)

which might be helpful for applications where $m_{j}$ is not a good quantum
number.

For illustration, we consider two spin-$1/2$ fermions in the $m_{j}=0$
subspace, and omit the notation of $m_{j}$ hereafter unless specify otherwise.
The basis are denoted by
$\nu\equiv\\{m_{\ell},S,m_{S}\\}=\\{-1,1,1\\},\\{0,0,0\\},\\{1,1,-1\\}$. In
this order of the basis, the rotational matrix can be written out explicitely

$\underline{\mathcal{U}}=\left[\begin{array}[]{ccc}{\cos^{2}\left(\frac{\lambda\rho}{2}\right)}&{-\frac{\sin(\lambda\rho)}{\sqrt{2}}}&{-\sin^{2}\left(\frac{\lambda\rho}{2}\right)}\\\
{\frac{\sin(\lambda\rho)}{\sqrt{2}}}&{\cos(\lambda\rho)}&{\frac{\sin(\lambda\rho}{\sqrt{2}}}\\\
{-\sin^{2}\left(\frac{\lambda\rho}{2}\right)}&{-\frac{\sin(\lambda\rho)}{\sqrt{2}}}&{\cos^{2}\left(\frac{\lambda\rho}{2}\right)}\end{array}\right],$
(9)

where $\lambda=k_{{\rm SO}}/2$ is introduced for convenience. After rotation,
the SOC-induced energy shift is given by $\underline{\epsilon}={\rm
diag}[0,\hbar^{2}\lambda^{2}/\mu_{2b},0]$, where ${\rm diag[\cdot]}$
represents a diagonal matrix. The energy shift determines the pseudopotential
in the rotated frame as $\tilde{\underline{V}}(\rho,k)={\rm diag}[V_{1}^{{\rm
fs}}(\rho,k_{p}),V_{0}^{{\rm fs}}(\rho,k_{s}),V_{1}^{{\rm fs}}(\rho,k_{p})]$,
where $k_{s}=\sqrt{k^{2}+2\lambda^{2}}$ and $k_{p}=k$. We then apply Eq. (7)
to obtain the pseudopotential in the original frame. We can write the the
pseudopotential as a summation of $s$\- and $p$\- wave contribution:

$\underline{V}=-\frac{\hbar^{2}}{\mu_{2b}}\left\\{\frac{\delta(\rho-s)}{2\pi\rho}\left[\tan\delta_{s}(k_{s})\underline{O}_{s}+\frac{\tan\delta_{p}(k_{p})}{k_{p}^{2}}\underline{O}_{p}\right]\right\\}_{s\rightarrow
0},$ (10)

where

$\underline{O}_{s}=\left[\begin{array}[]{ccc}0&0&0\\\
-\frac{\lambda}{\sqrt{2}}\hat{O}_{0}\rho&\frac{1}{\rho}\hat{O}_{0}&-\frac{\lambda}{\sqrt{2}}\hat{O}_{0}\rho\\\
0&0&0\end{array}\right],$ (11)

and

$\underline{O}_{p}=\left[\begin{array}[]{ccc}\frac{1}{\rho}\hat{O}_{1}\left(1-\frac{\lambda^{2}\rho^{2}}{4}\right)&\frac{\lambda}{\rho\sqrt{2}}\hat{O}_{1}\rho&-\frac{1}{\rho}\hat{O}_{1}\frac{\lambda^{2}\rho^{2}}{4}\\\
\frac{\lambda}{\sqrt{2}}\hat{O}_{1}\left(1-\frac{\lambda^{2}\rho^{2}}{2}\right)&\lambda^{2}\hat{O}_{1}\rho&\frac{\lambda}{\sqrt{2}}\hat{O}_{1}\left(1-\frac{\lambda^{2}\rho^{2}}{2}\right)\\\
-\frac{1}{\rho}\hat{O}_{1}\frac{\lambda^{2}\rho^{2}}{4}&\frac{\lambda}{\rho\sqrt{2}}\hat{O}_{1}\rho&\frac{1}{\rho}\hat{O}_{1}\left(1-\frac{\lambda^{2}\rho^{2}}{4}\right)\end{array}\right].$
(12)

Figure 1: Results near $s$-wave resonances with $\lambda r_{{\rm vdW}}=0.01$.
The free-space scattering phase-shifts are obtained from a Lennard-Jones model
potential with parameter $r_{0}=0.58r_{{\rm vdW}}$ that leads to
$a_{s}(0)\approx 21.777r_{{\rm vdW}}$. (a) and (b) partial cross-sections. (c)
and (d) $K$ matrix elements. The curves represents analytical results, and the
symbols are obtained by numerical scattering calculations. Figure 2: Results
near $p$-wave resonances with $\lambda r_{{\rm vdW}}=0.01$. The free-space
scattering phase-shifts are obtained from a Lennard-Jones model potential with
parameter $r_{0}=0.552981r_{{\rm vdW}}$ that leads to
$A_{p}\approx-8.577\times 10^{5}r_{{\rm vdW}}^{2}$. (a) and (b) partial cross-
sections. (c) and (d) $K$ matrix elements. The curves represents analytical
results, and the symbols are obtained by numerical scattering calculations.

Terms of higher order of $\rho$ can be ignored with the consideration that the
pseudopotential only contribute to $K$ matrix with terms proportional to
$F_{\nu^{\prime}\tau^{\prime}}^{*}(s)V_{\nu\nu^{\prime}}(s,k)F_{\nu\tau}(s)_{s\rightarrow
0}$ and
$F_{\nu^{\prime}\tau^{\prime}}^{*}(s)V_{\nu\nu^{\prime}}(s,k)G_{\nu\tau}(s)_{s\rightarrow
0}$ [see Eqs. (15) below]. Comparing with the the free-space pseudopotential
${\rm diag}[V_{1}^{{\rm fs}}(k,\rho),V_{0}^{{\rm fs}}(k,\rho),V_{1}^{{\rm
fs}}(k,\rho)]$, there are two important differences. One is the SOC-induced
energy shift leads to a different $s$-wave phase shift $\delta_{s}(k_{s})$.
The other is the non-diagonal terms rised from the rotational transformation,
which describes the intrinsical partial waves mixing at short distances
induced by SOC. As we will see, both of these differences play significant
roles in producing scattering observables correctly.

To verify the validity of the pseudopotential, we apply the Lippmann-Schwinger
equation to calulate the $K$ matrix. The Lippmann-Schwinger equation is the
integral form of the Schrödinger equation
$\Psi_{\tau}(\mathbf{r})=\Psi_{0}(\mathbf{r})+\int
G(\mathbf{r},\mathbf{r^{\prime}})V(\mathbf{r^{\prime}})\Psi_{\tau}(\mathbf{r}^{\prime})d\mathbf{r}^{\prime}$,
or equivalently in the matrix form

$\underline{R}^{{\rm
SO}}(\rho)=\underline{F}(\rho)+\int\underline{\mathcal{G}}(\rho,\rho^{\prime})\underline{V}(\rho^{\prime})\underline{R}(\rho^{\prime})\rho^{\prime}d\rho^{\prime}.$
(13)

Here $\underline{\mathcal{G}}(\rho,\rho^{\prime})$ is the matrix
representation of the Green’s function
$\mathcal{G}(\mathbf{r},\mathbf{r^{\prime}})=\sum_{\nu\nu^{\prime}}A_{\nu}(\Omega)\underline{\mathcal{G}}_{\nu\nu^{\prime}}(\rho,\rho^{\prime})A_{\nu^{\prime}}^{*}(\Omega)$,
which is given by

$\underline{\mathcal{G}}\left(\rho,\rho^{\prime}\right)=\pi\left\\{\begin{array}[]{ll}{\underline{F}(\rho)G^{\dagger}\left(\rho^{\prime}\right),}&{\rho<\rho^{\prime}}\\\
{\underline{G}(\rho)\underline{F}^{\dagger}\left(\rho^{\prime}\right),}&{\rho>\rho^{\prime}}\end{array}.\right.$
(14)

The Green’s function approach becomes very helpful when the potential can be
replaced by a pseudopotential
$\underline{V}(\rho^{\prime})\propto\delta(\rho^{\prime})$, and the $K$ matrix
can then be obtained by
$\underline{\mathcal{K}}=(I+\underline{\mathcal{B}})^{-1}\mathcal{\underline{A}}$,
where

$\begin{gathered}\underline{\mathcal{A}}=-\pi\int\rho^{\prime}d\rho^{\prime}\underline{F}(\rho^{\prime})\underline{V}(\rho^{\prime})\underline{F}(\rho^{\prime}),\\\
\underline{\mathcal{B}}=-\pi\int\rho^{\prime}d\rho^{\prime}\underline{F}(\rho^{\prime})\underline{V}(\rho^{\prime})\underline{G}(\rho^{\prime}).\end{gathered}$
(15)

Noticing it is a special property of 2D that $\mathcal{B}$ does not vanish. As
a direct consequence, the $K$ matrix in general cannot be written as a
summation of $s$\- and $p$-wave contribution in contrast to the 3D case as
shown in Eq. (11) of Ref. (Guan and Blume, 2017).

For illustration, we focus on the case of two spin-$1/2$ fermions in the
$m_{j}=0$ subspace. The regular solution $\underline{F}$ and irregular
solution $\underline{G}$ can be determined by the coefficient $C_{\nu\tau}$ in
a matrix form

$\underline{C}=\left[\begin{array}[]{ccc}{-1/2}&{-1/2}&{1/\sqrt{2}}\\\
{-1/\sqrt{2}}&{1/\sqrt{2}}&{0}\\\ {1/2}&{1/2}&{1/\sqrt{2}}\end{array}\right],$
(16)

where the the column index $\tau$ corresponding to cannonical momentum
$\\{k_{\tau}\\}\equiv\\{k_{1},k_{2},k_{3}\\}=\\{k_{b}+\lambda,k_{b}-\lambda,k\\}$
and normalization
$\\{\hbar^{2}N_{\tau}^{2}/\mu_{2b}\\}=\left\\{1/k_{b},1/k_{b},1/k\right\\}$
where $k_{b}=\sqrt{\lambda^{2}+k^{2}}$. One can identify $\tau=\\{1,2,3\\}$
corresponds to three different configurations $|-,-\rangle$, $|+,+\rangle$ and
$|-,+\rangle$, where $-$ ($+$) indicates the helicity, i.e. whether the spin
is anti-parallel/parallel to the direction of current (Duan  _et al._ , 2013).

Inserting $\underline{F}$ and $\underline{G}$ into Eq. (15) gives
$\underline{\mathcal{A}}$ and $\underline{\mathcal{B}}$ that determines
$\underline{\mathcal{K}}$. We find that the $K$ matrix is block-diagonal and
can be expressed as

$\mathcal{\underline{K}}=\left[\begin{array}[]{cc}\underline{\mathcal{K}}^{(+)}&0\\\
0&\underline{\mathcal{K}}^{(-)}\end{array}\right],$ (17)

where $\underline{\mathcal{K}}^{(-)}=\tan[\delta_{p}(k_{p})]$, and
$\underline{\mathcal{K}}^{(+)}$ is a $2\times 2$ matrix. The block-diagonal
structure can be understood by studying the $\mathcal{\mathcal{PT}}$-symmetry
of $\sum_{\nu}F_{\nu\tau}(\rho)A_{\nu}(\Omega)$ in the $m_{j}=0$ subspace,
where $\mathcal{P}$ is defined as $\phi\rightarrow\phi+\pi$ and $\mathcal{T}$
is defined as
$\left|s_{n},m_{n}\right\rangle\rightarrow\left|s_{n},-m_{n}\right\rangle$ for
both $n=1,2$. Defining $\mathcal{P}\mathcal{T}[\sum
F_{\nu\tau}(\rho)A_{\nu}(\Omega)]=\Pi_{\tau}[\sum
F_{\nu\tau}(\rho)A_{\nu}(\Omega)]$, one finds that $\Pi_{\tau}=+1$ /$-1$ for
$\tau=\\{1,2\\}$/$\\{3\\}$, leading to the block-diagonal structure. [For
$m_{j}\neq 0$ subspaces, ${\rm sign}(m_{j})$ and $\mathcal{P}\mathcal{T}$
cannot simutaneously be good quantum numbers]. The expressions of matrix
elements of $\underline{\mathcal{K}}^{(+)}$ are analytical but quite
cumbersome, and hence we only give the full expression in the supplemental
materials and illustrate them in Fig. 1 and Fig. 2 as two numerical examples
near $s$\- and $p$-wave resonances respectively. In these two examples, the
energy-dependent phase-shifts $\delta_{s}(k)$ and $\delta_{p}(k)$ are obtained
from a free-space scattering calculation with Lenard-Jones potential
$U(\rho)=-\frac{C_{6}}{\rho^{6}}\left(1-\frac{r_{0}^{6}}{\rho^{6}}\right)$,
where $C_{6}$ defines a length scale $r_{{\rm
vdW}}=(2\mu_{2b}C_{6}/\hbar^{2})^{1/4}/2$ and $r_{0}$ controls short-range
physics and is used to tune zero-energy scattering phase shifts. The
analytical results shows a good agreement with a full numerical scatteirng
calculation with the representation of SOC, using a similar procedure as Ref.
(Wang  _et al._ , 2018). The technical details are shown in the Supplemental
Material.

Near $s$-wave resonance, $p$-wave scattering is negligible and the scattering
matrix can be simplified as

$\underline{\mathcal{K}}^{(s)}=T_{0}\left[\begin{array}[]{cc}k_{1}&-k\\\
-k&k_{2}\end{array}\right],$ (18)

where $T_{0}=\tan[\delta_{s}(k_{s})]/2k_{b}\alpha_{s}$ and
$\alpha_{s}=1+\frac{2}{\pi}\tan[\delta_{s}(k_{s})]\left(\log\frac{k}{k_{s}}+\frac{\lambda}{k_{b}}\tanh^{-1}\frac{\lambda}{k_{b}}\right)$.
The $S$ matrix can be obtained via
$\underline{\mathcal{S}}^{(s)}=(I+i\mathcal{\underline{\mathcal{K}}}^{(s)})(I-i\mathcal{\underline{K}}^{(s)})^{-1}$
and determines the scattering cross-section
$\sigma_{\tau\tau^{\prime}}^{(s)}=2|\mathcal{\underline{S}}_{\tau\tau^{\prime}}^{(s)}-\delta_{\tau\tau^{\prime}}|^{2}/k_{\tau}$
that are read as
$\sigma_{11}^{(s)}=\sigma_{21}^{(s)}=k_{1}8T_{0}^{2}/(1+4T_{0}^{2}k_{b}^{2})$
and
$\sigma_{12}^{(s)}=\sigma_{22}^{(s)}=k_{2}8T_{0}^{2}/(1+4T_{0}^{2}k_{b}^{2})$.
As a result, $\sigma_{21}^{(s)}/\sigma_{12}^{(s)}=k_{1}/k_{2}>1$ indicates
that particles are preferentially scattered into the lower energy helicity
“$-$” state. The validity of Eq. (18) near $s$-wave resonances is varified in
Fig. 1. In the zero-energy limit, $\lambda\sigma_{11}^{(s)}\rightarrow
4/\left\\{1+[\gamma+\frac{2}{\pi}\log\left(\lambda
a_{s}(\sqrt{2}\lambda)\right)]^{2}\right\\}$, where $a_{s}(k)$ is the
generalized energy-dependent $s$-wave scattering length defined by
$\cot[\delta_{s}(k)]=\frac{2}{\pi}\log[k_{s}a_{s}(k)]+\gamma$. The rescaled
cross-section therefore reach maximum when $a_{s}(k_{s})$ equals to $a_{{\rm
res}}\equiv e^{-\pi\gamma/2}/\lambda$. In comparison, if we replace $U(\rho)$
directly by ${\rm diag}[0,V_{0}^{{\rm fs}}(k,\rho),0]$, the free-space
pseudopotential with $s$-wave only, and apply the Lippmann-Schwinger equation,
the obtained $K$ matrix will obey the same formula Eq. (18) with $k_{s}$
replaced by $k$. Consequently, the rescaled cross-section reaches maximum when
$a_{s}(0)=a_{{\rm res}}$. Figure 3 shows such comparison, where one can see
the SOC-induced energy shift that leads to $k_{s}\neq k$ is crucial to
chareterize the two-body scattering correctly, especially near the maximum of
$\lambda\sigma_{11}$.

Figure 3: Scattering results at zero scattering energy for the same Lennard-
Jones potential of Fig. 1. (a) Scaled partial cross-section
$\lambda\sigma_{11}$ as a function of $\lambda$. (b) $K$ matrix element
$\mathcal{K}_{11}$ as a function of $\lambda$. The blue solid curves are the
analytical results, and the red dashed curves are determined by the $s$-wave
only approximation Eq. (18), which are indistinguishable to the solid curves
on the scale shown. The purple dash-dotted curves are calculated using the
free-space pseudopotential directly. The green crosses are results from a
numerical calculation using the same Lennard-Jones potential with the presence
of SOC. (c) The blue solid curve shows $r_{{\rm vdW}}/[a_{s}(k_{s})-a_{{\rm
res}}]$ as a function of $\lambda$, while the purple dash-dotted curve shows
$r_{{\rm vdW}}/[a_{s}(0)-a_{{\rm res}}]$. Figure 4: Scattering results at zero
scattering energy for the same Lennard-Jones potential of Fig. 2. (a) Scaled
partial cross-section $\lambda\sigma_{11}$ as a function of $\lambda$. (b) $K$
matrix element $\mathcal{K}_{11}$ as a function of $\lambda$. The blue solid
curves are the analytical results, and the yellow dashed curves are determined
by the $s$-wave only approximation Eq. (18). The red crosses are results from
a numerical calculation using the same Lennard-Jones potential with the
presence of SOC.

Near $p$-wave resonances, the $p$-wave phase shift can no longer be neglected.
Neverthuless, a simplified fromula can be obtained in the low-energy limit
$k\rightarrow 0$, where $\tan[\delta_{s}(k_{s})]\rightarrow-
A_{s}\equiv-\tan[\delta_{s}(\sqrt{2}\lambda)]$ and
$\tan[\delta_{p}(k_{p})]\rightarrow-A_{p}k^{2}$, and the $K$ matrix is given
by,

$\lim_{k\rightarrow
0}\underline{\mathcal{K}}^{(+)}=\frac{1}{d}\left[\begin{array}[]{cc}b_{11}&b_{12}k/\lambda\\\
b_{21}k/\lambda&b_{22}k^{2}/\lambda^{2}\end{array}\right].$ (19)

Here, $d$ and $b_{\tau^{\prime}\tau}$ are all constants, which is given by
$d=1-2(\log\sqrt{2})A_{s}/\pi-2\lambda^{2}A_{p}/\pi+4\lambda^{2}A_{p}A_{s}[(\log\sqrt{2})-1]/\pi^{2}$,
$b_{11}=-A_{s}(1/2-\lambda^{2}A_{p}/\pi)$,
$b_{12}=b_{21}=A_{s}(1/2+\lambda^{2}A_{p}/\pi)$ and
$b_{22}=-A_{s}/4-\lambda^{2}A_{p}-\lambda^{2}A_{p}A_{s}(3/2-\log 2)/\pi$. When
$|\lambda^{2}A_{p}|\gtrsim 1$, $p$-wave scattering gives a significant
contribution as shown in Fig. 4, where Eq. (18) is no longer valid.
Neverthuless, the threshold laws for cross-section and $K$ matrix elements are
valid for all situations as shown in Figs. 1 and 2. Interestingly, the elastic
scattering rate $\propto k_{1}\sigma_{11}$ that determines thermalization
remains constant in the zero-energy limit, in contrast to the vanishing
$1/(\log k)^{2}$ rate without the presence of SOC. We also remark here that,
using the free-space pseudopotential including the $p$-wave contribution ${\rm
diag}[V_{1}^{{\rm fs}}(k,\rho),V_{0}^{{\rm fs}}(k,\rho),V_{1}^{{\rm
fs}}(k,\rho)]$ will wrongly give a vanishingly small $K$ matrix due to a
$\log(s)|_{s\rightarrow 0}$ term in the denumerator of all the matrix elements
of $\underline{\mathcal{K}}^{(+)}$, reflecting the importance of the non-
diagonal terms in the pseudopotnetial.

In summary, we have derived a pseudopotential in the COM frame with the
presence of SOC in 2D using a frame-transformation approach. Different than
the free-space pseudopotential, the s -wave scattering phase-shift changes due
to a SOC-induced energy shift. The frame-transformation also introduces non-
diagonal terms, which are also essential to reproduce two-body scattering
observables. We applied this pseudopotential with the Lippmann-Schwinger
equation to obtain the analytical scattering matrix and compare it with a
numerical scattering calculation with finite-range potential. Our
pseudopotential is valid even near $s$\- or $p$-wave resonances as long as
$\lambda r_{{\rm vdW}}\ll 1$, which is usually well satisfied in ultracold
quantum gases. Our results indicate that, if we consider, if we consider
$s$-wave only (which usually implies near $s$-wave resonances), and the
energy-dependency of $a_{s}(k)$ is very weak (which usually implies a very
broad resonance) so that $a_{s}(\sqrt{2}\lambda)\approx a_{s}(0)$, the free-
space pseudopotential can give a good approximation, which gives the valid
regime of previous studies in Refs. (Zhang  _et al._ , 2012c, 2013b). On the
other hand, if the energy-dependency of $a_{s}(k)$ is strong or $p$-wave
interaction is nonnegligible, our pseudopotential has to be adopted to
reproduce two-body scattering. Our approach can also be easily applied in 3D
and reproduce Eq. (11) of Ref. (Guan and Blume, 2017), which we will pursuit
elsewhere. Our results are also useful for investigating universal relations
and Tan’s contacts for SOC quantum gases in 2D (Cai-Xia Zhang, 2019) and might
eventually be applied in many-body physics studies.

## Supplemental Material

###  .1 Full analytical expression of $K$ matrix

Here, we give the full analytical expression of
$\underline{\mathcal{K}}^{(+)}$,

$\underline{\mathcal{K}}^{(+)}=\left[\begin{array}[]{cc}\mathcal{K}_{11}&\mathcal{K}_{12}\\\
\mathcal{K}_{21}&\mathcal{K}_{22}\end{array}\right]=\frac{1}{D}\left[\begin{array}[]{cc}B_{11}&B_{12}\\\
B_{21}&B_{22}\end{array}\right],$ (20)

where

$\begin{gathered}D=1+\frac{2}{\pi}t_{s}\left(\log\frac{k}{k_{s}}+\frac{\lambda}{k_{b}}\tanh^{-1}\frac{\lambda}{k_{b}}\right)+\frac{2}{\pi}t_{p}\left(\frac{\lambda^{2}}{k^{2}}-\frac{\lambda}{k_{b}}\tanh^{-1}\frac{\lambda}{k_{b}}\right)+\frac{4}{\pi^{2}}t_{s}t_{p}\left[\log
k_{2}\log k_{3}-(\log k)^{2}\right]\\\
+\frac{4}{\pi^{2}}t_{s}t_{p}\left[\frac{\lambda^{2}}{k^{2}}\left(\log\frac{k}{k_{s}}+\frac{\lambda}{k_{b}}\tanh^{-1}\frac{\lambda}{k_{b}}-1\right)+\left(2-\log\frac{k}{k_{s}}\right)\left(\frac{\lambda}{k_{b}}\tanh^{-1}\frac{\lambda}{k_{b}}\right)\right],\end{gathered}$
(21)
$B_{11}=t_{s}\frac{k_{1}}{2k_{b}}+t_{p}\frac{k_{2}}{2k_{b}}+\frac{t_{s}t_{p}}{\pi}\left(\frac{\lambda}{k_{b}}+\frac{\lambda}{k_{2}}-\frac{\lambda}{k_{b}}\log\frac{k}{k_{s}}-\log\frac{kk_{s}}{k_{2}^{2}}\right),$
(22)
$B_{12}=-t_{s}\frac{k}{2k_{b}}+t_{p}\frac{k}{2k_{b}}+\frac{2}{\pi}t_{s}t_{p}\frac{k}{2k_{b}}\left(\log\frac{k}{k_{s}}+\frac{\lambda^{2}}{k^{2}}\right),$
(23)

$B_{21}=-t_{s}\frac{k}{2k_{b}}+t_{p}\frac{k}{2k_{b}}+\frac{2}{\pi}t_{s}t_{p}\frac{k}{2k_{b}}\left(\log\frac{k}{k_{s}}+\frac{\lambda^{2}}{k^{2}}\right),$
(24)

$B_{22}=t_{s}\frac{k_{2}}{2k_{b}}+t_{p}\frac{k_{1}}{2k_{b}}-\frac{t_{s}t_{p}}{\pi}\left(\frac{\lambda}{k_{b}}+\frac{\lambda}{k_{1}}-\frac{\lambda}{k_{b}}\log\frac{k}{k_{s}}+\log\frac{kk_{s}}{k_{1}^{2}}\right),$
(25)

with the notations $t_{s}\equiv\tan\left[\delta_{s}(k_{s})\right]$,
$t_{p}\equiv\tan[\delta_{p}(k_{p})]$, $k_{s}=\sqrt{k^{2}+2\lambda^{2}}$,
$k_{p}=k$, $k_{b}=\sqrt{k^{2}+\lambda^{2}}$, $k_{1}=k_{b}+\lambda$ and
$k_{2}=k_{b}-\lambda$.

###  .2 Numerical method

We carry out a numerical calculation with finite range Lennard-Jones
potentials $U(\rho)$ to verify our analytical results, using a similar
procedure as Ref. (Wang  _et al._ , 2018). We expand the rescaled Hamiltonian
$h\equiv\rho^{1/2}H_{{\rm rel}}\rho^{-1/2}$ with the channel functions that
leads to a set of coupled differential equations for
$\underline{u}=\rho^{1/2}\underline{R}^{{\rm SO}}$:
$\underline{h}\underline{u}=E\underline{u}$, where
$\underline{h}=\underline{h}^{{\rm fs}}+\underline{V}^{{\rm SO}}$ has matrix
elements

$h_{\nu^{\prime}\nu}^{{\rm
fs}}=\left[-\frac{\hbar^{2}}{\mu_{2b}}\left(\partial_{\rho}^{2}-\frac{m_{\ell}^{2}-\frac{1}{4}}{\rho^{2}}\right)+U(\rho)\right]\delta_{\nu\text{\textquoteright}\nu}$
(26)

and

$V_{\nu^{\prime}\nu}^{{\rm
SO}}=\frac{\hbar^{2}k_{\mathrm{soc}}}{2\mu_{2b}}\left[\Sigma_{\chi^{\prime}\chi}^{(+)}p_{m_{\ell}^{\prime}m_{\ell}}^{(-)}+\Sigma_{\chi^{\prime}\chi}^{(-)}p_{m_{\ell}^{\prime}m_{\ell}}^{(+)}\right],$
(27)

where $p_{m_{\ell}^{\prime}m_{\ell}}^{(\pm)}=\mp(\partial_{\rho}-1/2\rho\pm
m_{\ell}/\rho)\delta_{m_{\ell}^{\prime},m_{\ell}\pm 1}$. Here we define
$\Sigma_{\chi^{\prime}\chi}^{(\pm)}=\left\langle\chi^{\prime}\left|s_{1}^{\pm}-s_{2}^{\pm}\right|\chi\right\rangle/2$
with $s_{n}^{+}$ ( $s_{n}^{-}$ ) beign the raising (lowering) operator for the
spin state of the $n$’th particle, which can be obtained explicitly as

$\begin{array}[]{cc}\Sigma_{\chi^{\prime}\chi}^{(\pm)}=&\sum_{m_{1}m_{2}}C_{s_{1}m_{1};s_{2}m_{2}}^{Sm_{S}}(a_{1\pm}C_{s_{1}m_{1}\pm
1;s_{2}m_{2}}^{S^{\prime}m_{S}^{\prime}}\\\
&-a_{2\pm}C_{s_{1}m_{1};s_{2}m_{2}\pm
1}^{S^{\prime}m_{S}^{\prime}})/2\end{array}$ (28)

with $a_{n\pm}=\sqrt{s_{n}(s_{n}+1)-m_{n}(m_{n}\pm 1)}$ and
$C_{s_{1}m_{1};s_{2}m_{2}}^{Sm_{S}}$ being Clebsch-Gordan coefficients. For
numerical calculations, one can solve the multichannel radial Schrödinger
equations by propagating wavefunction matrix $\underline{u}$ or equivalently
the logarithmic derivative matrix
$\underline{\mathcal{L}}=\underline{u}^{\prime}\underline{u}^{-1}$ to a large
enough distances $\rho_{{\rm max}}\gg\rho_{0}$. The $K$ matrix can be obtained
via
$\underline{\mathcal{K}}=(\underline{\mathcal{L}}\underline{g}-\underline{g}^{\prime})^{-1}(\underline{\mathcal{L}}\underline{f}-\underline{f}^{\prime})$,
where $\underline{f}=\sqrt{\rho}\underline{F}$ and
$\underline{g}=\sqrt{\rho}\underline{G}$. The numerical results are shown as
symbols in all figures in the main text, showing perfect agreement in the
range of $k$ and $\lambda$ considered.

## References

  * Fermi (1934) E. Fermi, Nuovo Cimento 11, 157 (1934).
  * Huang and Yang (1957) K. Huang and C. N. Yang, Phys. Rev. 105, 767 (1957).
  * Bolda  _et al._ (2002) E. L. Bolda, E. Tiesinga, and P. S. Julienne, Phys. Rev. A 66, 013403 (2002).
  * Chin  _et al._ (2010) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
  * Lin  _et al._ (2011) Y.-J. Lin, K. Jiménez-García, and I. B. Spielman, Nature(London) 471, 83 (2011).
  * Stock  _et al._ (2005) R. Stock, A. Silberfarb, E. L. Bolda, and I. H. Deutsch, Phys. Rev. Lett. 94, 023202 (2005).
  * Derevianko (2005) A. Derevianko, Phys. Rev. A 72, 044701 (2005).
  * Idziaszek and Calarco (2006) Z. Idziaszek and T. Calarco, Phys. Rev. Lett. 96, 013201 (2006).
  * Roth and Feldmeier (2001) R. Roth and H. Feldmeier, Phys. Rev. A 64, 043603 (2001).
  * Bloch  _et al._ (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * Paredes  _et al._ (2004) B. Paredes, A. Widera, V. Murg, O. Mandel, S. Fölling, I. Cirac, G. V. Shlyapnikov, T. W. Hänsch, and I. Bloch, Nature(London) 429, 177 (2004).
  * Kinoshita  _et al._ (2004) T. Kinoshita, T. Wenger, and D. S. Weiss, Science 305, 1125 (2004).
  * Hadzibabic  _et al._ (2004) Z. Hadzibabic, P. Krüger, M. Cheneau, B. Battelier, and J. Dalibard, Science 305, 1125 (2004).
  * Cladé  _et al._ (2009) P. Cladé, C. Ryu, A. Ramanathan, K. Helmerson, and W. D. Phillips, Phys. Rev. Lett. 102, 170401 (2009).
  * Rath  _et al._ (2010) S. P. Rath, T. Yefsah, K. J. Günter, M. Cheneau, R. Desbuquois, M. Holzmann, W. Krauth, and J. Dalibard, Phys. Rev. A 82, 013609 (2010).
  * Nishida  _et al._ (2013) Y. Nishida, S. Moroz, and D. T. Son, Phys. Rev. Lett. 110, 235301 (2013).
  * Gao  _et al._ (2015) C. Gao, J. Wang, and Z. Yu, Phys. Rev. A 92, 020504 (2015).
  * Kanjilal and Blume (2006) K. Kanjilal and D. Blume, Phys. Rev. A 73, 060701 (2006).
  * Zhang  _et al._ (2014) J. Zhang, H. Hu, X.-J. Liu, and H. Pu, Annu. Rev. Cold At. Mol. 2, 81 (2014).
  * Zhai (2015) H. Zhai, Rep. Prog. Phys. 78, 026001 (2015).
  * Dalibard  _et al._ (2011) J. Dalibard, F. Gerbier, G. Juzeliūnas, and P. Öhberg, Rev. Mod. Phys. 83, 1523 (2011).
  * Goldman  _et al._ (2014) N. Goldman, G. Juzeliunas, P. Öhberg, and I. B. Spielman, Rep. Prog. Phys. 77, 126401 (2014).
  * Struck  _et al._ (2012) J. Struck, C. Ölschläger, M. Weinberg, P. Hauke, J. Simonet, A. Eckardt, M. Lewenstein, K. Sengstock, and P. Windpassinger, Phys. Rev. Lett. 108, 225304 (2012).
  * Cheuk  _et al._ (2012) L. W. Cheuk, A. T. Sommer, Z. Hadzibabic, T. Yefsah, W. S. Bakr, and M. W. Zwierlein, Phys. Rev. Lett. 109, 095302 (2012).
  * Wang  _et al._ (2012) P. Wang, Z.-Q. Yu, Z. Fu, J. Miao, L. Huang, S. Chai, H. Zhai, and J. Zhang, Phys. Rev. Lett. 109, 095301 (2012).
  * Zhang  _et al._ (2012a) J.-Y. Zhang, S.-C. Ji, Z. Chen, L. Zhang, Z.-D. Du, B. Yan, G.-S. Pan, B. Zhao, Y.-J. Deng, H. Zhai, S. Chen, and J.-W. Pan, Phys. Rev. Lett. 109, 115301 (2012a).
  * Qu  _et al._ (2013) C. Qu, C. Hamner, M. Gong, C. Zhang, and P. Engels, Phys. Rev. A 88, 021604 (2013).
  * Khamehchi  _et al._ (2017) M. A. Khamehchi, K. Hossain, M. E. Mossman, Y. Zhang, T. Busch, M. M. Forbes, and P. Engels, Phys. Rev. Lett. 118, 155301 (2017).
  * Huang  _et al._ (2016) L. Huang, Z. Meng, P. Wang, P. Peng, S.-L. Zhang, L. Chen, D. L. andQi Zhou, and J. Zhang, Nat. Phys. 12, 540 (2016).
  * Meng  _et al._ (2016) Z. Meng, L. Huang, P. Peng, D. Li, L. Chen, Y. Xu, C. Zhang, P. Wang, and J. Zhang, Phys. Rev. Lett. 117, 235304 (2016).
  * Cui (2012) X. Cui, Phys. Rev. A 85, 022705 (2012).
  * Zhang  _et al._ (2012b) P. Zhang, L. Zhang, and Y. Deng, Phys. Rev. A 86, 053608 (2012b).
  * Yu (2012) Z. Yu, Phys. Rev. A 85, 042711 (2012).
  * Zhang  _et al._ (2013a) L. Zhang, Y. Deng, and P. Zhang, Phys. Rev. A 87, 053626 (2013a).
  * Duan  _et al._ (2013) H. Duan, L. You, and B. Gao, Phys. Rev. A 87, 052708 (2013).
  * Wang and Greene (2015) S.-J. Wang and C. H. Greene, Phys. Rev. A 91, 022706 (2015).
  * Guan and Blume (2016) Q. Guan and D. Blume, Phys. Rev. A 94, 022706 (2016).
  * Wang  _et al._ (2018) J. Wang, C. R. Hougaard, B. C. Mulkerin, and X.-J. Liu, Phys. Rev. A 97, 042709 (2018).
  * Zhang  _et al._ (2012c) P. Zhang, L. Zhang, and W. Zhang, Phys. Rev. A 86, 042707 (2012c).
  * Zhang  _et al._ (2013b) L. Zhang, Y. Deng, and P. Zhang, Phys. Rev. A 87, 053626 (2013b).
  * Guan and Blume (2017) Q. Guan and D. Blume, Phys. Rev. A 95, 020702 (2017).
  * Cai-Xia Zhang (2019) K. J. Cai-Xia Zhang, Shi-Guo Peng, “Universal relations for spin-orbit-coupled fermi gases in two and three dimensions,” (2019), arXiv:1907.11433 .


# Linear Cosmological perturbations in $f(Q)$ Gravity

Shambel Sahlu1,2, Endalkachew Tsegaye2 1 Astronomy and Astrophysics Research
and Development Department,
Entoto Observatory and Research Centre.  2 Department of Physics, College of
Natural and Computational Science, Wolkite University, Wolkite, Ethiopia

###### Abstract

In this manuscript, we studied the accelerated expansion history of the
universe and the formations of large-scale structures using $f(Q)$ gravity
model. The expansion rate of the universe within distance modulus and redshift
has been performed in $f(Q)$ gravity. This work is also devoted to
investigating the linear cosmological perturbations in the $f(Q)$ gravity
model using the $1+3$ covariant formalism. The scalar and harmonic
decomposition techniques are applied to find the evolution equation. The
growth of matter density contrasts is analyzed with redshift $z$. Using the
nonmetricity modified gravity model which is $f(Q)=Q+\alpha Q^{n}$, the
accelerating expanding universe and the formations of the large-scale
structures have been explained. The quasi-static approximation technique is
applied to analyze the growth of matter contrasts. The growth of density
contrast is studied within the dust and radiation-dominated universe in the
presence of nonmetricity gravity.

PACS numbers

04.50.Kd, 98.80.Jk, 98.80.-k, 95.36.+x, 98.80.Cq

##  1 Introduction

The modification of General Relativity (GR) [1, 2, 3, 4] has been gained much
attention to explain the late-time cosmic accelerated expansion history of the
universe[5, 6, 7] without the need for the dark-energy scenario.The symmetric
teleparallel gravity so called $f(Q)$ gravity theory is one of a modified
theory gravity which has vanishing curvature and torsion connection. $f(Q)$
gravity is described by a nonmetricity scalar $Q$ and it represents one of the
geometrical equivalent formulations of GR[8, 9, 10]. Some bservational
constraints has been investigated in $f(Q)$ model using SNIa data. And the
$f(Q)$ model is a good candidate as $\Lambda$CDM to describes the accelerated
expansion of universe.

The cosmological solutions and growth index of matter perturbations in $f(Q)$
gravity has been studied in [9] the metric formalisim approach and here in
this manuscript, we study the linear cosmologica perturbations using 1+3
covariant formalisim as presented in [11, 12, 13, 14] to investigate the
evolution of structure and to compute how the gravitational forces causing
small perturbations to grow and eventually seed the formation of stars,
quasars, galaxies and clusters[15]. The theory of cosmological perturbations
has become a cornerstone of modern cosmology since it is the framework which
provides the link between the models of early Universe such as the
inflationary universe scenario which yield causal mechanisms for the
generation of density fluctuations and allows us to connect.

This paper mainly focus on to investigate the linear cosmological
perturbations in $f(Q)$ gravity model using the 1+3 covariant formalism. The
power law model $f(Q)=Q+\alpha Q^{n}$ in[9, 16] is considered to find the
density contrast through cosmic time. For the case of $f(Q)=Q$, the symmetric
teleparallel gravity equivalent to GR and the field equation also reduced to
GR. The growth of density contrast is investigated and clearly observe the
contribution of nonmetricity scalar to study the formation of large-scale as
presented in [9, 16, 17].

##  2 $f(Q)$ Gravity Cosmology

Curvetureless and torsion-free connection describes the symmetric teleparallel
equivalen of general relativity (STEGR) we call it concident general
relativity. The covariant form of the field equation of $f(Q)$ gravity
generates based on the fundamental object is the nonmetricity tensor given by

$Q_{cab}=\triangledown_{c}g_{ab},$ (1)

where $g_{ab}$ is the metric. The field equation for $f(Q)$ gravity obtained
from the action presented in[9]. We consider a modified gravity theory in
which the fundamental object is the nonmetricity scalar given by

$S=\int\frac{1}{2}f(Q)\sqrt{-g}d^{4}x+\int L_{m}\sqrt{-g}d^{4}x,$ (2)

where $f(Q)$ is an arbitrary function of the nonmetricity $Q$, $\sqrt{-g}$ is
the determinant of the metric $g_{ab}$ and $L_{m}$ is the matter lagrangian
density. It is also useful to introduce the superpotential

$4P^{c}_{ab}=-Q^{c}_{ab}+2Q^{c}_{(ab)}-Q^{c}g_{ab}-\tilde{Q}^{c}g_{ab}-\delta^{c}_{(a}{(Q_{b})}.$
(3)

Varying the action with respect to the metric, one obtains the gravitational
field equation given by

$\frac{2}{\sqrt{-g}}\triangledown_{c}(\sqrt{-g}f^{\prime}P^{c}\hskip
0.28436pt_{ab})+\frac{1}{2}g_{ab}f+f_{Q}(P_{abd}Q_{b}\hskip
0.28436pt^{cd}-2Q_{cda}P^{c\beta}\hskip 0.28436pt_{b})=-T_{ab}.$ (4)

The field equation in (4) we can also rewrite as

$f^{\prime}G_{ab}+\frac{1}{2}g_{ab}(f^{\prime}Q-f)+2f^{\prime\prime}\nabla_{c}Q^{c}_{ab}=T_{ab},$
(5)

where $g_{ab}$ is metric, $f^{\prime\prime}=\frac{d^{2}f}{dQ^{2}}$,
$f^{\prime}=\frac{df}{dQ}$ and $T_{ab}=-\frac{2}{\sqrt{-g}}\frac{\delta
S_{m}}{\delta g_{ab}}$ is the matter field of EMT. The energy momentum tensor
of a perfect fluid is defind as

$T_{ab}=(\rho+p)u_{a}u_{b}+pg_{ab},$ (6)

where $u_{a}$ is the four-velocity satisfying the normalization condition
$u_{a}u_{b}=-1$, $\rho$ and p are the energy density and pressure of a perfect
fluid respectively. Under the homogeneous and isotropic universe described by
the Friedmann-Lemaitre-Robertson-Walker metric

$ds^{2}=-dt^{2}+a^{2}(t)g_{ab}dx^{a}dx^{b},$ (7)

where a(t) is the scale factor of the universe.

From the above modified field equation (5), the modified Friedmann equations
is presented as[9] for the case of $f(Q)=Q+F(Q)$, and its corresponding
thermodynamic quantities yield as Where $\rho$ and p are the density and
pressure respectively given by

$\displaystyle\rho_{eff}=\rho_{m}+\rho_{Q},$ (8) $\displaystyle
p_{eff}=p_{m}+p_{Q}.$ (9)

And the corresponding nonmetricity energy density and pressure become

$\displaystyle\rho_{Q}=\frac{F}{2}-QF^{\prime},$ (10) $\displaystyle
p_{Q}=-\rho_{Q}+\frac{\theta\dot{Q}}{3Q}(2QF^{\prime\prime}+F^{\prime}).$ (11)

and the modified Friedmann equations given as

$\displaystyle H^{2}=\frac{1}{3}(\rho_{m}+\frac{F}{2}-QF^{\prime}),$ (12)
$\displaystyle 2\dot{H}+3H^{2}=-\frac{\rho_{eff}}{3}-p_{eff}.$ (13)

From equations (10), (11), (12) and (13) we can write the energy density and
pressure of the universe respectively as follow:

$\displaystyle\frac{\dot{H}}{H^{2}}=-\dfrac{3}{2\mathcal{X}}+\frac{1}{2}\left(-\frac{\Omega_{m}}{\mathcal{X}}+\frac{\Re}{3\mathcal{X}}\right),$
(14)

where
$\mathcal{X}=(1+2QF^{\prime\prime}+F^{\prime}),\qquad\Omega_{m}=\frac{\rho}{3H^{2}},\qquad\Re=\frac{F-2QF^{\prime}}{6H^{2}}$.
In this manuscript, we we consider the power-law $F(Q)=\alpha Q^{n}$ gravity
model, and $\alpha$ becomes

$\alpha=\frac{(1-\Omega_{m})6H_{0}^{2}}{Q^{n}(1-n)}.$ (15)

The normalize Hubble parameter $h(z)$ in the presence of redshift the power-
law $f(Q)$ model expressed as

$\displaystyle
h(z)=\sqrt{(1+z)^{3}\Omega_{m}+\frac{(1-\Omega_{m})(1-2n)}{1-n}}.$ (16)

For the case of $n=0$ equation (16) reduced to the $\Lambda CDM$ model as
presented [18]

$h(z)=\sqrt{(1+z)^{3}\Omega_{m}+\Omega_{\Lambda}},\quad$ (17)

where $h=\frac{H}{H_{0}}$ and $\frac{a_{0}}{a}=1+z$, $a_{0}$ is the scale
factor at the present time which we can normalize to 1. The best-known way to
trace the evolution of the universe observationally is to look into the
redshift - luminosity distance relation. Now we can calculate the distance
modulus $\mu$ for (16) is given by

$\mu=m-M=25-5\times\log_{10}\left\\{3000\bar{h}^{-1}(1+z)\int_{0}^{z}\frac{dz^{\prime}}{h(z^{\prime})}\right\\},$
(18)

where apparent magnitude m and absolute magnitude M are logarithmic measure of
flux and luminosity respectively as well as the line-of-sight comoving
distance takes a value, $f(z)=3000\bar{h}^{-1}km/sMPc$ and the Hubble
uncertainty parameter is given interms of Hubble parameter by
$\bar{h}=\frac{H(z)}{100km/sMPc}$.

Figure 1: Some cosmological constraints in $f(Q)$ gravity model using SNIa
data. We used 100 random walkers and 50000 iterations.

Here we us MCMC simulation fits at the average values of free constrained
parameters we had the matter density, $\Omega_{m}$ and for the Hubble
uncertainty parameter, $\bar{h}$ and the free parameter $n$ as presented in
Fig. 1. From the plots we observe that $f(Q)$ gravity to explore the late-
history of the universe as the well known model $\Lambda$CDM.

##  3 Linear Cosmological Perturbations in $f(Q)$ Gravity

The Raychaudhuri equation of the fluid for GR limits as presented in [12, 13,
14, 17]

$\dot{\theta}=-\frac{1}{3}\theta^{2}-\frac{1}{2}(\rho+3p)+\nabla^{a}\dot{u}_{a}.$
(19)

For the case of $f(Q)$ gravity model Eq. (19) yields as

$\dot{\theta}=-\frac{1}{3}\theta^{2}-\frac{1}{2}(1+3w)\rho_{m}-\frac{1}{2}\left(-2\rho_{Q}+\frac{\theta\dot{Q}}{Q}(2QF^{\prime\prime}+F^{\prime})\right)+\nabla^{a}\dot{u}_{a}.$
(20)

This is a key equation to find the evolution equations. At the cosmological
background level we work on a homogeneous and isotropic expanding (FLRW)
univerese based on spatial gradients of gauge-invariant variables such as
$D^{m}_{a}$ and $Z_{a}$ represents energy density and volume expansion of the
fluid respectively are the basic tools to extract the evolution equation for
matter density fluctuations shows as

$\displaystyle D^{m}_{a}=\frac{a}{\rho_{m}}\tilde{\nabla}_{a}\rho_{m},$ (21)
$\displaystyle Z_{a}=a\tilde{\nabla}_{a}\theta.$ (22)

Based on the nonmetricity fluid for $f(Q)$ gravity by introducing new terms
$\mathcal{W}_{a}$ and $\mathcal{L}_{a}$ resulting from spatial gradients of
gauge-invariant quantities to characterize the fluctuations in the
nonmetricity density and momentum respectively defined as

$\displaystyle\mathcal{W}_{a}=a\tilde{\nabla}_{a}Q,$ (23)
$\displaystyle\mathcal{L}_{a}=a\tilde{\nabla}_{a}\dot{Q}.$ (24)

The conservation equations for non interacting fluid are given by

$\displaystyle\dot{\rho}_{m}=-\theta(\rho_{m}+p_{m})+(\rho_{m}+p_{m})\tilde{\nabla}^{a}\Psi_{a}$
(25)
$\displaystyle(\rho_{m}+p_{m})\dot{u}_{a}+)\tilde{\nabla}_{a}p_{m}+\dot{\Psi}_{a}-(3c^{2}_{s}-1)\frac{\theta}{3}\Psi_{a}+\Pi_{a}=0,$
(26)

where

$\Psi_{a}=\frac{q_{a}}{(\rho_{m}+p_{m})},\hskip
28.45274pt\Pi_{a}=\frac{\tilde{\nabla}^{b}\pi_{ab}}{(\rho_{m}+p_{m})}.$ (27)

And also $\Psi_{a}$ and $\Pi_{a}$ is vanish for perfect fluid, in our case we
consider the perfect fluid. From equations (21) and (22) we can get the first
order matter fluid perturbation evolution equations takes a form

$\displaystyle\dot{D}^{m}_{a}=-(1+w)Z_{a}+w\theta D^{m}_{a}\;,$ (28)
$\displaystyle\dot{Z}_{a}=-\left[\frac{2\theta}{3}+\frac{\dot{Q}F^{\prime}}{2Q}-\dot{Q}F^{\prime\prime}\right]Z_{a}+\left[\frac{w\theta^{2}}{3(1+w)}+\frac{w(1+3w)}{2(1+w)}\rho_{m}-\frac{wF}{2(1+w)},\right.$
$\displaystyle\left.-\frac{QF^{\prime}w}{1+w}+\frac{\theta\dot{Q}F^{\prime\prime}w}{1+w}+\frac{\theta\dot{Q}F^{\prime}w}{2Q(1+w)}-\frac{1}{2}(1+3w)\rho_{m}-\tilde{\nabla}^{2}\frac{w}{1+w}\right]D^{m}_{a}+\left[\frac{1}{2}F^{\prime},\right.$
$\displaystyle\left.-QF^{\prime\prime}-F^{\prime}+\theta\dot{Q}F^{\prime\prime\prime}-\frac{\theta\dot{Q}F^{\prime\prime}}{2Q}+\frac{\theta\dot{Q}F^{\prime}}{2Q^{2}}\right]\mathcal{W}_{a}-\left[\frac{\theta
F^{\prime}}{2Q}\mathcal{-}\theta F^{\prime\prime}\right]\mathcal{L}_{a}.$ (29)
$\displaystyle\dot{\mathcal{W}}_{a}=\mathcal{L}-\frac{w\dot{Q}}{1+w}D^{m}_{a},$
(30)
$\displaystyle\dot{\mathcal{L}}_{a}=\frac{\dddot{Q}}{\dot{Q}}\mathcal{W}-\frac{w\ddot{Q}}{1+w}D^{m}_{a}.$
(31)

By appling the scalar and harmonic decomposition technique as [12, 13, 14, 17,
19], the second-order evolution equaion in $f(Q)$ gravity is yield as

$\displaystyle\ddot{\Delta}_{m}^{k}=-\left[\frac{2\theta}{3}+\frac{\dot{Q}F^{\prime}}{2Q}+\dot{Q}F^{\prime\prime}-w\theta\right]\dot{\Delta}_{m}^{k}-\left[-wF-\theta\dot{Q}F^{\prime\prime}w+\frac{\theta\dot{Q}F^{\prime}w}{2Q},\right.$
$\displaystyle\left.-\frac{(1+3w)\rho_{m}}{2}(1-w)-\tilde{\nabla}^{2}w\right]\vartriangle_{m}^{k}-\left[\frac{1}{2}F^{\prime}-QF^{\prime\prime}-F^{\prime}+\theta\dot{Q}F^{\prime\prime\prime},\right.$
$\displaystyle\left.-\frac{\theta\dot{Q}F^{\prime\prime}}{2Q}+\frac{\theta\dot{Q}F^{\prime}}{2Q^{2}}\right](1+w)\mathcal{W}^{k}+\left[\frac{\theta
F^{\prime}}{2Q}-\theta F^{\prime\prime}\right](1+w)\dot{\mathcal{W}}^{k}.$
(32)
$\displaystyle\ddot{\mathcal{W}}^{k}=\frac{\dddot{Q}}{\dot{Q}}\mathcal{W}^{k}-\frac{2w\ddot{Q}}{1+w}\Delta_{m}-\frac{w\dot{Q}}{1+w}\dot{\Delta}_{m}^{k}\;,$
(33)

where

$k=\frac{2\pi a}{\lambda}\;,$ (34)

$k$ being the comoving wave-number and $\lambda$ the wavelength of the
perturbations. We apply the quasi- static approximation technique which is the
fluctuations of nonmetricity density $\mathcal{L}$ and $\mathcal{W}$ are
assumed to be constant with time, i.e.
$\dot{\mathcal{W}}^{k}=\ddot{\mathcal{W}}^{k}=\dot{\mathcal{L}}^{k}\approx 0$,
and the second order scalar linear equation reduced to

$\displaystyle\ddot{\vartriangle}_{m}^{k}=-\left[\frac{2\theta}{3}+\frac{\dot{Q}F^{\prime}}{2Q}-\dot{Q}F^{\prime\prime}-w\theta+\left(\frac{1}{2}F^{\prime}-QF^{\prime\prime}-F^{\prime}+\theta\dot{Q}F^{\prime\prime\prime}-\frac{\theta\dot{Q}F^{\prime\prime}}{2Q},\right.\right.$
$\displaystyle\left.\left.+\frac{\theta\dot{Q}F^{\prime}}{2Q^{2}}\right)\frac{w\dot{Q}^{2}}{\dddot{Q}}\right]\dot{\Delta}_{m}^{k}-\left[-wF-\theta\dot{Q}F^{\prime\prime}w+\frac{\theta\dot{Q}F^{\prime}w}{2Q}-\frac{(1+3w)\rho_{m}}{2}(1-w),\right.$
$\displaystyle\left.+\frac{k^{2}}{a^{2}}w+\left(\frac{1}{2}F^{\prime}-QF^{\prime\prime}-F^{\prime}+\theta\dot{Q}F^{\prime\prime\prime}-\frac{\theta\dot{Q}F^{\prime\prime}}{2Q}+\frac{\theta\dot{Q}F^{\prime}}{2Q^{2}}\right)\frac{2w\dot{Q}\ddot{Q}}{\dddot{Q}}\right]\vartriangle_{m}^{k}.$
(35)

As mentioned earlier, we admit here the power-law $f(Q)=Q^{n}$ gravity model
for illustrative purpose. We also use the the redshift-space transformation
technique as presented in [17, 14, 19], i.e.,for any time derivative functions
$f$ transform into a redshift derivative as $\dot{f}=H\frac{df}{dN}$, where
$N=\ln{a}$. Thereofore, the above second order evolution equation (35) for the
$f(Q)$ model reduced to

$\displaystyle\Delta^{\prime\prime}_{m}=\Bigg{[}2-3w-\frac{3}{2}(1+w)+\frac{\mathcal{G}n(1-\Omega_{m})}{(1-n)}-\frac{\mathcal{G}n(1-\Omega_{m})}{2}-\Big{(}\frac{(1-\Omega_{m})n}{2(1-n)}+(1-\Omega_{m})n$
$\displaystyle-\frac{(1-\Omega_{m})}{(1-n)}n-\mathcal{G}(1-\Omega_{m})n(n-2)+\frac{\mathcal{G}(1-\Omega_{m})n(n-2)}{2}+\frac{\mathcal{G}(1-\Omega_{m})n}{2(1-n)}\Bigg{)}\frac{2w}{3(1+w)^{2}}\Bigg{]}\frac{1}{(1+z)}{\Delta}^{\prime}_{m}$
$\displaystyle-\Bigg{[}-w\frac{(1-\Omega_{m})6}{(1-n)}+\frac{6\mathcal{G}(1-\Omega_{m})}{(1-n)}nw+\frac{3\mathcal{G}(1-\Omega_{m})}{(1-n)}n-\frac{(1+3w)3\Omega_{m}}{2}(1-w)+\frac{k^{2}}{H^{2}a^{2}}w$
$\displaystyle+\Bigg{(}\frac{(1-\Omega_{m})n}{2(1-n)}+(1-\Omega_{m})n-\frac{(1-\Omega_{m})}{(1-n)}n-\mathcal{G}(1-\Omega_{m})n(n-2)+\frac{\mathcal{G}(1-\Omega_{m})n}{2}w$
$\displaystyle+\frac{\mathcal{G}(1-\Omega_{m})n}{2(1-n)}\Bigg{)}\frac{4w}{(1+w)^{2}}\Bigg{]}\frac{1}{(1+z)^{2}}\Delta_{m}\;,$
(36)

where

$\mathcal{G}=\frac{-6-6\Omega_{m}}{6+6\Omega_{m}n+12n^{2}-12n^{2}\Omega_{m}}.$

Eq. (35) is the a closed system equation and we compute the growth of density
contrast from this equation in different approches. For the case of $f(Q)=Q$,
equation (35) reduced to $\Lambda CDM$ model as presented in [17, 14, 19]. We
admit the defined normalized energy density contrast for matter fluid as
presented in [17] as

$\delta^{k}(z)=\frac{\Delta_{m}^{k}(z)}{\Delta_{in}}.$ (37)

where the subscript in refers initial value of $\Delta_{m}(z)$ at redshift
$z_{in}$. For GR also we have

$\delta_{GR}=\frac{\Delta_{m}^{k}(z)(n=0)}{\Delta_{in}(z_{in})}.$ (38)

The normalized energy density can be defined as presented in [17]

$\delta_{[m,ch]}(z)=\frac{\Delta_{[m,ch]}(z)}{\Delta(z_{in})}\;,$ (39)

where $\Delta_{in}$ is the initial value of $\Delta_{[m,ch]}(z)$ at
$z_{in}=1100$, since the variation of CMB temperature detected observational
in the order of $10^{-3}$ [20] at $z\approx 1100$.

##  4 The growth density contrast in $f(Q)$ gravity

In this section, we are going to present the solutions and numerical results
of the density contrast ($\delta(z)$) from Eq. (42) in dust ($w=0$) and
radiation ($w=\frac{1}{3}$) dominanted univerese.

###  4.1 Dust Dominated Universe

The evolution equation in the presence of dust fluid and nonmetricity scalar
as

$\displaystyle{\Delta}^{\prime\prime}_{d}=\frac{1}{1+z}\left(\frac{1}{2}+\frac{\mathcal{G}n(1-\Omega_{m})}{(1-n)}-\frac{\mathcal{G}n(1-\Omega_{m})}{2}\right){\Delta^{\prime}}_{d}+\frac{3\Omega_{m}}{2(1+z)^{2}}\Delta_{d}.$
(40)

The solutions of the above Eq. (41) are yield as

$\Delta(z)=C_{1}\left(1+z\right)^{\frac{\alpha}{2}+\frac{1}{2}+\frac{\sqrt{\alpha^{2}+2\alpha+4\beta+1}}{2}}+{C_{2}}\left(1+z\right)^{\frac{\alpha}{2}+\frac{1}{2}-\frac{\sqrt{\alpha^{2}+2\alpha+4\beta+1}}{2}}$
(41)

where

$\alpha=\frac{1}{2}+\frac{Gn\left(1-\Omega\right)}{1-n}+\frac{Gn\left(1-\Omega\right)}{2}\qquad\mbox{and}\qquad\beta=\frac{2}{3}$

The numerical results of Eq. (41) is presented in Fig. 2. For illustrative
purpose, we use the constrained valuse of $n=1.0^{+07}_{-07}$ from Fig. 1.
From the plots at $n=0$, the results exactly reduced to GR limit as presented
in [17].

Figure 2: The growth of matter contrasts in the dust dominated universe in
$f(Q)$ gravity.

###  4.2 Radiation Dominated Universe

The evolution equation in the presence of radiation fluid and nonmetricity
scalar as

$\displaystyle\Delta^{\prime\prime}_{m}=\Bigg{[}-\frac{3}{2}+\frac{\mathcal{G}n(1-\Omega_{m})}{(1-n)}-\frac{\mathcal{G}n(1-\Omega_{m})}{2}-\Bigg{(}\frac{(1-\Omega_{m})n}{2(1-n)}+(1-\Omega_{m})n$
$\displaystyle-\frac{(1-\Omega_{m})}{(1-n)}n-\mathcal{G}(1-\Omega_{m})n(n-2)+\frac{\mathcal{G}(1-\Omega_{m})n(n-2)}{2}+\frac{\mathcal{G}(1-\Omega_{m})n}{2(1-n)}\Bigg{)}\frac{1}{8}\Bigg{]}\frac{1}{(1+z)}{\Delta}^{\prime}_{m}$
$\displaystyle-\Bigg{[}-\frac{(1-\Omega_{m})6}{3(1-n)}+\frac{6\mathcal{G}(1-\Omega_{m})}{3(1-n)}n+\frac{3\mathcal{G}(1-\Omega_{m})}{(1-n)}n-\frac{4\Omega_{m}}{3}+\frac{k^{2}}{3H^{2}a^{2}}$
$\displaystyle+\Bigg{(}\frac{(1-\Omega_{m})n}{2(1-n)}+(1-\Omega_{m})n-\frac{(1-\Omega_{m})}{(1-n)}n-\mathcal{G}(1-\Omega_{m})n(n-2)+\frac{\mathcal{G}(1-\Omega_{m})n}{6}$
$\displaystyle+\frac{\mathcal{G}(1-\Omega_{m})n}{2(1-n)}\Bigg{)}\frac{3}{4}\Bigg{]}\frac{1}{(1+z)^{2}}\Delta_{m}\;.$
(42)

###  4.3 Long Wavelength Mode

For the case of longwave length, $\frac{k^{2}}{3a^{2}H^{2}}<<1$, all
cosmological fluctuations begin and remain inside the Hubble horizon. The
solutions of Eq. (42) in longwave mode yield as

$\displaystyle\Delta_{r}=C_{1}\left(1+z\right)^{-\frac{l}{16}+\frac{Gn\left(1-\Omega\right)}{2-2n}-\frac{Gn\left(1-\Omega\right)}{4}+\Gamma}+C_{2}\left(1+z\right)^{-\frac{l}{16}+\frac{Gn\left(1-\Omega\right)}{2-2n}-\frac{Gn\left(1-\Omega\right)}{4}-\Gamma}$
(43)

where

$\displaystyle\Gamma=\frac{\sqrt{\Psi}}{16}\;,$
$\displaystyle\Psi=l^{2}-16l\left(-1+\frac{Gn\left(1-\Omega\right)}{1-n}-\frac{Gn\left(1-\Omega\right)}{2}\right)+64\left(-1+\frac{Gn\left(1-\Omega\right)}{1-n}-\frac{Gn\left(1-\Omega\right)}{2}\right)^{2}$
$\displaystyle-256\eta-16l-64+\frac{128Gn\left(1-\Omega\right)}{1-n}-64Gn\left(1-\Omega\right)-192\zeta$
$\displaystyle\eta=-2\Omega-\frac{2\left(1-\Omega\right)}{1-n}+\frac{5Gn\left(1-\Omega\right)}{1-n}\;,$
$\displaystyle
l=\frac{n\left(1-\Omega\right)}{2-2n}+n\left(1-\Omega\right)-\frac{\left(1-\Omega\right)n}{1-n}-\frac{G\left(1-\Omega\right)n\left(n-2\right)}{2}+\frac{Gn\left(1-\Omega\right)}{2-2n}\;,$
$\displaystyle\zeta=\frac{n\left(1-\Omega\right)}{2-2n}+n\left(1-\Omega\right)-\frac{\left(1-\Omega\right)n}{1-n}-G\\!\left(1-\Omega\right)n\left(n-2\right)+\frac{Gn\left(1-\Omega\right)}{2}+\frac{Gn\left(1-\Omega\right)}{2-2n}\;.$

Then the numerical results os density contrast for longwave mode is presented
in Fig Eq. 3 from Eq. (43). From the plots we observe the amplitude of density
contrast is decaying through redshift.

Figure 3: The growth of energy contrasts in the radiation dominated universe
in longwave length mode for different values of $n$.

###  4.4 Short Wavelength Mode

For the case of short wavelength mode $\frac{k^{2}}{3a^{2}H^{2}}>>1$, the
solutions of Eq. (42) yield as

$\displaystyle\Delta_{r}(z)=C_{3}\left(1+z\right)^{\frac{1}{2}+\frac{\alpha}{2}}BesselJ\Bigg{(}{\frac{\sqrt{\alpha^{2}+2\alpha+4\beta+1}}{4}},\frac{\sqrt{3}\,\pi}{3\lambda\left(1+z\right)^{2}}\Bigg{)}$
(44)
$\displaystyle+C_{4}\left(1+z\right)^{\frac{1}{2}+\frac{\alpha}{2}}BesselY\left({\frac{\sqrt{\alpha^{2}+2\alpha+4\beta+1}}{4}},\frac{\sqrt{3}\,\pi}{3\lambda\left(1+z\right)^{2}}\right)$

From Eq.(34), we rewrite as
$\frac{k^{2}}{3a^{2}H^{2}}=\frac{16\pi^{2}}{3\lambda^{2}(1+z)^{4}}$ as
presented in [17], and the numerical results of the above solutions are
presented in Figs. 5 and 5 for different values of wavelength $\lambda$. From
the plot. we can see the amplitude of density contrast is oscillates in the
short wavelength. for different values of $n$.

Figure 4: The growth of matter contrasts in the radiation dominated universe
in short wave length.

Figure 5: The growth of matter contrasts in the radiation dominated universe
in short wave length.

##  5 Summary

In this work, we considered the symmetric teleparallelism, it is curvatureless
and torsion-free gravity and fully described by nonmetricity scalar $Q$ and
studied the modified gravity in the regime of most general class of symmetric
teleparallel gravity theories whose action is quadratic in the nonmetricity
tensor. The analysis is done by using linear cosmological perturbations using
1+3 covariant gauge-invariant formalism. We have also considered the power law
$f(Q)$ gravity model $f(Q)=Q+\alpha Q^{n}$. In which the presented $f(Q)$
gravity model approach provides gravitational alternatives to the dark energy
to explain the late and early phase of universe. For both regimes, the
symmetric teleparallel gravity is reduced to GR at $f(Q)=Q$ gravity model.

## References

## References

  * [1] M. Elmardi, A. Abebe, and A. Tekola, “Chaplygin-gas solutions of $f(R)$ gravity,” International Journal of Geometric Methods in Modern Physics, vol. 13, no. 10, p. 1650120, 2016.
  * [2] S. Sahlu, J. Ntahompagaze, M. Elmardi, and A. Abebe, “The chaplygin gas as a model for modified teleparallel gravity?,” The European Physical Journal C, vol. 79, no. 9, pp. 1–31, 2019.
  * [3] H. B. Benaoum, “Modified chaplygin gas cosmology,” Advances in High Energy Physics, vol. 2012, 2012.
  * [4] S. Sahlu, J. Ntahompagaze, A. Abebe, and D. F. Mota, “Accelerating universe in modified teleparallel gravity theory,” Proceedings of the International Astronomical Union, vol. 15, no. S356, pp. 397–399, 2019.
  * [5] E. V. Linder, “Exploring the expansion history of the universe,” Physical Review Letters, vol. 90, no. 9, p. 091301, 2003.
  * [6] T. Harko, T. S. Koivisto, F. S. Lobo, G. J. Olmo, and D. Rubiera-Garcia, “Coupling matter in modified $f(Q)$ gravity,” Physical Review D, vol. 98, no. 8, p. 084043, 2018.
  * [7] B. S. Haridasu, V. V. Luković, M. Moresco, and N. Vittorio, “An improved model-independent assessment of the late-time cosmic expansion,” Journal of Cosmology and Astroparticle Physics, vol. 2018, no. 10, p. 015, 2018\.
  * [8] J. B. Jiménez, L. Heisenberg, T. Koivisto, and S. Pekar, “Cosmology in $f(q)$ geometry,” Physical Review D, vol. 101, no. 10, p. 103507, 2020\.
  * [9] W. Khyllep, A. Paliathanasis, and J. Dutta, “Cosmological solutions and growth index of matter perturbations in f (q) gravity,” Physical Review D, vol. 103, no. 10, p. 103521, 2021.
  * [10] R. Solanki, S. Pacif, A. Parida, and P. Sahoo, “Cosmic acceleration with bulk viscosity in modified $f(Q)$ gravity,” Physics of the Dark Universe, vol. 32, p. 100820, 2021.
  * [11] M. Bruni, P. K. Dunsby, and G. F. Ellis, “Cosmological perturbations and the physical meaning of gauge-invariant variables,” The Astrophysical Journal, vol. 395, pp. 34–53, 1992.
  * [12] A. Abebe, M. Abdelwahab, Á. De la Cruz-Dombriz, and P. K. Dunsby, “Covariant gauge-invariant perturbations in multifluid ${f(R)}$ gravity,” Class. Quant. Grav., vol. 29, no. 13, p. 135011, 2012.
  * [13] J. Ntahompagaze, A. Abebe, and M. Mbonye, “A study of perturbations in scalar–tensor theory using 1+ 3 covariant approach,” Int. J. Mod. Phys. D, vol. 27, no. 3, p. 1850033, 2018.
  * [14] H. Sami, S. Sahlu, A. Abebe, and P. K. Dunsby, “Covariant density and velocity perturbations of the quasi-newtonian cosmological model in $f(T)$ gravity,” The European Physical Journal C, vol. 81, no. 10, pp. 1–17, 2021\.
  * [15] J. N. Fry, “The galaxy correlation hierarchy in perturbation theory,” The Astrophysical Journal, vol. 279, pp. 499–510, 1984.
  * [16] I. Ayuso, R. Lazkoz, and V. Salzano, “Observational constraints on cosmological solutions of $f(Q)$ theories,” Physical Review D, vol. 103, no. 6, p. 063505, 2021.
  * [17] S. Sahlu, J. Ntahompagaze, A. Abebe, Á. de la Cruz-Dombriz, and D. F. Mota, “Scalar perturbations in $f(T)$ gravity using the $1+3$ covariant approach,” The European Physical Journal C, vol. 80, no. 5, pp. 1–19, 2020.
  * [18] M. Moresco, L. Pozzetti, A. Cimatti, R. Jimenez, C. Maraston, L. Verde, D. Thomas, A. Citro, R. Tojeiro, and D. Wilkinson, “A 6% measurement of the hubble parameter at z 0.45: direct evidence of the epoch of cosmic re-acceleration,” Journal of Cosmology and Astroparticle Physics, vol. 2016, no. 05, p. 014, 2016.
  * [19] J. Ntahompagaze, S. Sahlu, A. Abebe, and M. R. Mbonye, “On multifluid perturbations in scalar–tensor cosmology,” International Journal of Modern Physics D, vol. 29, no. 16, p. 2050120, 2020.
  * [20] G. F. Smoot et al., “Structure in the cobe differential microwave radiometer first-year maps,” The Astrophysical Journal, vol. 396, pp. L1–L5, 1992.