# Exploring scale invariance in the expansion of a spherical unitary Fermi gas

Lu Wang1,2, Xiangchuan Yan1, Jing Min1,2 Dali Sun1<EMAIL_ADDRESS>Xin
Xie1,2, Shi-Guo Peng1 Mingsheng Zhan1 Kaijun Jiang1<EMAIL_ADDRESS>1State
Key Laboratory of Magnetic Resonance and Atomic and Molecular Physics, Wuhan
Institute of Physics and Mathematics, Innovation Academy for Precision
Measurement Science and Technology, Chinese Academy of Sciences, Wuhan 430071,
China  2University of Chinese Academy of Sciences, Beijing 100049, China

###### Abstract

A unitary Fermi gas in an isotropic harmonic trap is predicted to show scale
and conformal symmetry that have important consequences in its thermodynamic
and dynamical properties. By experimentally realizing an isotropic harmonic
trap, we study the expansion of a unitary Fermi gas and demonstrate its
universal expansion dynamics along different directions and at different
temperatures. We show that as a consequence of SO(2,1) symmetry, the measured
release energy is equal to that of the trapping energy. In addition, away from
resonance when scale invariance is broken, we determine the effective exponent
$\gamma$ that relates the chemical potential and average density along the
BEC-BCS crossover, which qualitatively agrees with the mean field predictions.
This work opens the possibility of studying non-equilibrium dynamics in a
conformal invariant system in the future.

Strongly interacting Fermi gases are created by tuning the interaction
strength between atoms of different spin states via Feshbach resonance
Chin2010 . The unitary Fermi gas, realized when the $s$-wave scattering length
is tuned to infinity, is of special interest Kinast2005 ; Nascimbene2010 ;
Ku2012 ; Sidorenkov2013 ; Bardon2014 ; Patel2020 ; Li2022 , as it is not only
strongly correlated but also an example of scale-invariant quantum many-body
system. One of the basic tools used to explore the properties of unitary Fermi
gas is the expansion dynamics Menotti2002  and much insight has been obtained
about the role of interactions. For instance, anisotropic expansions OHara2002
, virial theorem Thomas2005 , conformal symmetry breaking Elliott2014 , and
Efimovian expansion Deng2016  of Fermi gases have been demonstrated by imaging
the expansion of the atomic cloud.

The strongly interacting Fermi gas at finite temperature is described by a
hydrodynamic theory (see Eq. 1), where the transport behaviors are determined
by viscosities Cao2011a ; Levin2011viscosity ; Shafer2017viscosity . At
unitarity, the bulk viscosity $\zeta_{B}$ vanishes, and the friction force
arise from shear viscosity $\eta$. For a unitary Fermi gas in an anisotropic
trap studied previously, the conformal symmetry is broken and the shear
viscosity plays an dominant role, which allowed its extraction from expansion
dynamics Cao2011 ; Thomas2014anomalousViscosity ;
Thomas2015superfluidViscosity . On the other hand, for a spherical unitary
Fermi gas, the transverse relative motion of the atomic cloud is absent
($\sigma_{ii}=0$), and consequently the effect of the shear viscosity can be
neglected. The system without viscosity contribution would have universal
properties regardless of the interaction details. Contrary to the anisotropic
system, the spherical unitary Fermi gas has a hidden SO(2,1) symmetry
Werner2006 , which is predicted to has several universal properties, such as
the exact relations between the trapping potential energy and total energy.
However, the preparation and exploration of the universal properties of a
spherical unitary Fermi gas are yet to be demonstrated experimentally.

In this letter, we produce a spherical Fermi gas in an optical dipole trap
(ODT) and study its expansion behaviors in the strongly interacting regimes.
By tuning the interaction strength to unitarity, the expansion of the system
shows the scale invariance along different directions and at different
temperatures, which is absent in an anisotropic system. We find that the
trapping potential energy equals to the half of the total energy of the
system, and verify the virial theorem for the unitary Fermi gas Werner2006 .
To the best of our knowledge, this is the first experiment on the 3D ultracold
quantum gases with the SO(2, 1) symmetry. Furthermore, we explore expansion
dynamics away from unitarity when scale invariance is broken, and measure the
effective exponent $\gamma$ of $\mu(n)\propto n^{\gamma}$ where $\mu$ is the
chemical potential and $n$ is the average density. The measured values of
$\gamma$ qualitatively agree with the mean-field calculations.

The expansion of strongly interacting Fermi gases is described by the
hydrodynamic theory Cao2011 ; Elliott2014 ; Supplemental2023 ,

$\begin{split}\frac{d^{2}}{dt^{2}}\frac{m{\langle}x_{i}^{2}{\rangle}}{2}=&{\langle}x_{i}\cdot{\frac{\partial
U}{\partial x_{i}}}{\rangle}_{0}+\frac{1}{N}\int{d^{3}r[\Delta p-\Delta
p_{0}]}\\\
&-\frac{1}{N}\int{d^{3}r(\eta\sigma_{ii}+\zeta_{B}\sigma^{\prime})},\end{split}$
(1)

where ${\langle}x_{i}^{2}{\rangle}$ represents the mean square cloud radius
along the $i$th axis ($i=x,y,z$), $U$ is the trapping potential, $t$ is the
expansion time, the subscript $(_{0})$ denotes the initial condition in the
trap at $t=0$, and $\Delta p=p-(2/3)\varepsilon$ is the scale-invariance
breaking pressure, where $\varepsilon$ is the energy density Ho2004 . The last
term on the right describes the friction forces arising from shear viscosity
$\eta$ and bulk viscosity $\zeta_{B}$. Here,
$\sigma_{ii}=2\dot{b}_{i}/b_{i}-(2/3)\sum_{j}\dot{b}_{j}/b_{j}$ represents the
transverse relative motion and $\sigma^{\prime}=\sum_{i}\dot{b}_{i}/b_{i}$ for
the dilation process, where $b_{i}$ denotes the expansion scale factor. In the
unitary regime, both $\Delta p$ and $\zeta_{B}$ vanish Cao2011 ; Son2007 ;
Escobedo2009 ; Dusling2013 . The value of $\sigma_{ii}$ depends on the
geometry or symmetry of the atomic cloud. For the system in an anisotropic
trap, $\sigma_{ii}$ is non zero. Only for a spherical gas, the relative motion
is absent with $\sigma_{ii}=0$, and in this case, we obtain the expansion
behavior,

$\begin{split}{\langle}x_{i}^{2}{\rangle}&={\langle}x_{i}^{2}{\rangle}_{0}+\frac{t^{2}}{m}{\langle}x_{i}\cdot{\frac{\partial
U}{\partial x_{i}}}{\rangle}_{0}.\end{split}$ (2)

Eq. (2) shows the ballistic expansion analogous to a non-interacting ideal
gas, and the interaction is included in the in-situ atomic cloud size
${\langle}x_{i}^{2}{\rangle}_{0}$.

The scale-invariant expansion can be tested by determining the value,

$\begin{split}\tau_{i}^{2}(t)=\frac{{\langle}x_{i}^{2}{\rangle}_{t}-{\langle}x_{i}^{2}{\rangle}_{0}}{{\langle}x_{i}^{2}{\rangle}_{0}\omega^{2}},\end{split}$
(3)

where $\omega_{x}=\omega_{y}=\omega_{z}=\omega$ is the trapping frequency of
the isotropic trap. According to Eq. (2), $\tau^{2}(t)=t^{2}$.

In the experiment, we prepare a spherical Fermi gas with variable trapping
frequencies and tunable interactions based on our previous works Yan2021 ;
Yan2022 . We initially prepare a 6Li atomic degenerate Fermi gas with two spin
states $|F=1/2,m_{F}=\pm 1/2\rangle$ in an elongated ODT and at the Feshbach-
resonance magnetic field of $834$ G. The experimental setup of the isotropic
trap is schematically displayed in Fig. 1(a), where some special techniques
are applied. Firstly, a magnetic field with a gradient $B^{\prime}_{z}=1.05$
G/cm is applied along $z$ axis to simultaneously compensate the gravity force
of the two spin states. This is valid for 6Li atoms because the hyperfine
interaction is much smaller than the Zeeman shift at the applied magnetic
field. Secondly, two elliptic optical beams with a cross-sectional aspect
ratio of $\sqrt{2}$, propagating perpendicularly in the horizontal plane, form
the isotropic trap. Under these conditions, the trapping frequency can be
varied by adjusting the optical power (See Supplemental Material
Supplemental2023  for details). We transfer the unitary Fermi gas to the
isotropic trap with an efficiency of more than ${90\%}$, by lowering the power
of the elongated ODT and simultaneously increasing that of the isotropic ODT
in a period of 25 ms. After performing the evaporative cooling in the
isotropic trap, we slowly increase the optical power to $3.8$ W in about 75
ms. The temperature is adjusted by controlling the depth of the evaporative
cooling. The atom number is $N=2.9(3)\times 10^{4}$. The trapping frequencies
are $(\omega_{x},\omega_{y},\omega_{z})=2\pi\times(1234(6),1165(11),1204(3))$
Hz, which are nearly the same along three axes Supplemental2023 .

Figure 1: (Colour online). Production of a spherical Fermi gas and scale-
invariant expansion along different directions. (a) Schematics of the
experimental setup. Two elliptic optical beams, propagating along $x$ and $y$
axis, respectively, form the isotropic trap. The ratio of the beam width along
$z$ direction to that in the horizontal plane is $\sqrt{2}$. A pair of anti-
Helmholtz coils produce a gradient magnetic field along vertical direction to
compensate the gravity. Another pair of Helmholtz coils (not shown) produce a
homogeneous magnetic field along vertical direction to tune the interactions.
(b) Exemplary atomic images are taken in vertical and horizontal directions,
respectively. The expansion time is 0.5 ms, 0.8 ms, 1.1 ms, 1.4 ms and 1.7 ms.
$y^{\prime}$ direction is between $x$ and $y$ axes. (c) Values of
$\tau_{i}^{2}(t)$ versus the expansion time $t$ along different directions
($i=x,y,z,y^{\prime}$). The error bar is the standard deviation of several
measurements. The black solid curve represents the scale theory
$\tau^{2}(t)=t^{2}$. $T/T_{\textrm{F}}=$ $0.36(3)$ where $T_{\textrm{F}}$ is
the Fermi temperature of the non-interacting Fermi gas.

We switch off the isotropic ODT and measure the cloud width versus the
expansion time $t$ at the magnetic field $B=834$ G. Two laser beams with a
frequency difference of 76 MHz, propagating along vertical and horizonal
directions, respectively, are applied to detect two spin states. Exemplary
atomic images during the expansion are shown in Fig. 1(b), indicating an
isotropic expansion in direct contradiction to an elongated Fermi gas
OHara2002 . We use a fringe-removal algorithm Supplemental2023 ; Ockeloen2010
to reduce the imaging noise, which is helpful in accurately determining the
size of the atomic cloud. The cloud radius ${\langle}x_{i}^{2}{\rangle}_{t}$
is obtained by fitting a Gaussian distribution to the atomic density profile.
In the unitary regime, the cloud radius in the trap
${\langle}x_{i}^{2}{\rangle}_{0}$ can be theoretically calculated, and the
temperature is determined by analyzing the atomic density distribution
Kinast2005 ; Yan2021 ; Supplemental2023 . At the temperature
$T/T_{\textrm{F}}=$ $0.36(3)$, values of $\tau^{2}(t)$ are calculated
according to Eq. (3). As shown in Fig. 1(c), the expansion behaviors along
different directions all obey the scale theory $\tau^{2}(t)=t^{2}$, which
indicates the absence of the effect of viscosity. This scale-invariant
expansion along different directions is unique for a spherical Fermi gas in
the unitary regime. For an anisotropic Fermi gas, only the sum of sizes along
three axes shows the scale-invariant expansion Elliott2014 .

Figure 2: (Colour online). Scale-invariant expansion at different
temperatures. (a) Atomic cloud size ${\langle}x^{2}(t){\rangle}$ versus the
expansion time $t$ at the unitary regime ($a\rightarrow\infty$).
$T/T_{\textrm{F}}=0.36(3)$ (black), $0.62(2)$ (red) and $0.81(2)$ (blue),
respectively. For comparison, the non-interacting Fermi gas ($a=0$) is also
shown with $T/T_{\textrm{F}}=$ $0.91(3)$ (green). The solid curves represent
the calculations with Eq. (2). The error bar is the standard deviation of
several measurements. (b) Values of $\tau^{2}(t)$ versus the expansion time
$t$ at different temperatures. The solid black curve denotes the scale theory
$\tau^{2}(t)=t^{2}$.

In Fig. 2, we measure the atomic expansion at different temperatures. Only
expansion along the $x$-axis is displayed for simplicity. Due to the finite-
temperature effect, the atomic cloud size ${\langle}x^{2}(t){\rangle}$ shows
an obvious difference, as shown in Fig. 2(a). System at a higher temperature
has a larger in-situ cloud radius ${\langle}x_{i}^{2}{\rangle}_{0}$, leading
to a faster expansion, which agrees well with the theoretical prediction of
Eq. (2). While values of $\tau^{2}(t)$ at different temperatures are
consistent, all obeying the scale theory $\tau^{2}(t)=t^{2}$ (see Fig. 2(b)).
For comparison, the expansion behavior of a non-interacting Fermi gas ($a=0$,
where $a$ is the $s$-wave scattering length) is also shown. The Fermi gas in
the unitary ($a=\infty$) regime has the same scaled expansion behavior with
that of the non-interacting ($a=0$) Fermi gas.

Figure 3: Verifying the virial theorem for a unitary Fermi gas. The trapping
potential $U$ shows a linear dependence on the release energy
$E_{\textrm{rel}}$. The error bar of the trapping potential comes from the
uncertainty of the trapping frequency. The error bar of the release energy
from the standard deviation of several measurements is smaller than the mask
size. Both energies are normalized by the Fermi energy
$E_{\textrm{F}}=(3N)^{1/3}\hbar\omega$, where $N$ is the atom number and
$\omega$ is the trapping frequency. The solid line represents the relation
$U=E_{\textrm{rel}}$.

According to the exact property predicted by the hidden SO(2,1) symmetry, the
total energy should be twice the trapping potential energy Werner2006 . This
is called the virial theorem for a unitary Fermi gas, which can be verified
using the expansion method. The total energy of the trapped gas is the sum of
trapping potential energy $U$, kinetic energy $E_{\textrm{kin}}$ and
interaction energy $E_{\textrm{int}}$, i.e.,
$E_{\textrm{tot}}=U+E_{\textrm{kin}}+E_{\textrm{int}}$. After switching off
the trapping potential ($U\rightarrow 0$), the release energy
$E_{\textrm{rel}}=E_{\textrm{kin}}+E_{\textrm{int}}$ remains constant during
the expansion process Stringari1999RMP ; Cooper1997PRL ;
liExpansionDynamicsSpherical2019  and will be completely converted to the
kinetic energy in the long-time expansion. So we only need to demonstrate the
equality relation $U=E_{\textrm{rel}}$. By fitting the slope of atomic cloud
radii respect to the expansion time, we obtain the release energy
$E_{\textrm{rel}}=(3/2)mv_{x}^{2}$, where $v_{x}$ is the expansion velocity
along $x$ axis. We can also determine the trapping potential energy
$U=(3/2)m\omega^{2}{\langle}x^{2}{\rangle}_{0}$. The experimental results are
shown in Fig. 3, where the atomic temperature varies across the Fermi
degeneracy. The trapping potential energy $U$ is equal to the release energy
$E_{\textrm{rel}}$ over a wide range of energies.

Away from Feshbach resonance, the scattering length is finite with $\Delta
p\neq 0$ and the scale invariance will be broken. We assume a power-law
dependence of the chemical potential, $\mu(n)=n^{\gamma}$, where $n$ is the
average atomic density and $\gamma=n/\mu\cdot d\mu/dn$ is the effective
exponent Stringari2008RMP ; huCollectiveModesBallistic2004 ; Heiselberg2004PRL
. By imposing that the total energy variation vanishes to first order, one
gets the energy relation in the BEC-BCS crossover Stringari2008RMP , $3\gamma
E_{\textrm{rel}}=2U$. Considering that the bulk viscosity is negligibly small
Elliott2014  and the effect of the shear viscosity in a spherical system is
zero, we obtain the solution of the expansion scale factors
huCollectiveModesBallistic2004 ,

$\ddot{b}_{i}-(\omega^{2}/b_{i})\Gamma^{-\gamma}=0,$ (4)

where $\Gamma=b^{3}_{i}$. Eq. (4) is a decoupled equation for each single
direction, where $b_{i}$ can be easily calculated if knowing the value of
$\gamma$. Especially in the unitary regime with $\gamma=2/3$, Eq. (4) has an
analytical solution the same as Eq. (2). In the BCS-BEC crossover,
$\tau^{2}(t)$ can be calculated as

$\displaystyle\tau^{2}(t)=\frac{b^{2}_{i}-1}{\omega^{2}}.$ (5)

To measure $\gamma$ at different interactions, we adiabatically ramp the
magnetic field from 834 G to the desired value within 300 ms. The power of the
optical trap is decreased to about 100 mW to perform the experiment at low
temperature, where the trapping frequencies are
$(\omega_{x},\omega_{y},\omega_{z})=2\pi\times(188(1),176(2),187(1))$ Hz. From
Eq. (4), the value of $\gamma$ can be obtained according to an iterative
method. We only consider the expansion along $x$ axis. The value of $\gamma$
at zero temperature huCollectiveModesBallistic2004  is initially input into
Eq. (4) to calculate $b_{x}(t)$. Then we can obtain the cloud radius in the
trap from the expansion data,
${\langle}x^{2}{\rangle}_{0}={\langle}x(t)^{2}{\rangle}_{t}/b(t)^{2}_{x}$,
which is used to determine the trapping potential energy $U$. We also measure
the release energy $E_{\textrm{rel}}$ from the long-time expansion. Using the
energy relation, we obtain a new value of $\gamma$ and input it into Eq. (4)
again for the next iterative calculation. We repeat the calculation until
$(\gamma_{i+1}-\gamma_{i})/\gamma_{i}\leq 10^{-5}$, where $i$ denotes the
number of iterations.

Figure 4: (Colour online). Scale invariance breaking in the BEC-BCS crossover.
(a) Values of $\tau^{2}(t)$ versus the expansion time $t$ at different
interactions. The expansion behaviors are represented for BEC
($1/k_{\textrm{F}}a=1.1$, red squares), unitary ($1/k_{\textrm{F}}a=0$, black
circles) and BCS ($1/k_{\textrm{F}}a=-0.68$, blue diamonds) regions,
respectively. The solid curves denote the calculations of Eq. (5) with
experimentally obtained $\gamma$. In the unitary regime,
$T/T_{\textrm{F}}=0.21(2)$. (b) Values of $\gamma$ in the BEC-BCS crossover.
The solid curve depicts the mean-field calculation at zero temperature as in
Ref. huCollectiveModesBallistic2004 . The dashed line denotes $\gamma=1$ valid
deeply in the BEC regime, and the dot-dashed line denotes $\gamma=2/3$ valid
at the unitary and deeply in the BCS regimes.

The obtained $\gamma$ in the BEC-BCS crossover is shown in Fig. 4(b). In the
unitary regime ($1/k_{\textrm{F}}a=0$), $\gamma\approx 2/3$, indicating the
scale-invariant expansion. On the BEC side ($1/k_{\textrm{F}}a>0$), $\gamma$
increases towards the molecule condensate limit with $\gamma=1$. On the BCS
side, $\gamma$ decreases to some extent. The experimental measurements have
the same variation trend with the mean-field calculation at zero temperature
huCollectiveModesBallistic2004 ; wangOscillatorylikeExpansionFermionic2020 .
The value of $\gamma$ in the unitary regime does not change with temperature.
But due to the finite-temperature effect (in the experiment,
$T/T_{\textrm{F}}=0.21(2)$), $\gamma$ is smaller than the zero-temperature
calculation on the BEC side, and on the BCS side it is larger. This can be
reasonably understood that, as temperature increases, $\gamma$ will change
towards the thermal gas with $\gamma=2/3$. On the BCS side, there is a shallow
dip in the zero-temperature calculation. In current condition, we could not
distinguish the non-monotonic behaviors on the BCS side due to the finite-
temperature effect and experimental fluctuation. The measurement of $\gamma$
is helpful to obtain the equation of state, frequencies of collective modes
and expansion behavior in the BEC-BCS crossover Stringari2008RMP ;
huCollectiveModesBallistic2004 ; Heiselberg2004PRL .

With obtained $\gamma$, we can determine the cloud radius,
${\langle}x^{2}{\rangle}_{0}=\gamma E_{\textrm{rel}}/m\omega^{2}$. Three
exemplary expansion behaviors in the BEC-BCS crossover are shown in Fig. 4(a),
which display the obvious deviation from that in the unitary regime
($1/k_{\textrm{F}}a=0$) when the interaction strength is tuned away from the
resonance. The expansion is fast on the BCS side ($1/k_{\textrm{F}}a=-0.68$)
and slow on the BEC side ($1/k_{\textrm{F}}a=1.1$), which can be well
calculated with Eq. (5).

In conclusion, we experimentally produce a spherical Fermi gas with variable
trapping frequencies and tunable interactions. We observe the unique feature
of the scale invariance induced by the coexistence of the spherical symmetry
and unitary interaction. The scale invariance is broken when the system is
tuned away from the resonant interaction. The virial theorem for the unitary
Femi gas has been verified. We also measure the effective exponent $\gamma$ in
the equation of state along the BEC-BCS crossover. The unitary Fermi gas in an
isotropic trap, a system with the hidden SO(2,1) symmetry, has several exact
properties, such as mapping between the trapped problem and the free-space
zero-energy problem, relations between the moments of the trapping potential
energy and the moments of the total energy, and undamped breathing mode with
twice the trapping frequency Werner2006 ; PitaevskiiPRA1997SOsymmetry . We
also have the opportunity to study the non-equilibrium dynamics in the
presence of conformal symmetry ZhoufeiPRA2019Conformal ;
ZhoufeiPRA2020Conformal ; ZhoufeiPRL2022Conformal .

We thank Shizhong Zhang and Xi-Wen Guan for carefully reading and revising the
paper, and Georgy Shlyapnikov for favorite discussions. This work has been
supported by the NKRDP (National Key Research and Development Program) under
Grant No. 2022YFA1404102, NSFC (Grant No. 12004398, 12121004, 11974384 and
12374250), CAS under Grant No. YJKYYQ20170025, and Hubei province under Grant
No. 2021CFA027.

Lu Wang, Xiangchuan Yan and Jing Min contributed equally to this work.

## References

  * (1) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga.  Feshbach resonances in ultracold gases.  Rev. Mod. Phys. 82, 1225 (2010).
  * (2) J. Kinast, A. Turlapov, J. E. Thomas, Q. Chen, J. Stajic, and K. Levin.  Heat capacity of a strongly interacting Fermi gas.  Science 307, 1296 (2005).
  * (3) S. Nascimbène, N. Navon, K. J. Jiang, F. Chevy, and C. Salomon.  Exploring the thermodynamics of a universal Fermi gas.  Nature 463, 1057 (2010).
  * (4) M. J. H. Ku, A. T. Sommer, L. W. Cheuk, and M. W. Zwierlein.  Revealing the superfluid lambda transition in the universal thermodynamics of a unitary Fermi gas.  Science 335, 563 (2012).
  * (5) L. A. Sidorenkov, M. K. Tey, R. Grimm, Y.-H. Hou, L. Pitaevskii, and S. Stringari.  Second sound and the superfluid fraction in a Fermi gas with resonant interactions.  Nature 498, 78 (2013).
  * (6) A. B. Bardon, S. Beattie, C. Luciuk, W. Cairncross, D. Fine, N. S. Cheng, G. J. A. Edge, E. Taylor, S. Zhang, S. Trotzky, and J. H. Thywissen.  Transverse demagnetization dynamics of a unitary Fermi gas.  Science 344, 722 (2014).
  * (7) P. B. Patel, Z. Yan, B. Mukherjee, R. J. Fletcher, J. Struck, and M. W. Zwierlein.  Universal sound diffusion in a strongly interacting Fermi gas.  Science 370, 1222 (2020).
  * (8) X. Li, X. Luo, S. Wang, K. Xie, X.-P. Liu, H. Hu, Y.-A. Chen, X.-C. Yao, and J.-W. Pan.  Second sound attenuation near quantum criticality.  Science 375, 528 (2022).
  * (9) C. Menotti, P. Pedri, and S. Stringari.  Expansion of an interacting Fermi gas.  Phys. Rev. Lett. 89, 250402 (2002).
  * (10) K. M. O’Hara, S. L. Hemmer, M. E. Gehm, S. R. Granade, and J. E. Thomas.  Observation of a strongly interacting degenerate Fermi gas of atoms.  Science 298, 2179 (2002).
  * (11) J. E. Thomas, J. Kinast, and A. Turlapov.  Virial theorem and universality in a unitary Fermi gas.  Phys. Rev. Lett. 95, 120402 (2005).
  * (12) E. Elliott, J. A. Joseph, and J. E. Thomas.  Observation of conformal symmetry breaking and scale invariance in expanding Fermi gases.  Phys. Rev. Lett. 112, 040405 (2014).
  * (13) S. Deng, Z.-Y. Shi, P. Diao, Q. Yu, H. Zhai, R. Qi, and H. Wu.  Observation of the efimovian expansion in scale-invariant Fermi gases.  Science 353, 371 (2016).
  * (14) C. Cao, E. Elliott, H. Wu, and J. E. Thomas.  Searching for perfect fluids: Quantum viscosity in a universal Fermi gas.  New J. Phys. 13, 075007 (2011).
  * (15) H. Guo, D. Wulin, C.-C. Chien, and K. Levin.  Microscopic approach to shear viscosities of unitary Fermi gases above and below the superfluid transition.  Phys. Rev. Lett. 107, 020403 (2011).
  * (16) M. Bluhm, J. Hou, and T. Schäfer.  Determination of the density and temperature dependence of the shear viscosity of a unitary Fermi gas based on hydrodynamic flow.  Phys. Rev. Lett. 119, 065302 (2017).
  * (17) C. Cao, E. Elliott, J. Joseph, H. Wu, J. Petricka, T. Schäfer, and J. E. Thomas.  Universal quantum viscosity in a unitary Fermi gas.  Science 331, 58 (2011).
  * (18) E. Elliott, J. A. Joseph, and J. E. Thomas.  Anomalous minimum in the shear viscosity of a Fermi gas.  Phys. Rev. Lett. 113, 020406 (2014).
  * (19) J. A. Joseph, E. Elliott, and J. E. Thomas.  Shear viscosity of a unitary Fermi gas near the superfluid phase transition.  Phys. Rev. Lett. 115, 020401 (2015).
  * (20) F. Werner and Y. Castin.  Unitary gas in an isotropic harmonic trap: Symmetry properties and applications.  Phys. Rev. A 74, 053604 (2006).
  * (21) See Supplemental Material for the formation of an isotropic trap, measurement of the trapping frequency, determination of the atomic temperature and in-situ atomic cloud size, fringe-removal analysis of the atomic images, and hydrodynamic description of the atomic expansion from an isotropic trap.
  * (22) T.-L. Ho.  Universal thermodynamics of degenerate quantum gases in the unitarity limit.  Phys. Rev. Lett. 92, 090402 (2004).
  * (23) D. T. Son.  Vanishing bulk viscosities and conformal invariance of the unitary Fermi gas.  Phys. Rev. Lett. 98, 020604 (2007).
  * (24) M. A. Escobedo, M. Mannarelli, and C. Manuel.  Bulk viscosities for cold Fermi superfluids close to the unitary limit.  Phys. Rev. A 79, 063623 (2009).
  * (25) K. Dusling and T. Schäfer.  Bulk viscosity and conformal symmetry breaking in the dilute Fermi gas near unitarity.  Phys. Rev. Lett. 111, 120603 (2013).
  * (26) X. Yan, D. Sun, L. Wang, J. Min, S. Peng, and K. Jiang.  Production of degenerate Fermi gases of 6Li atoms in an optical dipole trap.  Chin. Phys. Lett. 38, 056701 (2021).
  * (27) X. Yan, D. Sun, L. Wang, J. Min, S. Peng, and K. Jiang.  Observation of the BEC-BCS crossover in a degenerate Fermi gas of lithium atoms.  Chin. Phys. B 31, 016701 (2022).
  * (28) C. F. Ockeloen, A. F. Tauschinsky, R. J. C. Spreeuw, and S. Whitlock.  Detection of small atom numbers through image processing.  Phys. Rev. A 82, 061606 (2010).
  * (29) F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari.  Theory of Bose-Einstein condensation in trapped gases.  Rev. Mod. Phys. 71, 463 (1999).
  * (30) M. J. Holland, D. S. Jin, M. L. Chiofalo, and J. Cooper.  Emergence of interaction effects in Bose-Einstein condensation.  Phys. Rev. Lett. 78, 3801 (1997).
  * (31) R. Li, T. Gao, D. Zhang, S. Peng, L. Kong, X. Shen, and K. Jiang.  Expansion dynamics of a spherical Bose–Einstein condensate.  Chin. Phys. B 28, 106701 (2019).
  * (32) S. Giorgini, L. P. Pitaevskii, and S. Stringari.  Theory of ultracold atomic Fermi gases.  Rev. Mod. Phys. 80, 1215 (2008).
  * (33) H. Hu, A. Minguzzi, X.-J. Liu, and M. P. Tosi.  Collective modes and ballistic expansion of a Fermi gas in the BCS-BEC crossover.  Phys. Rev. Lett. 93, 190403 (2004).
  * (34) H. Heiselberg.  Collective modes of trapped gases at the BEC-BCS crossover.  Phys. Rev. Lett. 93, 040402 (2004).
  * (35) X. Wang, Y. Wu, X. Liu, Y. Wang, H. Chen, M. Maraj, Y. Deng, X.-C. Yao, Y.-A. Chen, and J.-W. Pan.  Oscillatory-like expansion of a fermionic superfluid.  Sci. Bull. 65, 7 (2020).
  * (36) L. P. Pitaevskii and A. Rosch.  Breathing modes and hidden symmetry of trapped atoms in two dimensions.  Phys. Rev. A 55, R853 (1997).
  * (37) J. Maki and F. Zhou.  Quantum many-body conformal dynamics: Symmetries, geometry, conformal tower states, and entropy production.  Phys. Rev. A 100, 023601 (2019).
  * (38) J. Maki and F. Zhou.  Far-away-from-equilibrium quantum-critical conformal dynamics: Reversibility, thermalization, and hydrodynamics.  Phys. Rev. A 102, 063319 (2020).
  * (39) J. Maki, S. Zhang, and F. Zhou.  Dynamics of strongly interacting Fermi gases with time-dependent interactions: Consequence of conformal symmetry.  Phys. Rev. Lett. 128, 040401 (2022).

## Supplemental materials

##  I Production of an isotropic optical trap

###  I.1 Experimental setup of the isotropic optical trap

Experimental setup of the isotropic trap is schematically displayed in Fig.
S1, which is also shown as Fig. 1(a) in the main text. Here we will display
more details about the setup. We prepare 6Li atomic degenerate Fermi gas with
two spin states $|F=1/2,m_{F}=\pm 1/2\rangle$. A pair of Helmholtz coils (not
shown in the figure) produce a homogeneous magnetic field in vertical
direction, which is used to tune the interactions. The Feshbach-resonance
magnetic field for 6Li atoms is about $834$ G. Some special techniques are
applied to produce the isotropic trap. Firstly, A pair of anti-Helmholtz coils
produce a gradient magnetic field in vertical direction to simultaneously
compensate the gravity force of the two spin states, which is valid for 6Li
atoms because the hyperfine interaction is much smaller than the Zeeman shift
at the strong magnetic field applied in the experiment. The magnetic gradient
is $B^{\prime}_{z}=1.05$ G/cm. Secondly, two elliptic optical beams with a
cross-sectional aspect ratio of $\sqrt{2}$, propagating perpendicularly to
each other, form the isotropic trap. Under this condition, the trapping
frequency can be varied by the optical power. The waist radius of the optical
beam is 60 $\mu$m and the wavelength is 1064 nm. Two laser beams with a
frequency difference of 76 MHz, propagating in vertical and horizonal
directions, respectively, are applied to detect two spin states separately.

Figure S1:  Experimental setup of the isotropic trap. Two elliptic optical
beams with a cross-sectional aspect ratio of $\sqrt{2}$, propagating
perpendicularly to each other, form the isotropic trap. A pair of anti-
Helmholtz coils produce a gradient magnetic field in vertical direction to
compensate the gravity. Another pair of Helmholtz coils (not shown) produce a
homogeneous magnetic field in vertical direction to tune the interactions. Two
laser beams with a frequency difference of 76 MHz, propagating in vertical and
horizonal directions, respectively, are applied to detect two spin states
separately.

###  I.2 Production of the elliptic beam with a cross-sectional aspect ratio
of $\sqrt{2}$

To form an isotropic optical trap, we should produce two elliptical beams
whose cross-sectional aspect ratio is $\sqrt{2}$. The optical configuration to
produce one beam is shown in Fig. S2. A 1064 nm laser beam outputs from a
polarization-maintaining fiber. The Glan-Taylor prism is used to purify the
optical polarization. A set of three cylindrical lens increases the radius in
$x$ direction, while maintains that in $z$ direction. The focus lengths and
distances of the lens are carefully selected to obtain the desired aspect
ratio of the cross-section, $w_{z}:w_{x}=1:\sqrt{2}$. After passing the final
achromatic doublets, the aspect ratio reverses, i.e.,
$w_{z}:w_{x}=\sqrt{2}:1$. In order to increase the optical stability, the
optical devices are constructed with stainless steel mounts, and the
experimental setup rests on an air-floating platform.

Figure S2: Optical configuration for preparing the elliptical beam. The
optical beam propagates in $y$ direction. The first row shows the optical
devices, the second row schematically displays the shape of the cross-section
on $x-z$ plane, and the third row indicates the value of the cross-sectional
aspect ratio $w_{z}:w_{x}$. $w_{x,z}$ is the waist radius in $x$ ($z$)
direction. After the output of the optical fiber, the cross-section has a
circular shape, i.e., $w_{z}:w_{x}=1:1$. A set of three cylindrical lens
increases the radius in $x$ direction, while maintains that in $z$ direction,
resulting an elliptical cross-section with $w_{z}:w_{x}=1:\sqrt{2}$. The focus
lengths of the lens are -50 mm, 50 mm, and 70 mm, respectively. The distances
between lens are 50 mm and 30 mm, respectively. The final achromatic doublets
with a focus length of 250 mm reverses the aspect ratio, i.e.,
$w_{z}:w_{x}=\sqrt{2}:1$.

###  I.3 Gravity compensation with a gradient magnetic field

To study interactions between two spin states in an isotropic optical trap
with variable trapping frequencies, the gravity force of the two spin states
should be simultaneously compensated. For 6Li atoms at the Feshbach resonant
magnetic field of $834$ G, the hyperfine interaction ($\approx 228$ MHz) is
much smaller than the Zeeman shift ($\approx 1.2$ GHz). In this condition, the
hyperfine quantum number F is no longer a good quantum number. The energy
shift due to the Zeeman effect at this strong magnetic field is

$\begin{split}\Delta
E=\mu_{\textrm{B}}(g_{J}m_{J}+g_{I}m_{I})B_{z},\end{split}$ (S1)

where $g$ is the $Land\acute{e}$ factor and $\mu_{\textrm{B}}$ is the Bohr
magneton. For the two lowest spin states $|F=1/2,m_{F}=\pm 1/2\rangle$, $L=0$
and $m_{J}=m_{S}=-1/2$. With $g_{I}\ll g_{S}$, the nuclear contribution can be
neglected. Then the energy shift is well approximated by

$\begin{split}\Delta E\simeq\mu_{\textrm{B}}g_{S}m_{S}B_{z},\end{split}$ (S2)

where $g_{S}=2$. The atoms in two spin states have the same magnetic moment
$\mu=m_{S}g_{S}\mu_{\textrm{B}}=-\mu_{\textrm{B}}$. To compensate the gravity
force of the two spin states, the magnetic field gradient $B^{\prime}_{z}$ can
be calculated,

$\begin{split}\mu_{\textrm{B}}B_{z}^{\prime}=mg.\end{split}$ (S3)

In the experiment, we use a pair of anti-Helmholtz coils in $z$ direction to
generate a quadrupole magnetic field, which is combined with the Feshbach
magnetic field to create a linear magnetic field in $z$ direction. According
to Eq. (S3), $B_{z}^{\prime}=1.05$ G/cm.

###  I.4 Residual Confinement of the Feshbach Magnetic Field

We probe the atomic expansion at a magnetic field of 834 G. So it is important
to analyze the residual confinement of the Feshbach magnetic field. A pair of
Helmholtz coils in $z$ direction are used to tune the Feshbach resonance. The
Feshbach magnetic field has a curvature which gives rise to an additional
trapping potential,

$\begin{split}U_{\textrm{mag}}=\frac{1}{2}\mu_{\textrm{B}}B_{z}^{\prime\prime}z^{2}.\end{split}$
(S4)

Then the trap frequency due to the magnetic field can be calculated,

$\begin{split}\omega_{\textrm{mag}}=\sqrt{\frac{\mu_{\textrm{B}}B_{z}^{\prime\prime}}{m}}.\end{split}$
(S5)

Through the parameters of the Feshbach coils, the curvature of the magnetic
field in the vertical direction is calculated to be $B^{\prime\prime}_{z}=0.1\
\textrm{G/cm}^{2}$ , while on the horizonal plane, it is negligibly small,
i.e., $B^{\prime\prime}_{\perp}\rightarrow 0\ \textrm{G/cm}^{2}$. Then the
residual trapping frequency originating from the Feshbach magnetic field is
only $\omega_{\textrm{mag}}=2\pi\times 1.54$ Hz, which is much smaller than
that of the optical dipole trap. According to Ref. ZhoufeiPRA2020Conformal ,
the scale factor $b(t)$ of the atomic cloud after released from the optical
trap in z-direction is

$\begin{split}b(t)=\left[\cos^{2}(\omega_{\textrm{mag}}t)+\frac{\omega^{2}_{\textrm{opt}}}{\omega^{2}_{\textrm{mag}}}\sin^{2}(\omega_{\textrm{mag}}t)\right]^{1/2}\end{split}$
(S6)

where $\omega_{\textrm{opt}}$ is the trapping frequency of the optical trap
and $t$ is the expansion time in the magnetic field. When
$\omega_{\textrm{mag}}=0$, the expansion evolves according to
$b_{0}(t)=(1+\omega^{2}_{\textrm{opt}}t^{2})^{1/2}$, which corresponds to the
scale invariant process. The effect of the residual magnetic confinement can
be defined by $\delta_{b}=(b_{0}(t)-b(t))/b_{0}(t)$, which represents the
shift of the atomic cloud size.

In the experiment, the trapping frequency of the optical trap is
$\omega_{\textrm{opt}}\approx 2\pi\times 1201$ Hz. For $t=2$ ms, the longest
expansion time in the experiment, $\delta_{b}\sim 6\times 10^{-5}$. So the
residual confinement effect of the Feshbach magnetic field can be ignored.

###  I.5 Theoretical analysis on how to form an isotropic optical dipole trap.

####  I.5.1 One Gaussian Beam

We first consider one focused Gaussian beam. Suppose that the optical beam
propagates along $z$ axis. Then the optical intensity is

$\begin{split}I(r,z)=\frac{I_{0}}{1+(z/z_{0})^{2}}\exp\left[\frac{-2r^{2}}{w_{0}^{2}}\right],\end{split}$
(S7)

where $I_{0}$ is the peak intensity, $z_{0}=\pi w_{0}^{2}/\lambda$ is the
Rayleigh length, and $w_{0}$ is waist radius (the $1/e^{2}$ intensity radius
of the beam at the focus). Then the trapping potential is given by

$\begin{split}U(r,z)=-\frac{U_{0}}{1+(z/z_{0})^{2}}\exp\left[\frac{-2r^{2}}{w_{0}^{2}}\right],\end{split}$
(S8)

where

$\begin{split}&U_{0}=\frac{3\pi
c^{2}\Gamma}{2\omega_{0}^{3}\Delta\omega}I_{0},\\\ &I_{0}=\frac{2P}{\pi
w_{0}^{2}}.\end{split}$ (S9)

Here $c$ is the speed of light, $\omega_{0}=2\pi c/\lambda_{0}$, $\omega=2\pi
c/\lambda$, and $\Delta\omega=\omega_{0}-\omega$. In the experiment, the
natural line width is $\Gamma$=$2\pi\times$ 5.87 MHz for 6Li atom,
$\lambda_{0}=671$ nm and $\lambda=1064$ nm. To determine the trapping
frequency of the trap, we expand the trapping potential of Eq. (S8) into
Taylor series around the center $(x,y,z)=(0,0,0)$,

$\begin{split}U(r,z)\simeq-
U_{0}(1-\frac{z^{2}}{z_{0}^{2}})(1-\frac{2r^{2}}{w_{0}^{2}})+\ldots\\\ \simeq-
U_{0}+\frac{U_{0}}{z_{0}^{2}}z^{2}+2\frac{U_{0}}{w_{0}^{2}}r^{2}+\ldots.\end{split}$
(S10)

Then

$\begin{split}&\frac{U_{0}}{z_{0}^{2}}z^{2}=\frac{1}{2}m\omega_{\textrm{axial}}^{2}z^{2},\\\
&2\frac{U_{0}}{w_{0}^{2}}r^{2}=\frac{1}{2}m\omega_{\textrm{radial}}^{2}r^{2},\end{split}$
(S11)

where $\omega_{\textrm{axial}}$ and $\omega_{\textrm{radial}}$ are the
trapping frequencies of the optical beam in axial and radial directions,
respectively.

$\begin{split}&\omega_{\textrm{axial}}=\sqrt{\frac{2U_{0}}{mz_{o}^{2}}},\\\
&\omega_{\textrm{radial}}=\sqrt{\frac{4U_{0}}{mw_{o}^{2}}}.\end{split}$ (S12)

In the experiment, the waist radius $w_{o}$ is 60 $\mu$m. Then
$\omega_{\textrm{axial}}/\omega_{\textrm{radial}}\simeq 4.0\times 10^{-3}$.
The trapping effect in axial direction of the optical beam is negligibly
small.

####  I.5.2 Two Orthogonal Beams with the Circular Cross-Section

We consider two identical Gaussian beams with the circular cross-section. The
two optical beams propagate along $x$ and $y$ axes, respectively. Then the
trapping potential can be expressed as

$\begin{split}U_{\textrm{OPT}}=U_{\textrm{OPT1}}+U_{\textrm{OPT2}}=-\frac{U_{0}}{1+(x/x_{0})^{2}}\exp\left[-\frac{2(y^{2}+z^{2})}{w_{0}^{2}}\right]-\frac{U_{0}}{1+(y/y_{0})^{2}}\exp\left[-\frac{2(x^{2}+z^{2})}{w_{0}^{2}}\right],\end{split}$
(S13)

where $x_{0}=y_{0}$ is the Rayleigh length and $w_{0}$ is the waist radius.

We expansion Eq. (S13) into Tailor series at the point $(x,y,z)=(0,0,0)$,

$\begin{split}U_{OPT}\simeq-2U_{0}+\left(\frac{U_{0}}{x_{0}^{2}}+2\frac{U_{0}}{w_{0}^{2}}\right)x^{2}+\left(\frac{U_{0}}{y_{0}^{2}}+2\frac{U_{0}}{w_{0}^{2}}\right)y^{2}+4\frac{U_{0}}{w_{0}^{2}}z^{2}+\ldots.\end{split}$
(S14)

As mentioned above, the trapping effect in the axial direction of the optical
beam can be ignored, and we only consider the trapping effect in the radial
direction. Then the trapping frequencies can be calculated as

$\begin{split}\omega_{x}=\omega_{y}\simeq\omega_{\textrm{radial}},\\\
\omega_{z}\simeq\sqrt{2}\omega_{\textrm{radial}}.\end{split}$ (S15)

According to Eq. (S15), the trapping frequency in the vertical direction is
about $\sqrt{2}$ times that in the horizontal direction.

####  I.5.3 Two Orthogonal Beams with the Elliptical Cross-Section

According to the above analysis, we couldn’t form an isotropic trap using the
beams with the circular cross-section. Here we consider two beams with the
elliptical cross-section. The two optical beams propagate along $x$ and $y$
axes, respectively. The waist radius in $x$ and $y$ direction is different
from that in $z$ direction, i.e., $w_{x0}=w_{y0}=w_{0}$, $w_{x0}\neq w_{z0}$.
Ignoring the trapping effect in the axial direction of the optical beam as
mentioned above, then the trapping potential is calculated as

$\begin{split}U_{\textrm{OPT}}=&U_{\textrm{OPT1}}+U_{\textrm{OPT2}}\approx-
U_{0}\exp\left[-2\left(\frac{y^{2}}{w_{0}^{2}}+\frac{z^{2}}{w_{z0}^{2}}\right)\right]-U_{0}\exp\left[-2\left(\frac{x^{2}}{w_{0}^{2}}+\frac{z^{2}}{w_{z0}^{2}}\right)\right].\end{split}$
(S16)

Ignoring the trapping effect in the axial direction of the optical beam, we
expand Eq. (S16) into Tailor series at the point $(x,y,z)=(0,0,0)$,

$\begin{split}U_{\textrm{OPT}}\simeq-2U_{0}+2\frac{U_{0}}{w_{0}^{2}}x^{2}+2\frac{U_{0}}{w_{0}^{2}}y^{2}+4\frac{U_{0}}{w_{z0}^{2}}z^{2}+\ldots.\end{split}$
(S17)

Then the trapping frequencies can be calculated,

$\begin{split}\omega_{x}=\omega_{y}\simeq\omega_{\textrm{radial}},\\\
\frac{\omega_{z}}{\omega_{x}}=\frac{\sqrt{2}w_{0}}{w_{z0}}.\end{split}$ (S18)

According to Eq. (S18), if $w_{z0}=\sqrt{2}w_{0}$, for any optical power, the
trapping frequencies along three orthogonal axes are the same,

$\begin{split}\omega_{x}=\omega_{y}=\omega_{z}.\end{split}$ (S19)

In conclusion, in our experiment, the gravity force is compensated with a
gradient magnetic field, and the residual confinement of the Feshbach magnetic
field is negligibly small. Then we can form an isotropic optical trap using
two orthogonal beams with the elliptical cross-section. The aspect ratio of
the cross-section is $\sqrt{2}$. The trapping frequency can be varied by the
power of the optical beam.

##  II Measurement of the trapping frequency

In the experiment, the trapping frequency of the trap is determined by
measuring the center-of-mass oscillation of the atomic cloud. We perturb the
position of the atomic cloud in the spherical trap by controlling a pulse of
the elongated optical trap. We tune the relative position of the elongated
trap to the spherical trap, simultaneously shifting positions of the atomic
cloud in three axes. After switching off the pulse of the elongated trap, the
atoms will oscillate in the trap. After different waiting time in the trap, we
probe atoms with a time-of-flight (TOF) of 1 ms. The atomic temperature is
above the superfluid temperature $T_{\textrm{c}}$, and the atomic position
$x_{i}$ ($i\rightarrow x,y,z$) is obtained by fitting the density profile
using a Gaussian distribution. We can determine the oscillation of the center-
of-mass $\Delta x_{i}(t)=x_{i}(t)-\bar{x}_{i}$, where $t$ is the waiting time
in the trap, and $\bar{x}_{i}$ is the mean value of $x_{i}(t)$. Then we can
fit $\Delta x_{i}(t)$ using a sinusoidal function,

$\begin{split}\Delta x_{i}(t)=A_{i}\sin(\omega_{i}t+\phi_{i}),\end{split}$
(S20)

where $A_{i}$ is the oscillation amplitude, and $\omega_{i}$ is the trapping
frequency in $x_{i}$ direction.

Fig. S3 shows the measurement results, where the optical power is 3.8 W. The
trapping frequencies are
$(\omega_{x},\omega_{y},\omega_{z})=2\pi\times(1234(6),1165(11),1204(3))$ Hz,
which are almost the same along three axes.

Figure S3: Center-of-mass oscillations of the atomic cloud. $\Delta x_{i}$
($i\rightarrow x,y,z$) denotes the distance away from the mean value of the
atomic position. The error bar is the standard deviation of three
measurements. The black solid line denotes the sinusoidal fitting.

##  III Measurement of the temperature and in-situ atomic cloud size

###  III.1 How to obtain the temperature of the unitary Fermi gas

The temperature of the unitary Fermi gas is obtained using the method similar
to J. E. Thomas Kinast2005  and C. J. Vale’s groups
veeravalliBraggSpectroscopyStrongly2008a . The normalized one-dimensional
profile of a non-interacting Fermi gas is

$n(r_{i}(t),T)=-\frac{3N}{\sqrt{\pi}\sigma_{\textrm{F}i}(t)}\left(\frac{T}{T_{\textrm{F}}}\right)^{5/2}\textrm{Li}_{5/2}\left[-\exp\left(\frac{\frac{\mu}{E_{\textrm{F}}}-\frac{r_{i}^{2}}{\sigma_{\textrm{F}i}^{2}(t)}}{T/T_{\textrm{F}}}\right)\right],$
(S21)

where $N$ is the total atom number, $\textrm{Li}_{5/2}$ is the polylogarithm
function, $\mu$ is the chemical potential and $\sigma_{\textrm{F}i}(t)$ is the
Thomas-Fermi radius of the atomic cloud after released from the trap.

The fitting procedure for the unitary Fermi gas is similar to that of the non-
interacting gas as mentioned above, except that the Fermi radius
$\sigma_{\textrm{F}i}$ and Fermi temperature $T_{\textrm{F}}$ should be
replaced by $\sigma_{\textrm{F}i}^{*}$ and $T_{\textrm{F}}^{*}$ ,
respectively. At unitarity,
$\sigma_{\textrm{F}i}^{*}=(1+\beta)^{1/4}\sigma_{\textrm{F}i}$ and
$T_{\textrm{F}}^{*}=(1+\beta)^{1/2}T_{\textrm{F}}$ , where $\beta$ is a
universal constant, and the fitting profile becomes

$n(r_{i}(t),T)=-\frac{3N}{\sqrt{\pi}\sigma_{\textrm{F}i}^{*}(t)}\left(\tilde{T}\right)^{5/2}\textrm{Li}_{5/2}\left[-exp\left(q-\frac{r_{i}^{2}(t)}{(\sigma_{\textrm{F}i}^{*}(t))^{2}\tilde{T}}\right)\right],$
(S22)

where $q={\mu}/(E_{\textrm{F}}\tilde{T})$ ,
$E_{\textrm{F}}=\hbar\omega_{0}(3N)^{1/3}$ and $\tilde{T}$ is the empirical
temperature given by

$\tilde{T}=\frac{T}{T_{\textrm{F}}\sqrt{1+\beta}}.$ (S23)

The expansion behavior is known in the unitary regime.
$\sigma_{\textrm{F}i}^{*}(t)$ can be calculated by the atom number and
trapping frequency. Another way to obtain $\sigma_{\textrm{F}i}^{*}(t)$ is to
fit the atomic density profile using the zero-temperature Thomas-Fermi
distribution at the lowest temperature. We fix $\sigma_{\textrm{F}i}^{*}(t)$
constant for the fits at all higher temperatures, leaving only $\tilde{T}$ and
$q$ as the free parameters. Then we obtain the value of $T/T_{\textrm{F}}$
from Eq. (S23), where $\beta=-0.56$ luo2009thermodynamic .

###  III.2 How to obtain the in-situ atomic cloud size of the unitary Fermi
gas in finite temperature

At finite temperature, for a non-interacting Fermi gas, the in-situ mean
square size $\left\langle r_{i}^{2}(T)\right\rangle_{0}$ is given by

$\left\langle
r_{i}^{2}(T)\right\rangle_{0}=\frac{\sigma_{\textrm{F}i}^{2}}{8}\frac{E}{E_{0}}\left(\frac{T}{T_{F}}\right)$
(S24)

where $\sigma_{\textrm{F}i}=\sqrt{{2E_{F}}/{m\omega_{0}^{2}}}$ is the Fermi
radius of the atomic cloud in the trap, $E_{0}$ is the ground energy.

For the unitary Fermi gas,

$\left\langle
r_{i}^{2}(\tilde{T})\right\rangle_{0}=\frac{(\sigma_{\textrm{F}i}^{*})^{2}}{8}\frac{E}{E_{0}}(\tilde{T}).$
(S25)

As seen in Ref. Kinast2005 , with the same value of $T/T_{F}$ and $\tilde{T}$
for the non-interacting Fermi gas and unitary Fermi gas,
$\frac{E}{E_{0}}(\frac{T}{T_{F}})=\frac{E}{E_{0}}(\tilde{T)}$. After knowing
the temperature $\tilde{T}$ of the unitary Fermi gas using the method
mentioned above, the value of $\frac{E}{E_{0}}(\tilde{T)}$ can be obtained.
According to Eq. (S25) and the relation
$\sigma_{\textrm{F}i}^{*}=(1+\beta)^{1/4}\sigma_{\textrm{F}i}$ , we can obtain
the in-situ size $\left\langle r_{i}^{2}\right\rangle_{0}$ of the unitary
Fermi gas in finite temperature.

##  IV Optimization of the Atomic image using the fringe-removal analysis

To obtain the density distribution of the gas, three atomic images should be
taken, respectively, as $P_{\textrm{abs}}$, $P_{\textrm{ref}}$ and
$P_{\textrm{bg}}$. $P_{\textrm{abs}}$ is the absorption image with atom and
light, $P_{\textrm{ref}}$ is the reference image without atom, and
$P_{\textrm{bg}}$ is the background image without light and atom. The optical
density (OD) distribution is obtained from

$P_{\textrm{od}}=\ln\left(\frac{P_{\textrm{ref}}-P_{\textrm{bg}}}{P_{\textrm{abs}}-P_{\textrm{bg}}}\right).$
(S26)

Due to changes in the intensity and spatial position of the imaging light,
$P_{\textrm{ref}}$ is different from the background of $P_{\textrm{abs}}$,
leading to fringes and other noises in $P_{\textrm{od}}$ (see Fig. S4(a)). As
known in Ref. Ockeloen2010 , these noises can be reduced by using an algorithm
to synthesize a new reference image $P_{\textrm{refn}}$, which is closest to
the background of $P_{\textrm{abs}}$. Replacing $P_{\textrm{ref}}$ with
$P_{\textrm{refn}}$ in Eq. (S26), we can optimize the atomic image
$P_{\textrm{od}}$.

A set of reference images compose a background library $R$, whose linear
superposition gives $P_{\textrm{refn}}$ by

$P_{\textrm{refn}}=RC.$ (S27)

The coefficient matrix $C$ is determined by setting the least square
difference between $P_{\textrm{refn}}$ and $P_{\textrm{abs}}$ in the regions
without absorption of atoms. Setting the partial derivative of the square
difference with respect to the coefficient to zero, a set of equations are
obtained, where the solutions give the coefficients.

Fig. S4 displays an example of the image optimization. Without optimization,
there are many fringes and noises in the background, as shown in Fig. S4(a).
After optimization, the fringes and noises are greatly reduced in Fig. S4(b).
In Fig. S4(c), we use a Gaussian function to fit the one dimensional OD,

$\textrm{OD}(x)=\textrm{OD}_{0}+A\exp\left(\frac{-x^{2}}{2\sigma^{2}_{x}}\right).$
(S28)

Without optimization, $\textrm{OD}_{0}=0.63\pm 0.13$, $\sigma_{x}=27.26\pm
0.87$, $\textrm{A}=10.73\pm 0.27$, and $\chi^{2}$ of the fitting is 0.45.
After optimization, $\textrm{OD}_{0}=0.20\pm 0.09$, $\sigma_{x}=25.95\pm
0.62$, $\textrm{A}=10.93\pm 0.20$, and $\chi^{2}$ of the fitting is 0.26. It
can be seen that, through image optimization, the background level and fitting
uncertainty are well reduced. The atomic cloud radius also changes, which
should be more accurate.

Figure S4: Optimization of atomic images. (a) is for the atomic image without
optimization and (b) with optimization. The integrated OD in $x$ direction is
shown in (c). The solid line is the fitting using a Gaussian function. The OD
distributions with and without optimization are displayed together for
comparison.

##  V Hydrodynamic description of the atomic expansion from an isotropic trap

A Gaussian distribution is used to fit the atomic density profile,
$n(x)=\textrm{A}e^{(-x^{2}/2\sigma_{x}^{2})}$, and the fitted value of
$\sigma_{x}$ is related to the mean square cloud size in $x$ direction by

${\langle}x^{2}{\rangle}=\frac{1}{N}\int{nx^{2}\text{d}x}=\sigma_{x}^{2}.$
(S29)

We consider how ${\langle}x_{i}^{2}{\rangle}$ ($i=x,y,z$) evolves with
expansion time $t$,

$\displaystyle\frac{\text{d}{\langle}x_{i}^{2}{\rangle}}{\text{d}t}$
$\displaystyle=\frac{1}{N}\int{\frac{\partial n}{\partial
t}x_{i}^{2}\text{d}x}.$ (S30)

Using the continuity equation of the hydrodynamic description for one
component fluid Cao2011 ; Cao2011a ; houScalingSolutionsTwofluid2013 ,
${\partial n}/{\partial t}+\nabla\cdot(n\mathbf{v})=0$, Eq. (S30) can be
written as

$\displaystyle\frac{\text{d}{\langle}x_{i}^{2}{\rangle}}{\text{d}t}$
$\displaystyle=\frac{1}{N}\int{[-\nabla\cdot(n\mathbf{v})]x_{i}^{2}\text{d}x}$
(S31)
$\displaystyle=\frac{1}{N}\int{[-\nabla\cdot(x_{i}^{2})]n\mathbf{v}\text{d}x}+\frac{1}{N}\int{[-\nabla\cdot(x_{i}^{2}n\mathbf{v})]\text{d}x}$
$\displaystyle=\frac{1}{N}\int{[-\nabla\cdot(x_{i}^{2})]n\mathbf{v}\text{d}x}$
$\displaystyle=\frac{1}{N}\int{2x_{i}n{v_{i}}\text{d}x}$
$\displaystyle=2{\langle}x_{i}{v_{i}}{\rangle}.$

Using the same procedure, we can write the evolution of
${\langle}x_{i}{v_{i}}{\rangle}$ with time:

$\displaystyle\frac{\text{d}{\langle}x_{i}v_{i}{\rangle}}{\text{d}t}$
$\displaystyle=\frac{1}{N}\int{\frac{\partial n}{\partial
t}x_{i}v_{i}\text{d}x}+\frac{1}{N}\int{nx_{i}\frac{\partial v_{i}}{\partial
t}\text{d}x}$ (S32)
$\displaystyle=\frac{1}{N}\int{[\nabla\cdot(x_{i}v_{i})]n\mathbf{v}\text{d}x}+\frac{1}{N}\int{nx_{i}\frac{\partial
v_{i}}{\partial t}\text{d}x}$
$\displaystyle={\langle}x(\partial_{t}+\mathbf{v}\cdot\nabla)v_{i}{\rangle}+{\langle}v_{i}^{2}{\rangle}.$

Combining Eq. (S31) and (S32), we obtain

$\displaystyle\frac{\text{d}^{2}}{\text{dt}^{2}}\frac{{\langle}x_{i}^{2}{\rangle}}{2}$
$\displaystyle={\langle}x_{i}(\partial_{t}+\mathbf{v}\cdot\nabla)v_{i}{\rangle}+{\langle}v_{i}^{2}{\rangle}.$
(S33)

The first term is similar to Euler’s equation,

$mn(\partial_{t}+\mathbf{v}\cdot\nabla)v=-\partial_{i}P-n\partial_{i}U+\sum_{j}\partial_{j}(\eta\sigma_{ij}+\zeta\sigma^{\prime}\delta_{ij}),$
(S34)

where $m$ is the atomic mass, $P$ is the scalar pressure, $U$ is the trapping
potential energy, and the last term on the right side denotes the friction
forces due to shear $\eta$ and bulk $\zeta$ viscosity. The viscosity can be
written as $\eta\equiv\alpha_{S}\hbar n$ and $\zeta\equiv\alpha_{B}\hbar n$.
$\sigma_{ij}\equiv{\partial v_{i}}/{\partial x_{i}}+{\partial v_{j}}/{\partial
x_{j}}-2/3\delta_{ij}\nabla\cdot\mathbf{v}$,
$\sigma^{\prime}=\nabla\cdot\mathbf{v}$. We then take the density averaged
product of the Euler’s equation with a position component, Eq. (S34) can be
written as

$\displaystyle\frac{1}{N}\int
nx_{i}(\partial_{t}+\mathbf{v}\cdot\nabla)v_{i}\text{d}^{3}r$
$\displaystyle=-\frac{1}{Nm}\int
x_{i}\partial_{i}P\text{d}^{3}r-\frac{1}{Nm}\int
x_{i}n\partial_{i}U\text{d}^{3}r+\frac{1}{Nm}\sum_{j}\int
x_{i}\partial_{j}(\eta\sigma_{ij}+\zeta\sigma^{\prime}\delta_{ij})\text{d}^{3}r$
(S35) $\displaystyle=\frac{1}{Nm}\int P\text{d}^{3}r-\frac{1}{Nm}\int
x_{i}n\partial_{i}U\text{d}^{3}r-\frac{1}{Nm}\int(\eta\sigma_{ii}+\zeta\sigma^{\prime})\text{d}^{3}r.$

The pressure, trapping potential and viscosity terms must be zero for
$x_{i}\rightarrow\pm\infty$. Then combining Eq. (S33) and (S35), we obtain

$\frac{\text{d}^{2}}{\text{dt}^{2}}\frac{{\langle}x_{i}^{2}{\rangle}}{2}=\frac{1}{Nm}\int
P\text{d}^{3}r-\frac{1}{m}{\langle}x_{i}\partial_{i}U{\rangle}-\frac{\hbar}{m}{\langle}\alpha_{S}\sigma_{ii}+\alpha_{B}\sigma^{\prime}{\rangle}+{\langle}v_{i}^{2}{\rangle}.$
(S36)

Eq. (S36) determines the evolution of the mean square cloud radius, which
depends on the conservative forces from the scale pressure and the trapping
potential, as well as the dissipative forces arising from the shear and bulk
viscosity. For a spherical system, expansion behaviors in all directions are
the same, i.e., ${\partial v_{i}}/{\partial x_{i}}={\partial v_{j}}/{\partial
x_{j}}={\partial v_{k}}/{\partial x_{k}}$. Then

$\displaystyle\sigma_{ii}$ $\displaystyle\equiv\frac{\partial v_{i}}{\partial
x_{i}}+\frac{\partial v_{i}}{\partial
x_{i}}-\frac{2}{3}\delta_{ii}\nabla\cdot\mathbf{v}$ (S37)
$\displaystyle=\frac{\partial v_{i}}{\partial x_{i}}+\frac{\partial
v_{i}}{\partial x_{i}}-\frac{2}{3}(\frac{\partial v_{i}}{\partial
x_{i}}+\frac{\partial v_{j}}{\partial x_{j}}+\frac{\partial v_{k}}{\partial
x_{k}})$ $\displaystyle=0,$ $\displaystyle\sigma^{\prime}$
$\displaystyle=\nabla\cdot\mathbf{v}=\frac{\partial v_{i}}{\partial
x_{i}}+\frac{\partial v_{j}}{\partial x_{j}}+\frac{\partial v_{k}}{\partial
x_{k}}=3\frac{\partial v_{i}}{\partial x_{i}}.$

For simplicity, we only consider the atomic expansion in $x$ direction.
Inserting Eq. (S37) into Eq. (S36), we can find that for a spherical system,
the evolution of the mean square cloud radius can be written as

$\frac{\text{d}^{2}}{\text{dt}^{2}}\frac{{\langle}x^{2}{\rangle}}{2}=\frac{1}{Nm}\int
P\text{d}^{3}r-\frac{1}{m}{\langle}x\partial_{x}U{\rangle}-\frac{\hbar}{m}{\langle}3\alpha_{B}\partial_{x}{v}_{x}{\rangle}+{\langle}v_{x}^{2}{\rangle}.$
(S38)

Now we need to eliminate ${\langle}v_{x}^{2}{\rangle}$ by using the energy
conservation equation:

$\frac{\text{d}}{\text{d}t}\int\text{d}^{3}r(n\frac{1}{2}m\mathbf{v}^{2}+\varepsilon+nU)=0.$
(S39)

At $t\geqslant 0^{+}$, atoms are released from the optical trap ($U=0$),

$\displaystyle t=0^{+}\qquad E$
$\displaystyle=\frac{1}{N}\int\text{d}^{3}r\varepsilon_{0},$ (S40)
$\displaystyle t>0^{+}\qquad E$
$\displaystyle=\frac{1}{N}\int\text{d}^{3}r\varepsilon+\frac{m}{2}{\langle}\mathbf{v}^{2}{\rangle}.$

Taking $\triangle P\equiv P-{2}/{3}\varepsilon$ into Eq. (S40),

$\displaystyle\frac{1}{N}\int\text{d}^{3}r\varepsilon_{0}$
$\displaystyle=\frac{3}{2N}\int\text{d}^{3}rP_{0}-\frac{3}{2N}\int\text{d}^{3}r\triangle
P_{0}$ (S41)
$\displaystyle=\frac{1}{N}\int\text{d}^{3}r\varepsilon+\frac{m}{2}{\langle}\mathbf{v}^{2}{\rangle}.$

Then we find

$\displaystyle\langle v^{2}\rangle$
$\displaystyle=\frac{2}{m}(\frac{3}{2N}\int\text{d}^{3}rP_{0}-\frac{3}{2N}\int\text{d}^{3}r\triangle
P_{0}-\frac{1}{N}\int\text{d}^{3}r\varepsilon)$ (S42)
$\displaystyle=\frac{3}{m}\langle r\cdot\nabla
U\rangle_{0}-\frac{3}{mN}\int\text{d}^{3}r\triangle
P_{0}-\frac{2}{mN}\int\text{d}^{3}r\varepsilon.$

For a spherical system, $\langle v_{x}^{2}\rangle=\langle v^{2}\rangle/3$.

Before release from the trap at $t=0^{-}$, $\mathbf{v}=0$. Eq. (S38) can be
written as

$\frac{1}{N}\int\text{d}^{3}rP_{0}=\langle{x}\partial_{x}U\rangle_{0},$ (S43)

where the subscript $()_{0}$ describe the initial condition in the trap.
Combining Eqs. (S42), (S43) and (S38), we obtain

$\displaystyle\frac{\text{d}^{2}}{\text{d}^{2}t}\frac{m{\langle}x^{2}{\rangle}}{2}$
$\displaystyle=\langle
x\partial_{x}U\rangle_{0}-\frac{1}{N}\int\text{d}^{3}r\triangle
P_{0}-\frac{2}{3N}\int\text{d}^{3}r\varepsilon$ (S44)
$\displaystyle+\frac{1}{N}\int\text{d}^{3}r\triangle
P+\frac{2}{3N}\int\text{d}^{3}r\varepsilon-\hbar{\langle}3\alpha_{B}\partial_{x}{v_{x}}{\rangle}$
$\displaystyle=\langle
x\partial_{x}U\rangle_{0}+\frac{1}{N}\int\text{d}^{3}r(\triangle P-\triangle
P_{0})-\hbar{\langle}3\alpha_{B}\partial_{x}{v_{x}}{\rangle}.$

This is the atomic expansion evolution from an isotropic harmonic trap. For a
unitary Fermi gas, $\triangle P=0$ and $\alpha_{B}=0$. Eq. (S44) can be
written as

$\begin{split}\frac{d^{2}}{dt^{2}}\frac{m{\langle}x^{2}{\rangle}}{2}={\langle}x{\frac{\partial
U}{\partial x}}{\rangle}_{0},\end{split}$ (S45)

${\langle}x^{2}{\rangle}={\langle}x^{2}{\rangle}_{0}+\omega_{x}^{2}{\langle}x^{2}{\rangle}_{0}t^{2}.$
(S46)

Eqs. (S45) and (S46) display the scale-invariant expansion of a unitary Fermi
gas from an isotropic trap, which is similar to a non-interacting Fermi gas.
While away from the resonant interaction, $\triangle P\neq 0$ and
$\alpha_{B}\neq 0$, the scale-invariant expansion will be broken.

## References

  * (1) J. Maki and F. Zhou.  Far-away-from-equilibrium quantum-critical conformal dynamics: Reversibility, thermalization, and hydrodynamics.  Phys. Rev. A 102, 063319 (2020).
  * (2) J. Kinast, J. E. Turlapov, A.and Thomas, Q. Chen, J. Stajic, and K. Levin.  Heat capacity of a strongly interacting Fermi gas.  Science 307, 1296 (2005).
  * (3) G. Veeravalli, E. Kuhnle, P. Dyke, and C. J. Vale.  Bragg spectroscopy of a strongly interacting Fermi gas.  Phys. Rev. Lett. 101, 250403 (2008).
  * (4) L. Luo and J. E. Thomas.  Thermodynamic measurements in a strongly interacting fermi gas.  J. Low Temp. Phys. 154, 1 (2009).
  * (5) C. F. Ockeloen, A. F. Tauschinsky, R. J. C. Spreeuw, and S. Whitlock.  Detection of small atom numbers through image processing.  Phys. Rev. A 82, 061606 (2010).
  * (6) C. Cao, E. Elliott, J. Joseph, H. Wu, J. Petricka, T. Schäfer, and J. E. Thomas.  Universal quantum viscosity in a unitary Fermi gas.  Science 331, 58 (2011).
  * (7) C. Cao, E Elliott, H. Wu, and J. E. Thomas.  Searching for perfect fluids: Quantum viscosity in a universal Fermi gas.  New J. Phys. 13, 075007 (2011).
  * (8) Y. Hou, L. P. Pitaevskii, and S. Stringari.  Scaling solutions of the two-fluid hydrodynamic equations in a harmonically trapped gas at unitarity.  Phys. Rev. A 87, 033620 (2013).


# Letter of Intent: Search for sub-millicharged particles at J-PARC

Suyong Choi  Korea University, Seoul, Korea  Jeong Hwa Kim  Korea University,
Seoul, Korea  Eunil Won  Korea University, Seoul, Korea  Jae Hyeok Yoo  Korea
University, Seoul, Korea  Matthew Citron  University of California, Santa
Barbara, California, USA  David Stuart  University of California, Santa
Barbara, California, USA  Christopher S. Hill  The Ohio State University,
Columbus, Ohio, USA  Andy Haas  New York University, New York, New York, USA
Jihad Sahili  Lebanese University, Hadeth-Beirut, Lebanon  Haitham Zaraket
Lebanese University, Hadeth-Beirut, Lebanon  A. De Roeck  CERN, Geneva,
Switzerland  Martin Gastal  CERN, Geneva, Switzerland

###### Abstract

We propose a new experiment sensitive to the detection of millicharged
particles produced at the $30$ GeV proton fixed-target collisions at J-PARC.
The potential site for the experiment is B2 of the Neutrino Monitor building,
$280$ m away from the target. With $\textrm{N}_{\textrm{POT}}=10^{22}$, the
experiment can provide sensitivity to particles with electric charge $3\times
10^{-4}\,e$ for mass less than $0.2$ $\textrm{GeV}/\textrm{c}^{2}$ and
$1.5\times 10^{-3}\,e$ for mass less than $1.6$ $\textrm{GeV}/\textrm{c}^{2}$.
This brings a substantial extension to the current constraints on the charge
and the mass of such particles.

##  1 Motivation for the experiment

Electric charge quantization is a long-standing question in particle physics.
While fractionally charged particles (millicharged particles hereafter) have
typically been thought to preclude the possibility for Grand Unified Theories
(GUTs), well-motivated dark-sector models have been proposed to predict the
existence of millicharged particles while preserving the possibility for
unification. Such models can contain a rich internal structure, providing
candidate particles for dark matter. One well motivated mechanism that leads
to millicharged particles is to introduce a new $U(1)$ in dark sector with a
massless dark-photon ($A^{\prime}$) and a massive dark-fermion ($\chi$) [1],

$\displaystyle\mathcal{L}_{\textrm{dark
sector}}=-\frac{1}{4}A^{\prime}_{\mu\nu}A^{\prime\mu\nu}+i\bar{\chi}\left(\not{\partial}+ie^{\prime}\not{A}^{\prime}+i\textrm{m}_{\chi}\right)\chi-\frac{\kappa}{2}A^{\prime}_{\mu\nu}B^{\mu\nu}$
(1)

where the last term shows that $A^{\prime}$ and $B$ kinetically mix with the
mixing parameter $\kappa$. After replacing $\not{A}^{\prime}$ with
$\not{A}^{\prime}+\kappa B_{\mu}$, the coupling between the dark fermion and
$B$ becomes $\kappa e^{\prime}$ ($~{}\kappa e^{\prime}\bar{\chi}\not{B}\chi$).
The charge of $\chi$ varies by the size of mixing, so $\chi$ can be a
millicharged particle. Hereafter, $\chi$ is used to denote millicharged
particles.

A number of experiments have searched for millicharged particles, including in
an electron fixed-target experiment [2], proton-proton colliders [3, 4, 5],
and neutrino experiments [6]. A comprehensive review is in Reference [7]. In
the parameter space of the charge ($Q$) and the mass ($m_{\chi}$), the region
of $m_{\chi}>0.1$ $\textrm{GeV}/\textrm{c}^{2}$ and $Q<10^{-2}e$ is largely
unexplored.

Proton fixed-target experiments provide a solid testing ground for $\chi$s.
The particle flux is much larger than the collider experiments and they can
reach higher energy regime than electron fixed-target experiments. The
sensitivity of such experiments to $\chi$s can reach beyond $Q_{\chi}\sim
10^{-3}e$ for a wide mass range from a few $\textrm{MeV}/\textrm{c}^{2}$ to a
couple of $\textrm{GeV}/\textrm{c}^{2}$. We propose an experiment, SUBMET
(SUB-Millicharge ExperimenT), which utilizes the $30$ GeV J-PARC proton beam
to search for $\chi$s in this unexplored region.

##  2 Production of millicharged particles at J-PARC

Figure 1: Expected number of $\chi$s that reach a $0.5\textrm{ m}^{2}$
detector located at 280 m from the target. $\textrm{N}_{\textrm{POT}}=10^{22}$
is assumed.

We consider $\chi$s with charge $Q$ produced from $\pi^{0},\eta$ and $J/\psi$
neutral meson decays. The $\Upsilon$ production is irrelevant because the
center-of-mass energy $\sqrt{s}$ is $7.5$ GeV for the collisions between the
$30$ GeV proton beam and the fixed target. The lighter mesons
($\mathfrak{m}=\pi^{0},\eta$) decay through photons
($\pi^{0},\eta\to\gamma\chi\bar{\chi}$), while $J/\psi$ mostly decay to
$\chi$s directly ($J/\psi\to\chi\bar{\chi}$). In both cases, $m_{\chi}$ up to
$m_{\mathfrak{m}}/2$ is kinematically allowed. The number of produced $\chi$s
$N_{\chi}$ can be calculated by the equation in [8],

$\displaystyle N_{\chi}\propto
c_{\mathfrak{m}}\epsilon^{2}\textrm{N}_{\textrm{POT}}\times
f(\frac{m^{2}_{\chi}}{m^{2}_{\mathfrak{m}}})$ (2)

where $c_{\mathfrak{m}}$ is the number of mesons produced per proton-on-target
(POT), $\textrm{N}_{\textrm{POT}}$ is the total number of POT, $\epsilon=Q/e$,
and $f$ is a phase space related integral that goes to $0$ as $m_{\chi}$
approaches $m_{\mathfrak{m}}/2$. Figure 1 shows the expected number of $\chi$s
that reach the detector with face area of $0.5~{}\textrm{m}^{2}$, located
$280$ m from the target. Assuming $\textrm{N}_{\textrm{POT}}=10^{22}$ that
corresponds to running the experiment for $\sim 3$ years, the expected number
of $\chi$s is in the order of $10^{16}$ at $\epsilon=1$ and it is in the order
of $10^{9}$ at $\epsilon=10^{-4}$.

##  3 Site and detector concept

In J-PARC, a $30$ GeV proton beam is incident on a graphite target to produce
hadrons that subsequently decay to a muon and muon neutrino in the decay
volume. The remaining hadrons are then dumped in the beam dump facility. Since
they are MIPs, muons penetrate the beam dump and are identified by the muon
monitor located behind the beam dump facility. The on-axis near detector
(INGRID) is inside the Neutrino Monitor (NM) building located at $280$ m from
the target. The space between the muon monitor and INGRID is filled with sand.

Figure 2: Illustration of the experimental site. $\chi$s are produced near the
target and reach SUBMET after penetrating the beam dump, the muon monitor and
the sand.

If $\chi$s are produced, because of their low interaction probability, they
will penetrate the beam dump facility, the muon monitor, and the sand without
a significant energy loss. Therefore, they can be detected at the NM building
if a detector sensitive to identifying such particles is installed. We
consider the area behind the V-INGRID on B2 ($\sim 30$ m underground) as a
potential detector site. The distance from the axis of the neutrino beam is
$\sim 5$ m. As described below, the signal acceptance is only slightly smaller
than the on-axis location.

Figure 3: Demonstration of the SUBMET detector. There are four layers of
scintillator bars (blue). At one end of each bar a PMT (black) is attached A
$\chi$ will penetrate 4 layers in a narrow time window.

The detector concept proposed here is based on a similar proposal made in [9],
namely, using a segmented detector with large scintillator bars as the core
detector technique. In order to be sensitive to charges down to $10^{-3}e$, a
thick sensitive volume is needed. It is advantageous to segment the large
volume because it reduces backgrounds due to dark currents and shower
particles from cosmogenic muons to a negligible level. It also allows for
utilizing the directionality of the incident $\chi$s to further suppress non-
pointing particles. The detector, as shown in Figure 3, is composed of 4
layers of stacked $5\times 5\times 80$ $\textrm{cm}^{3}$ plastic scintillator
bars. They are aligned such that the produced $\chi$s pass through all layers
in a narrow time window. In each layer there are $14\times 14$ scintillator
bars, so the area of the detector face is about $0.5\textrm{ m}^{2}$. A
prototype of a detector with a similar design has been installed at the LHC,
and shown robustness and sensitivity to $\chi$s [5].

At the end of each scintillator bar a photodetector is attached to convert the
photons to an electronic signal. There are several options, but we consider
Photomultipliers (PMTs), because of their large area coverage, low cost, and
low dark current. Since the size of a bar is $5\times 5$ $\textrm{cm}^{2}$, a
fast commercial $50$ mm diameter PMT is suitable. Including PMTs, the total
volume of the detector is approximately $0.7\times 0.7\times 4$
$\textrm{m}^{3}$.

Figure 4: Signal acceptance rate at 280 m from the target. One box corresponds
to $0.5\times 0.5$ $\textrm{m}^{2}$.

The signal acceptance rate, the fraction of $\chi$s that go into the detector
area of $0.5\textrm{ m}^{2}$ at $280$ m from the target, is calculated. Since
the photons are massless, the direction of the photons and that of $\chi$s
should be aligned. Therefore, the direction of photons from $\pi^{0}$’s is
used as a proxy for the direction of $\chi$s. Figure 4 shows the acceptance
rate in $8\times 8$ $\textrm{m}^{2}$ plane where the axis of the neutrino beam
is at the center. The acceptance rate does not depend on the location within
this range strongly, i.e., near the corners the rate is only $\sim 10$% lower
than the central region because of the large distance from the target. This
provides some flexibility in selecting the location of the detector.

The operation of SUBMET does not interfere with the operation of the neutrino
program and dedicated beam time will not be required. The detector is located
behind the V-INGRID detector, so we don’t expect any interference with the
measurements of the neutrino beam properties and the operation of the INGRID
detector. The detector will be triggered by its own trigger system and the
experiment will have an independent data stream. Therefore, the operation of
this new experiment does not have any conflict with the existing programs at
J-PARC.

##  4 Background sources

The $\chi$s that reach the detector will go through all layers within a $\sim
20$ ns time window (coincidence signal). In this section, the background
sources that can mimic this coincidence signal are discussed. They can be
divided into three categories, random coincidence, beam-induced, and cosmic-
induced backgrounds.

In PMTs, spurious current pulses can be produced by thermal electrons
liberated from the photocathode. If the rate of such pulses (dark count rate,
DCR) is high, the rate of random coincidence can be substantially large even
if the time window for coincidence signal is $20$ ns. The random coincidence
rate can be calculated by $R=nN^{n}\tau^{n-1}$ where $n$ is the number of
layers, $N$ is the DCR, and $\tau$ is the coincidence time window. Using a
typical DCR of $500$ Hz, $n=4$, and $\tau=20$ ns, the random coincidence rate
is $6\times 10^{-5}$ per year. Since there can be $14\times 14=196$ such
coincidence signals, the total coincidence rate is $\sim 0.01$ per year.

Muons are produced from the pion decays in the decay volume together with
neutrinos. The density of quartz, which typically takes up the largest
fraction of sand, is $2.65$ $\textrm{g}/\textrm{cm}^{3}$ and
$dE/dx=1.699\textrm{ MeV}/\textrm{cm}^{2}\textrm{g}$ [10], so the energy loss
of a MIP in $>100$ m of sand is much larger than $30$ GeV. Therefore, such
beam-induced muons can’t reach the detector. Although the muons from the pion
decays can’t reach the detector, neutrinos can and may interact with the
scintillator material to produce small signals. The number of neutrino
interaction events in INGRID is $\sim 1.5\times 10^{8}$ for
$\textrm{N}_{\textrm{POT}}=10^{22}$ [11]. Since a large fraction of INGRID
material is iron, we can use the rate of INGRID as a upper bound for SUBMET.
One layer of SUBMET is approximately $30$ times smaller, so the rate is $\sim
5\times 10^{6}$ for $\textrm{N}_{\textrm{POT}}=10^{22}$ in one layer of
SUBMET. Requiring coincidence in 4 layers, the expected number of this
background is negligible. The interaction of the neutrinos and the material of
the wall of the NM building in front of the detector can produce muons that go
through the detector. These muons can be identified and rejected by installing
scintillator plates between the wall and the detector or by using the very
large scintillation yield of a muon that can be separated from millicharge
signal.

Cosmic muons that penetrate the cavern or the materials above the detector can
produce a shower of particles that is large enough to hit all layers
simultaneously. In such events, the hits in multiple layers can be within the
coincidence time window and will look like a signal event. The particles in
the shower generate more photons than $\chi$s, so the signals from cosmic muon
showers can be rejected by vetoing large pulses. As done to tag the muons
produced in the wall of the NM building in front of the detector, scintillator
plates can be installed covering the whole detector to tag any ordinary-charge
particles or photons incident from top and sides of the detector. These
auxiliary components were proven to be effective in rejecting events with such
particles [5]. In addition, the cosmic shower penetrates the detector
sideways, leaving hits in multiple bars in the same layer, while $\chi$s will
cause a smaller number of hits. A cosmic shower and signals form radioactive
decays overlapping with dark current can be another source of background.
Since the rate of this background depends on the environment strongly, a
precise measurement can be performed in situ only.

To estimate the sensitivity of the experiment, we assume that the total
background over three years of running is $O(1)$ event.

##  5 Sensitivity

The probability of detecting a $\chi$ in a $n$-layer detector is given by
Poisson distribution $P=(1-e^{-\textrm{N}_{\textrm{PE}}})^{n}$ where
$\textrm{N}_{\textrm{PE}}$ is the number of photoelectrons.
$\textrm{N}_{\textrm{PE}}$ is proportional to the quantum efficiency (QE) of
PMT, $\epsilon^{2}$, and the number of photons that reach the end of the
scintillator ($N_{\gamma}$). The $\epsilon^{2}$ term comes from the fact that
the energy loss of a charged particle in matter is proportional to
$\epsilon^{2}$. In order to get $N_{\gamma}$ a GEANT4 [12] simulation is
performed. Using a $5\times 5\times 80$ $\textrm{ cm}^{3}$ Saint-Gobain BC-408
plastic scintillator with surface reflectivity of 98%, we get
$N_{\gamma}=6.3\times 10^{5}$. Using $\textrm{QE}=30$% and taking the area of
the PMT into account, we get $\textrm{N}_{\textrm{PE}}=1.3\times
10^{5}\epsilon^{2}$. We can now calculate the total number of signal events
$s=N_{\chi}P$ measured by the detector. Figure 5 shows the number of detected
signal events $s$ in the plane of $\epsilon$ and $m_{\chi}$.

Figure 5: Expected number of detected signal events $s$ in
$\textrm{N}_{\textrm{POT}}=10^{22}$. There is a sharp drop in
$\epsilon<10^{-3}$ due to small $\textrm{N}_{\textrm{PE}}$.

Figure 6 shows the $95$% CL exclusion curve for
$\textrm{N}_{\textrm{POT}}=10^{22}$. SUBMET provides exclusion down to
$\epsilon=3\times 10^{-4}$ in $m_{\chi}<0.2$ $\textrm{GeV}/\textrm{c}^{2}$ and
$\epsilon=1.5\times 10^{-3}$ in $m_{\chi}<1.6$ $\textrm{GeV}/\textrm{c}^{2}$.
We do not consider a systematic uncertainty on $b$ because it does not have a
significant impact on the exclusion limit (Table 1). Sensitivity of SUBMET is
comparable with that of FerMINI in $m_{\chi}<0.2$
$\textrm{GeV}/\textrm{c}^{2}$. Above that mass range the performance of SUBMET
is worse because the production rate of $J/\psi$ is much lower with the $30$
GeV proton beam.

As shown in Fig 5, the number of signal events recorded by the detector drops
rapidly in $\epsilon<10^{-3}$ due to small $\textrm{N}_{\textrm{PE}}$.
Therefore, one can expect that increasing $\textrm{N}_{\textrm{PE}}$ or
$N_{\chi}$ would not have a large impact on the sensitivity. This will be
discussed in a quantitative way in the next section.

Figure 6: Exclusion at $95$% CL for $\textrm{N}_{\textrm{POT}}=10^{22}$. For
comparison, the constraints from previous experiments are shown as shaded area
and the expected sensitivity of the proposed experiments are drawn in gray
dotted line. Note that FerMINI curve is from 3-layer detector design.

##  6 Alternative detector design

In this section, we discuss how sensitivity to sub-millicharged particles
depends on the configuration of the detector. If further optimization of the
detector is needed, we can use this quantitative study as a guidance.

Silicon Photomultipliers (SiPMs) can be considered instead of PMTs because of
its excellent sensitivity to single photons and small size. A challenge is
that they have a typical DCR of $5\textrm{ MHz}/\textrm{cm}^{2}$ at room
temperature, about three orders of magnitude larger than that of PMTs. This
leads to a significantly larger random coincidence rate. Since the dark
current is a thermal effect, it has a strong dependence on temperature.
Typical reduction rate is a factor of $\sim 2$ per $10^{\,\circ}$C. So, to
suppress the DCR to the similar level of PMT’s, the temperature should be
lowered by of order of $100^{\,\circ}$C. In the scenario where SiPMs are used,
the size of scintillator bars is chosen to be $1\times 1$ cm due to difficulty
in covering the large area of the scintillator. In this case the light
attenuation effect increases due to more reflections and the photon yields
decreases by a factor of $2$.

Searches in the sub-millicharge regime rely on the Poissonian fluctuation of
small $\textrm{N}_{\textrm{PE}}$. Approximating
$P\simeq\textrm{N}_{\textrm{PE}}$ to the first order for small $\epsilon$, we
arrive at the following relation

$\displaystyle\textrm{N}_{\epsilon=1,\chi}\epsilon^{2}(\textrm{N}_{\epsilon=1,\textrm{PE}}\epsilon^{2})^{n}\geq
s^{\textrm{95\%}}$ (3)

where the subscript $\epsilon=1$ refers to the values at $\epsilon=1$ and
$s^{\textrm{95\%}}$ is the number of signal events that provides 95% exclusion
limit. Reordering in terms of $\epsilon$, we recognize the exclusion limit at
95% CL is

$\displaystyle\epsilon=\left(\frac{\textrm{N}_{\epsilon=1,\textrm{PE}}^{n}\textrm{N}_{\epsilon=1,\chi}}{s^{\textrm{95\%}}}\right)^{-\frac{1}{2n+2}}.$
(4)

Not surprisingly, the $\epsilon$-power term induces a sharp cutoff in signal
at $\epsilon\sim O(10^{-4})$ as shown in Figure 5. This limits the sensitivity
to the same regime regardless of the detector configuration.

$N_{\chi}$(relative) |  $\textrm{N}_{\textrm{PE}}$(relative) | $b$ | Exclusion limit on $\epsilon$ for $m_{\chi}=10$ $\textrm{MeV}/\textrm{c}^{2}$
---|---|---|---
1 | 1 | 1 | $2.6\times 10^{-4}$
2 | 1 | 1 | $2.4\times 10^{-4}$
1 | 1 | 100 | $3.0\times 10^{-4}$
1 | 2 | 1 | $1.9\times 10^{-4}$
1 | 3 | 1 | $1.7\times 10^{-4}$
Table 1: Various detector configurations and their sensitivity. $N_{\chi}$
(relative) is the yields of signal events within the acceptance relative to
the baseline, $\textrm{N}_{\textrm{PE}}$ (relative) is the number of photon
electrons relative to the baseline $\textrm{N}_{\textrm{PE}}=1.3\times
10^{5}$, and $b$ is the background yields.

Table 1 shows different detector configurations and the corresponding
exclusion limit for $m_{\chi}=10$ $\textrm{MeV}/\textrm{c}^{2}$. The default
configuration is in the first row, $N_{\chi}$(relative)=1,
$\textrm{N}_{\textrm{PE}}$(relative)=1, $b=1$, and $n=4$ where “relative”
means relative to the baseline.

The improvements we can achieve by extending the duration of data collection
or making the detector area larger (adding more bars to each layer) is very
modest. Extending the duration or the detector area by a factor of 2 increases
$N_{\chi}$ by the same amount. Making an assumption that the background
remains same, the sensitivity is only improved by $10$% ($2^{\textrm{nd}}$ row
in Table 1).

The impact of $b$ is also small. If $b$ is increased by a factor of $100$, the
limit is degraded by $15$% ($3^{\textrm{rd}}$ row in Table 1).

As the $4^{\textrm{th}}$ and $5^{\textrm{th}}$ rows in Table 1 show, the most
effective component to enhance sensitivity is $\textrm{N}_{\textrm{PE}}$.
Increasing $\textrm{N}_{\textrm{PE}}$ by a factor of $2(3)$ improves the
sensitivity limit by $25(35)$%. This can be achieved by using an inorganic
crystal with large Z such as BGO or by making the scintillator bars longer.
GEANT4 simulation shows that a $5\times 5\times 20$ $\textrm{ cm}^{3}$ BGO bar
provides $N_{\gamma}=1.2\times 10^{6}$, about a factor of $2$ larger than the
$5\times 5\times 80$ $\textrm{ cm}^{3}$ BC-408 scintillator. Further
increasing the length of the bar to $40\textrm{ cm}$, one can achieve
$N_{\gamma}=1.7\times 10^{6}$. However, cost should be considered as well
because BGO is less cost-effective than plastic scintillators.

Though $\textrm{N}_{\textrm{PE}}$ plays the main role in enhancing sensitivity
to sub-millicharged particles, the exclusion limit is still in the range of
$\epsilon=1.5-3.0\times 10^{-4}$ with the variations considered. This
indicates that the sensitivity does not strongly depend on the configuration
of the detector.

With 3 layers the number of detected signal events is larger because
$P=(1-e^{-\textrm{N}_{\textrm{PE}}})^{3}$ instead of
$(1-e^{-\textrm{N}_{\textrm{PE}}})^{4}$. However, the background increases as
well; the random coincidence rate is approximately $3000$ with 3 layers. Using
$n=3$ and $b=3000$, the exclusion limit is slightly improved for
$m_{\chi}<0.2$ $\textrm{GeV}/\textrm{c}^{2}$ as shown in Fig 6. However, the
experience from the milliQan demonstrator [5] gives a lesson that the random
coincidence and the correlated backgrounds (induced by either cosmic showers
or radiation within the scintillator bars) overlapping with dark current are
potentially the major background sources. Adding one more layer (4-layer
configuration) helps understanding and rejecting these backgrounds so that the
experiment is carried out in a low background environment. This will make the
observation of a potential signal more robust.

##  7 Summary

We studied the feasibility of an experiment searching for millicharged
particles using 30 GeV proton fixed-target collisions at J-PARC. With
$\textrm{N}_{\textrm{POT}}=10^{22}$, the experiment provides a sensitivity to
$\chi$s with charge down to $\epsilon\simeq 3\times 10^{-4}$ in $m_{\chi}<0.2$
$\textrm{GeV}/\textrm{c}^{2}$ and $\epsilon\simeq 1.5\times 10^{-3}$ in
$m_{\chi}<1.6$ $\textrm{GeV}/\textrm{c}^{2}$; this is the regime largely
uncovered by the previous experiments. We also explored a few detector designs
to achieve an optimal sensitivity to $\chi$s. The photoelectron yields is the
main driver, but the sensitivity does not have a strong dependence on the
detector configurations in the sub-millicharge regime.

## Acknowledgements

We thank Tsutomu Mibe, Yoshiaki Fujii, Takeshi Nakadaira, and Toshifumi
Tsukamoto for the useful discussions on the detector site. In particular, we
thank Toshifumi Tsukamoto for taking photographs of the Neutrino Monitor
building so that we understand the spatial constraints inside the building. We
also thank Hong Joo Kim for the discussion on the property of various
scintillation materials.

## References

  * [1] Bob Holdom.  Two U(1)’s and Epsilon Charge Shifts.  Phys. Lett. B, 166:196, 1986.
  * [2] A Prinz, R Baggs, J Ballam, S D Ecklund, C Fertig, J A Jaros, K Kase, A Kulikov, W G J Langeveld, R Leonard, T P Marvin, T Nakashima, W R Nelson, A C Odian, M Pertsova, G Putallaz, and A” Weinstein.  Search for Millicharged Particles at SLAC.  Phys. Rev. Lett., 81:1175, Aug 1998.
  * [3] CMS Collaboration.  Search for fractionally charged particles in pp collisions at $\sqrt{s}=7$ TeV.  Phys. Rev. D, 87, 2013.
  * [4] CMS Collaboration.  Searches for long-lived charged particles in pp collisions at $\sqrt{s}$ = 7 and 8 TeV.  J. High Energy Phys., 07:122, 2013.
  * [5] A. Ball et al.  Search for millicharged particles in proton-proton collisions at $\sqrt{s}=13$ TeV.  5 2020.
  * [6] S. Davidson, B. Campbell, and David C. Bailey.  Limits on particles of small electric charge.  Phys. Rev. D, 43:2314, 1991.
  * [7] Timon Emken, Rouven Essig, Chris Kouvaris, and Mukul Sholapurkar.  Direct detection of strongly interacting sub-gev dark matter via electron recoils.  J. Cosmol. Astropart. Phys., 09:070, 2019.
  * [8] Kevin J. Kelly and Yu-Dai Tsai.  Proton fixed-target scintillation experiment to search for millicharged dark matter.  Phys. Rev. D, 100:015043, Jul 2019.
  * [9] Austin Ball et al.  A Letter of Intent to Install a milli-charged Particle Detector at LHC P5.  7 2016.
  * [10] Donald E. Groom, Nikolai V. Mokhov, and Sergei I. Striganov.  Muon stopping power and range tables 10-MeV to 100-TeV.  Atom. Data Nucl. Data Tabl., 78:183–356, 2001.
  * [11] K. Abe et al.  Measurements of the T2K neutrino beam properties using the INGRID on-axis near detector.  Nucl. Instrum. Meth. A, 694:211–223, 2012.
  * [12] S. Agostinelli et al.  GEANT4: A Simulation toolkit.  Nucl. Instrum. Meth. A, 506:250–303, 2003.