# Search for dark energy with neutron interferometry

Hauke Fischer  Atominstitut, Technische Universität Wien, Stadionallee 2,
A-1020 Vienna, Austria  Christian Käding  Atominstitut, Technische Universität
Wien, Stadionallee 2, A-1020 Vienna, Austria  Hartmut Lemmel  Atominstitut,
Technische Universität Wien, Stadionallee 2, A-1020 Vienna, Austria  Institut
Laue Langevin, 38000 Grenoble, France  Stephan Sponar  Atominstitut,
Technische Universität Wien, Stadionallee 2, A-1020 Vienna, Austria  Mario
Pitschmann<EMAIL_ADDRESS>Atominstitut, Technische Universität
Wien, Stadionallee 2, A-1020 Vienna, Austria

###### Abstract

We use previously obtained experimental results by neutron interferometry to
effectively constrain the parameter space of several prominent dark energy
models. This investigation encompasses the environment-dependent dilaton
field, a compelling contender for dark energy that emerges naturally within
the strong coupling limit of string theory, alongside symmetron and chameleon
fields. Our study presents substantial improvements over previous constraints
of the dilaton and symmetron fields, improving parameter constraints by
several orders of magnitude. However, the analysis does not yield any new
constraints on the chameleon field. Furthermore, we establish constraints for
the projected neutron split interferometer, which has recently concluded a
decisive proof-of-principle demonstration. Our symmetron simulations reveal
that depending on the parameter values there are multiple static solutions
with increasing number of nodes and increasing energy inside a cylindrical
vacuum chamber. This agrees with results obtained earlier in the literature
for infinitely parallel plates. Interestingly, while these multiple solutions
can correspond to domain walls forming inside the vacuum chamber, we also find
solutions that do not reach their vacuum expectation value inside the vacuum
chamber, but display multiple nodes nonetheless.

##  I Introduction

The origin of dark energy is one of the greatest unsolved problems in modern
physics. Observations on the cosmic scale unveil an ongoing acceleration in
the expansion of our Universe [1, 2, 3]. An unknown substance, called dark
energy, that fills the Universe is the most prominent explanation for this
acceleration. While introducing a cosmological constant into the framework of
general relativity can in principle explain an accelerated expansion, such an
approach grapples with intricate issues of fine-tuning [4]. For this reason,
new hypothetical scalar fields have been postulated, which couple to gravity
and can account for dark energy [5]. These scalar fields typically induce so-
called fifth forces. To avoid conflict with ongoing fifth force searches, they
must exhibit some screening mechanism in order to evade detection. Several
screening mechanisms have been theorized, including the chameleon [6, 7],
K-mouflage [8, 9], Vainshtein [10] and Damour-Polyakov mechanisms [11].

The analysis in this article focuses on the dilaton [12, 13], symmetron [14,
15, 16] and chameleon field [17].

The exponentially decreasing potential of the dilaton is expected from the
strong coupling limit of string theory [18, 19, 20] and has been suggested as
a dark energy candidate. First, experimental constraints from Lunar Laser
Ranging [21] and gravity resonance spectroscopy [22, 23] as well as projective
constraints from an experiment measuring the Casimir effect [24] have been
derived in [25, 26]. Furthermore, [27] has shown that additional constraints
can be obtained by studying the dilaton-induced open quantum dynamics based on
the formalism developed in [28, 29, 30]. Symmetrons employ spontanous symmetry
breaking similar to the Standard Model Higgs [31, 32]. In low density regions
symmetrons can obtain a non-vanishing vacuum expectation value (VEV) leading
to the presence of a fifth force, while the symmetry is restored in high
density regions. Recently, symmetron fifth forces have even been suggested as
alternative explanations for the observed effects otherwise attributed to
particle dark matter [33, 34, 35, 36]. The symmetron field has already been
constrained by several table top experiments, such as atomic interferometry
[37, 38, 39, 40], Eötwash experiments [41], with gravity resonance
spectroscopy [15, 42, 16, 43] and by atomic precision measurements [44].
Chameleons [45] are known to screen by increasing their masses in dense
environments and have already been investigated and constrained
comprehensively [17] (see e.g. also [46]). Furthermore, [47] proposed that
tightened constraints on both symmetrons and chameleons could be obtained from
Bose-Einstein condensate (BEC) interferometry in the future. Similarly, it was
suggested that popular scalar field models with screening mechanisms,
including those discussed in the present article, could be investigated in
BEC-based analogue gravity simulations [48].

A common feature among these models is that they are screened in dense
environments and that forces on large objects are suppressed since the fields
only couple to a thin-shell below the surface of an object [17]. As a result,
forces on macroscopic objects as e.g. in the Solar System are typically very
weak, while small objects such as neutrons in a vacuum chamber are ideal tools
to probe these fields. Neutrons propagating through a vacuum chamber would
experience a significant phase shift, and the absence thereof can be used to
set stringent constraints on these models. In [49, 50] neutron interferometry
has already been used to constrain the chameleon model for a few parameters
(albeit these constraint have been superseded [17]). We extend the analysis to
the dilaton and symmetron field and to a larger part of the chameleon
parameter space and derive improvements that can be made through a split
crystal interferometer [51]. It should be noted that ultra-cold neutrons
provide versatile tools in the search for new physics in addition to scalar
fields with screening mechanisms, see e.g. [52, 53, 54, 55, 56, 57, 58, 59,
60].

In section II, we provide the theoretical background for the investigated
scalar fields as well as the expected induced phase shift on neutrons in a
vacuum chamber. Details on the experimental setup are given in section III.
This is followed by sections IV, V and VI describing the dilaton, symmetron
and chameleon constraints, respectively. We conclude our findings in VII.
Finally, in appendices A and B, we provide details on the numerical solutions
of differential equations and derive the phase shift formula, respectively.

##  II Theoretical background

The effective potential of the scalar fields under consideration is given by

$\displaystyle V_{\text{eff}}(\phi;\rho)=V(\phi)+\rho A(\phi)\>,$ (1)

where $V(\phi)$ describes the self-interactions of the field and the Weyl
factor $A(\phi)$ couples to the ambient matter density $\rho$. The dilaton
(D), symmetron (S) and chameleon (C) models are defined by [17, 61]

$\displaystyle V_{D}(\phi)$
$\displaystyle=V_{0}\,e^{-\lambda_{D}\phi/m_{\text{Pl}}}\>,$ (2)
$\displaystyle V_{S}(\phi)$
$\displaystyle=-\frac{\mu^{2}}{2}\,\phi^{2}+\frac{\lambda_{S}}{4}\,\phi^{4}\>,$
(3) $\displaystyle V_{C}(\phi)$
$\displaystyle=\frac{\Lambda^{n+4}}{\phi^{n}}\>,$ (4)

together with

$\displaystyle A_{D}(\phi)$
$\displaystyle=1+\frac{A_{2}}{2}\frac{\phi^{2}}{m^{2}_{\text{Pl}}}\>,$ (5)
$\displaystyle A_{S}(\phi)$ $\displaystyle=1+\frac{\phi^{2}}{2M^{2}}\>,$ (6)
$\displaystyle A_{C}(\phi)$ $\displaystyle=e^{\phi/M_{c}}\simeq
1+\frac{\phi}{M_{c}}\>.$ (7)

The parameters of the dilaton are the energy scale $V_{0}$, related to the
dark energy of the Universe, the dimensionless parameter $\lambda_{D}$ as well
as the dimensionless coupling parameter $A_{2}$, while $m_{\text{Pl}}$ denotes
the reduced Planck mass. The symmetron parameters are given by the tachyonic
mass $\mu$, a dimensionless self-coupling constant $\lambda_{S}$ and $M$ as a
coupling constant to matter with the dimension of mass. For chameleon models,
$n\in\mathbb{Z}^{+}\cup 2\mathbb{Z}^{-}\setminus${-2} determines the power of
the self-interaction potential, $\Lambda$ defines an energy scale that is
sometimes related to the dark energy of the Universe and $M_{c}$ is a coupling
constant with the dimension of mass. To neglect possible couplings to matter
of higher order, we restrict our analysis to

$\displaystyle\frac{A_{2}}{2}\frac{\phi^{2}}{m^{2}_{\text{Pl}}},\frac{\phi^{2}}{2M^{2}},\frac{\phi}{M_{c}}\ll
1\>,$ (8)

which leads to sharp cut-offs in the displayed exclusion plots. For the
chameleon field, this restriction is actually not necessary. However, regions
where neutron interferometry is sensitive and this condition does not hold
have already been constrained by other experiments.

The experimental setup is such that the neutron beam is split into two paths,
in each the neutrons traverse a chamber containing some gas. For brevity, the
chamber containing gas of significantly lower density will be denoted as
”vacuum chamber”, while the other chamber will be referred to as ”filled with
air”. This leads to a relative phase difference as will be detailed in section
III. The presence of a scalar field inside a chamber induces an additional
phase shift, which to leading order is given by

$\displaystyle\delta\varphi_{X}=-\frac{m_{n}}{k_{0}}\int\limits_{\text{CFP}}U_{X}(\bm{x})\,ds\>,$
(9)

with $X\in\\{D,S,C\\}$, the integration extending over the classical flight
path (CFP) and

$\displaystyle U_{D}(\bm{x})$
$\displaystyle=\mathfrak{Q}_{D}\frac{A_{2}m_{n}}{2m^{2}_{\text{Pl}}}\left(\phi^{2}(\bm{x})-\phi^{2}_{\text{Air}}\right),$
(10) $\displaystyle U_{S}(\bm{x})$
$\displaystyle=\mathfrak{Q}_{S}\frac{m_{n}}{2M^{2}}\left(\phi^{2}(\bm{x})-\phi^{2}_{\text{Air}}\right),$
(11) $\displaystyle U_{C}(\bm{x})$
$\displaystyle=\mathfrak{Q}_{C}\frac{m_{n}}{M_{c}}\left(\phi(\bm{x})-\phi_{\text{Air}}\right).$
(12)

Here, $k_{0}$ is the wave number of the neutron, $m_{n}$ its mass,
$\phi_{\text{Air}}$ is the field value that minimizes the potential in air and
$\mathfrak{Q}_{X}$ the screening charge of the neutron. We provide a
derivation of these formulas and its applicability conditions in Appendix B.
Computing the field of the neutron and the chamber is in general a two-body
problem between the neutron and the vacuum or air chamber. We approximate this
problem by treating the neutron as a classical massive sphere with radius 0.5
fm in agreement with QCD (this approximation is refered to as ”fermi-
screening”), which allows the screening of the neutron to be taken into
account approximately. Thereby, the potential is multiplied with a screening
charge $\mathfrak{Q}$ that takes values between 0 (for fully screened spheres)
and 1 (for unscreened spheres). The analytical expressions used can be found
in [25] for the dilaton and in [15, 16] for the symmetron model. For the
chameleon field, we derived a screening charge by following analogous
reasoning as for the dilaton screening charge. We recall its derivation in
Appendix C. In [26] a second approximation (”micron-screening”) was
considered, by extracting a radius of the neutron from the vertical extent of
the wave function, which amounts to 5.9 $\mu$m in that experiment. The wave
function in the neutron interferometer has a very particular shape, being
extended over mm in one direction and nm in the others, cf. section III. Hence
we refrain from trying to extract a radius from it. However, it was verified
that assuming a neutron radius larger than 0.5 fm, as could be extracted from
any extent of the wave function, would lead to constraints containing the
fermi-screening ones. Hence, fermi-screening constraints are more
conservative.

Figure 1: Schematic illustration of neutron interferometric setups (not to
scale) to test various prominent dark energy models using vacuum and air
chambers in cylinder design. (a) Monolithic interferometer as used in [49] and
(b) proposed split-crystal interferometer setup with longer interaction region
and larger beam separation, allowing for larger cylinder diameters.

For each field model it is necessary to compute $\phi$ inside the vacuum and
air chamber, which amounts to solving the non-linear differential equation

$\displaystyle\Delta\phi(\bm{x})=V_{\text{eff},\phi}(\phi(\bm{x});\rho(\bm{x}))\>.$
(13)

As boundary condition we demand that $\phi$ reaches its potential minimum
$\phi_{M}$ inside the shell of the vacuum and air chamber, and restrict our
analysis to parameters where this assumption is justified. Due to the absence
of analytical solutions, the differential equations were solved numerically.
Details about the algorithm employed are provided in Appendix A.

Since every parameter combination requires a separate numerical solution, the
computational cost of constraining a 3D parameter space in this way is
unfeasible. Therefore, the rectangular vacuum and air chamber have been
approximated by cylinders, which amounts to solving

$\displaystyle\left[\frac{1}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial}{\partial r}\right)+\frac{\partial^{2}}{\partial
z^{2}}\right]\phi(r,z)=V_{\text{eff},\phi}(\phi;\rho)\>.$ (14)

These cylinders have the same length as the real chambers (which was 0.094 m
in the experiment, 0.5 m was assumed for the split interferometer) and either
just fit into the real chamber ($r=2$ cm in both chambers) or are just large
enough to encompass the real chamber ($r=\sqrt{8}$ cm in the vacuum chamber
and $r=\sqrt{3.25^{2}+2^{2}}$ cm in the air chamber). That the difference
between both geometries is negligible has been verified and hence the final
constraints have been computed assuming $r=2$ cm. For the future split
interferometer a chamber with cylinder symmetry and radius $r=4.75$ cm has
been assumed.

### Constraint calculation profile mode

In the experiment, two measurement modes were applied. First, in the profile
mode, the pressure in the vacuum chamber was fixed to $10^{-4}$ mbar and the
beam position varied. Initially, both beams passed through the center of
either chamber. Next, the beam position was changed to a position closer to
the chamber walls corresponding to a maximum displacement of $\sim$1.5 cm from
the center. Since scalar fields are more suppressed close to the chamber walls
than at the center of the chamber, the induced phase shift would be position
dependent. The phase difference is measured for both positions, and the
difference can be used to derive experimental constraints for scalar fields.
The experiment constrains phase shifts exceeding 5.5 degrees at a confidence
level of 95%. The phase shift was computed using

$\displaystyle\delta\varphi_{X,P}=-\frac{m_{n}}{k_{0}}\int_{-L/2}^{L/2}dz\,\big{(}U_{X}(0,z)-U_{X}(0.015\,\text{m},z)\big{)}\>,$
(15)

for both the air and vacuum chamber. Parameter constraints are obtained by
using

$\displaystyle\big{|}\delta\varphi_{X,P}(\text{Air})-\delta\varphi_{X,P}(\text{Vacuum})\big{|}\geq
5.5^{\circ}\>.$ (16)

Projective constraints for the split interferometer have been derived by the
same equations and experimental parameters with the exception that the second
beam position was assumed to be 4 cm displaced from the center.

### Constraint calculation pressure mode

In the pressure mode, both beams passed through the center of either chamber
but the pressure inside the vacuum chamber was varied. The phase shift
resulting from the highest pressure of $P_{0}=10^{-2}$ mbar was used as a
reference and the pressure lowered to $P_{1}=2.4\times 10^{-4}$ mbar. This
measurement mode constrains phase shifts exceeding 3.5 degrees with a
confidence level of 95%. The vacuum chamber was simulated for both pressures
and parameter constraints are obtained by using

$\displaystyle\frac{m_{n}}{k_{0}}\bigg{|}\int_{-L/2}^{L/2}dz\,\big{(}U_{X,P_{0}}(0,z)-U_{X,P_{1}}(0,z)\big{)}\bigg{|}\geq
3.5^{\circ}\>.$ (17)

The projective constraints for the split interferometer have been computed
using the same equations and experimental parameters.

##  III Experimental setup

In a neutron interferometer [62, 63] a beam of neutrons is split by amplitude
division into two beam paths which are superposed coherently after passing
through different regions of space (Mach-Zehnder geometry). One of the most
straightforward experiments is the measurement of the neutron index of
refraction of a sample material which allows to determine the coherent neutron
scattering length of the nuclei contained in the sample. In this sense, the
present experiment measures whether the index of refraction of vacuum indeed
vanishes or not.

The neutron beam has a typical diameter of a few mm to cm. The coherence
length of the incoming neutrons is in the order of nm, so the beam consists of
an incoherent ensemble of neutrons spread over the beam cross section. Each
individual neutron, however, is coherently split into the two beam paths,
separated by a few cm, and exhibits self-interference. Therefore, neutron
interferometry has also been used to study the foundations of quantum
mechanics, starting in the early 1970s with the verification of the 4$\pi$
spinor symmetry [64], up to recent experiments on entanglement and weak
values, reviewed e.g. in [56].

Beam splitters and mirrors in neutron interferometry are based on Bragg
diffraction on perfect silicon crystals in transmission (Laue) geometry. This
implies that the neutron wave function is not only coherently split into two
beam paths but also spread to a width of several mm within a single beam path
[63, 65]. It is still subject of discussion whether the full extension of the
wave function or the particle size of the neutron hast to be taken into
account when calculating the scalar field around the neutron. As already
pointed out, we have restricted our calculations to the latter since it gives
the more conservative constraints.

The beam splitters of an interferometer have to be aligned to each other with
nrad and sub-nm accuracy. Therefore, neutron interferometers are usually built
monolithically out of a single crystal, cf. Fig. 1 (a). This design allowed in
our last experiment [49] for a beam separation of 50 mm and a chamber length
of 94 mm. However, we would gain in sensitivity if the vacuum chambers were
longer and could have a larger diameter, allowing for a stronger scalar field
to build up. This could be realized with a split-crystal interferometer, cf.
Fig. 1 (b), which is currently under construction at our neutron
interferometry station S18 at the Institut Laue-Langevin (ILL) in Grenoble. We
envisage a beam separation of 10 cm and a chamber length of 0.5 m. A proof-of-
principle experiment with a small split-crystal interferometer was already
carried out [51].

In either setup the incident neutron beam with mean wavelength $\lambda_{\rm
n}=2.72$ Å ($\delta\lambda_{\rm n}/\lambda_{\rm n}\sim 0.043$, Bragg angle
45∘) and $4\times 8\,\mathrm{mm^{2}}$ beam cross section is split and passes
the air or vacuum chamber respectively. After recombination at the last
interferometer plate the intensity in the exiting beams (O in forward and H in
refracted direction) is measured by detectors with an efficiency above 99
percent. The intensity oscillates between the two exits as function of the
phase. By rotating the phase shifter flag a sinusoidal fringe pattern is
recorded. An additional phase created by the vacuum chambers can then be
determined by a shift of the fringes.

One chamber, labelled “vacuum” in Fig. 1, is evacuated to let the scalar field
build up while the other chamber, labelled “air”, is filled with some gas to
suppress the field. This configuration creates primarily a phase shift due to
the gas index of refraction, which is in the order of 100° to 1000° at ambient
pressure. By using a gas with a well-known composition, e.g. pure Helium, this
phase shift could be calculated and corrected for to high precision. However,
in our experiment we reduced the pressure in the air chamber from ambient 1000
mbar to 0.01 mbar. While this pressure is high enough to suppress the scalar
field, it is low enough to reduce the gas phase shift to below detection
limit. The pressure in the vacuum chamber is still lower by a few orders of
magnitude. In profile mode, the vacuum chambers are moved sidewards to vary
the distance between the neutron beam and the side walls, thereby probing the
shape of the scalar field within the chamber.

##  IV Derivation of dilaton constraints

The extremely high computational cost of solving the non-linear differential
equation for many parameter values requires to lower the amount of simulations
needed to find the constrained parameter volume (an example of a simulated
dilaton field is given in Fig. 2). While constraints from the profile mode are
stronger for long interaction ranges where the field varies significantly
throughout the chamber, the pressure mode is superior at probing short
interaction ranges. Fields with very short interaction ranges are more
difficult to simulate, since the slopes close to the chamber walls get
arbitrarily steep requiring an ever finer discretization.

Figure 2: Simulated dilaton field for $V_{0}=10$ MeV4, $A_{2}=10^{37}$ and
$\lambda_{D}=10^{-9}$ with $\rho_{V}=7.08\times 10^{-17}$ MeV4.

In this case, however, the following approximation is highly accurate inside
the vacuum chamber (obvious analogous expressions for the air chamber are
omitted)

$\displaystyle\int_{0}^{L}dz\left(\phi^{2}(r=0,z)-\phi^{2}_{\text{Air}}\right)\simeq
L\left(\phi^{2}_{V}-\phi^{2}_{\text{Air}}\right),$ (18)

since the field close to the chamber walls reaches its potential minimum
$\phi_{V}$ quickly. For parameter values corresponding to this case, the
profile mode cannot set constraints anymore, because the field effectively
looses its position dependence. To illustrate this point we define

$\displaystyle\delta_{\text{sim}}$
$\displaystyle:=-\mathfrak{Q}_{D}\frac{A_{2}m_{n}^{2}}{2k_{0}m^{2}_{\text{Pl}}}\int_{0}^{L}dz\left(\phi^{2}(0,z)-\phi^{2}_{\text{Air}}\right),$
(19) $\displaystyle\delta_{\text{approx}}$
$\displaystyle:=-\mathfrak{Q}_{D}\frac{A_{2}m_{n}^{2}}{2k_{0}m^{2}_{\text{Pl}}}\,L\left(\phi^{2}_{V}-\phi^{2}_{\text{Air}}\right).$
(20)

Here, $\delta_{\text{sim}}$ is the actual phase shift computed from a
simulation for a neutron propagating through the center of the chamber
($r=0$), while $\delta_{\text{approx}}$ is an approximation that is valid only
for very short-ranged fields. For the dilaton field, small interaction ranges
correspond to large values of $A_{2}$. We demonstrate in Fig. 3 that the error
in the approximation becomes continuously smaller for larger values of $A_{2}$
and that the field takes on its VEV for an ever larger region inside the
cylinder.

Figure 3: The plots are computed for $\lambda_{D}=10^{-9}$ and $V_{0}=10$
MeV4. Left: The relative error of the approximation as a function of $A_{2}$
is plotted. Right: Normalized field profiles along the $z$ axis for $r=0$. The
value 1 corresponds to $\phi(0,z)=\phi_{V}$, while 0 corresponds to
$\phi(0,z)=\phi_{M}$. The strong slopes occur at the chamber walls.

The constraints were computed as follows. For fixed $V_{0}$ and beginning from
the smallest allowed value of $A_{2}$ (this comes from the long-range cut-off
to ensure that the field decays to $\phi_{M}$ inside the shell of the vacuum
chamber) a search was performed for the contour of the constrained region
along the $\lambda_{D}$-axis with a step width of 0.1 in a logarithmic plot.
Next, $A_{2}$ was increased by a factor of 10 and the procedure repeated until
the regime was entered, where the pressure mode dominates and the
approximation becomes indistinguishable from the simulations. For even larger
values of $A_{2}$, only the approximation was used in order to compute the
rest of the constrained region. Next, the points just at the edge of the
constrained region were connected. An example of this procedure is shown in
Fig. 4.

Figure 4: The blue area corresponds to constraints given by the short-range
approximation Eq. 20. The bullet points are on the real contour of the
constrained region obtained from simulations with $V_{0}=10$ MeV4.

For still larger values of $A_{2}$, a limit is reached where the necessary
resolution of the finite element grid becomes so high that the field can
practically no longer be computed. Therefore, the approximation solves
multiple numerical challenges. The full constraints combining the pressure and
profile mode are shown in Fig. 5 together with already existing constraints
from gravity resonance spectroscopy (qBounce).

Figure 5: Here, the constrained regions obtained from neutron interferometry
for the dilaton field are plotted in color. The red lines correspond to
existing constraints from qBounce [26], assuming a radius of 0.5 fm for the
neutron as well. Dashed lines correspond to improvements expected for the
geometry of a split interferometer (cylindrical vacuum chamber with length of
0.5 meters and radius 4.75 cm). Left: In this part of the parameter space, the
constrained area does not change its shape for increasing $V_{0}$ but only
shifts towards lower $\lambda_{D}$. Right: The constrained areas shift
systematically towards lower values of $A_{2}$ for increasing $V_{0}$ without
changing their shape but become smaller due to a cut-off. This is indicated by
the arrows. The lower cut-off is a short-range cut-off to ensure that the
field decays to its potential minimum inside the cylinder shell, the remaining
cuts are due to the condition given in Eq. 8. The dilaton range is plotted for
$V_{0}=10\ \text{MeV}^{4}$ and $\rho_{V}=7.08\times 10^{-17}$ MeV4.

Lunar Laser Ranging constraints cover much smaller values of $A_{2}$ than
tabletop experiments and, therefore, these constraints are not shown [26].

##  V Derivation of symmetron constraints

The effective potential of the symmetron field [66]

$\displaystyle
V_{\text{eff}}(\phi)=\frac{1}{2}\left(\frac{\rho}{M^{2}}-\mu^{2}\right)\phi^{2}+\frac{\lambda_{S}}{4}\,\phi^{4}\>,$
(21)

allows for spontaneous symmetry breaking. In regions of high density, for
which $\frac{\rho}{M^{2}}-\mu^{2}>0$ holds, there is only one real minimum of
the effective potential at $\phi=0$. The field is in its symmetric phase and
fifth forces vanish. However, in regions of low densities, where
$\frac{\rho}{M^{2}}-\mu^{2}<0$ holds, the field is in its symmetry-broken
phase and obtains a non-vanishing VEV inducing a fifth force. Furthermore, due
to the $\phi\rightarrow-\phi$ symmetry of the potential the differential
equation Eq. 13 can have more than one solution (see [15, 16]). This property
of the field is elucidated in the next subsection. For the calculation of
constraints, we exclusively used solutions with a single node, which
correspond to the lowest energy solution in general (for an example see Fig.
8). We simulated the field numerically and found that for large enough
$\mu$-values no solution exists. This is in agreement with [67, 15], where it
was found that for large enough values of $\mu$ no solution can exist between
two infinitely extended mirrors. Therefore, with the given dimensions of the
vacuum and air chambers for symmetron ranges of $\sim$1 mm and larger no field
solution can exist anymore.

This natural restriction of the symmetron field limits the parameter space to
be probed to short-ranged fields. It was found that using the approximation
$\phi\simeq\phi_{V}$ in the integral for the phase shift is accurate enough in
the whole parameter space where new constraints can be set. In Fig. 6
constraints obtained from neutron interferometry are plotted alongside already
existing constraints published in [41, 44].

Figure 6: The constrained volume for the parameters of the symmetron field
obtained from neutron interferometry is shown in color, as well as existing
constraints from atomic precision tests and the Eöt-Wash experiment. The lower
cut-off comes from the condition given in Eq. 8. At the left edge the symmetry
of the potential is restored inside the vacuum region and the field vanishes
everywhere.

The improvements of the split interferometer are too small to be visible in
this log-log plot. For an alternative presentation with $\mu$ varying
continuously see Fig. 7.

Figure 7: Here, the constrained volume for the parameters of the symmetron
field obtained from neutron interferometry for several fixed values of $M$ is
depicted.

###  V.1 Domain walls and multiple node solutions

Since in the symmetry-broken phase there are two distinct VEVs, the
possibility arises that the symmetron settles to the positive VEV in some
regions of space and to the negative VEV in others. Regions where the field
takes on either VEV are called domains, while the boundaries connecting these
domains are the domain walls (see e.g. [68]). While the latter might be
unstable on cosmological time scales [69], it has been suggested that they
might lead to observable consequences for ultra-cold neutrons. This is because
neutrons would be accelerated towards the domain wall, resulting in a
deflection angle [70]. Furthermore, by draining gas in a vacuum chamber and
allowing the field to relax to either VEV, such domain walls might
deliberately be induced [71]. A search for static domain walls inside a
realistic cylindrical vacuum chamber was performed by trying a large number of
plausible initial guesses for Newton’s method in solving Eq. 14. The spectrum
of solutions found for a fixed parameter combination is shown in Fig. 8.

Figure 8: The field profiles of the three solutions from Table 1 are shown.
For the given parameters $\phi_{V}=0.61$ meV holds.

The top solution has only one domain and one local extremum at the center of
the chamber, while the solution in the middle has two domains with the field
at its positive VEV in the chamber region $z<0$ and at its negative VEV for
$z>0$. We note that the latter solution displays the anti-symmetry
$\phi(r,-z)=-\phi(r,z)$, which shows that in the symmetry-broken phase the
field does not have to obey the symmetries of the chamber. Instead of a
solution having three domains along the $z$-axis, the solution displayed at
the bottom has been found. While it does display multiple local extrema along
either axis, none of them approaches any VEV too closely. This implies that
static solutions can exist without the need for the field to take on its VEV
anywhere inside the vacuum chamber. Despite using many more seeds for Newton’s
method, no further solutions have been obtained. However, we suspect that also
solutions without cylinder symmetry exist, since the symmetron field solutions
do not necessarily obey the symmetry of its environment. Finding such
solutions would require solving the full 3D differential equation, which is
beyond the scope of this article. We further suspect that solutions containing
more domains can exist only for larger parameter values of $\mu$, where the
field can have more curvature. This behavior, including the absence of any
solutions for too small values of $\mu$, has already been observed in [15, 16]
for the case of a 1-dimensional setup containing two parallel mirrors. The
energies of these field solutions and their interaction with the matter
density can be obtained by the Hamilton density

$\displaystyle\mathcal{H}$
$\displaystyle=\frac{1}{2}\left(\bm{\nabla}\phi\right)^{2}+V_{\text{eff}}(\phi)$
$\displaystyle=\frac{1}{2}\left\\{\left(\frac{\partial\phi}{\partial
r}\right)^{2}+\left(\frac{\partial\phi}{\partial
z}\right)^{2}\right\\}+V_{\text{eff}}(\phi)\>,$ (22)

where cylinder symmetry is assumed such that derivatives with respect to
$\varphi$ vanish. This results in the corresponding energy $E$ inside the
cylinder

$\displaystyle
E=\int\limits_{\text{cylinder}}d^{3}x\,\mathcal{H}(\vec{x})=2\pi\int_{0}^{R}dr\,r\int_{-L/2}^{L/2}dz\,\mathcal{H}(r,z)\>,$
(23)

where $R$ and $L$ refer to the radius and length of the cylinder,
respectively. The obtained energies are listed in Table 1.

Mode |  $E$ [eV]
---|---
0 | -2.097
1 | -1.710
2 | -0.06
Table 1: Energies $E$ are given corresponding to the three obtained solutions
with the mode number enumerating the solutions with increasing energy. Hereby,
the following numerical values have been used $\rho_{M}=1.16\times 10^{-5}$
MeV4, $\rho_{V}=7.08\times 10^{-17}$ MeV4, $\lambda_{S}=10^{-2}$, $M=10^{3}$
GeV, and $\mu=6.1\times 10^{-5}$ eV.

##  VI Comments on chameleon constraints

For the chameleon field, each value of $n$ in the potential corresponds to a
separate model. Restricting the analysis to parameters that have been studied
frequently, no constraints improving on already existing ones are obtained for
$n=-4$ (in this case $V_{C}(\phi)$ is replaced by
$V_{C}(\phi)=\lambda_{C}\phi^{4}$, with a dimensionlass coupling constant
$\lambda_{C}$) and $n=1$ with varying $\Lambda$. Also setting $\Lambda=2.4$
meV to the dark energy scale and varying $n$ for small values does not provide
more stringent constraints than those available. Therefore, no exclusion plots
are provided herein. For comparison to existing constraints use has been made
of Ref. [17].

##  VII Discussion and conclusion

Since very short ranges of scalar fields induce short-range fifth forces,
classical experiments eventually are no longer able to detect or probe scalar
fields in this regime. However, for neutron interferometry very short ranges
of scalar fields lead effectively to a constant potential shift inside a
vacuum chamber after a steep gradient resulting in a phase shift. The relative
phase shift due to the different potential shifts induced in the two chambers
can still be measured even when classical experiments would fail to detect any
force. For this reason, constraints were obtained for extremely large values
of $A_{2}$ and $\mu$ (which corresponds to small ranges) for the dilaton and
symmetron field. However, since the neutron has a finite extent, eventually
the screening of the neutron itself sets in and prohibits probing arbitrarily
small ranges.

While a split-crystal interferometer opens up the possibility to increase the
chamber length by a small factor, this does not yield a large improvement.
However, lowering the vacuum pressure and hence the vacuum density by a few
orders of magnitude (the current setup uses a vacuum pressure of $10^{-4}$
mbar or higher, which might eventually be lowered to $10^{-9}$ mbar) seems to
provide a more powerful way to probe screened scalar fields, since they
generically are less suppressed in regions of lower mass density. Finally, it
was found that for the symmetron inside cylindrical vacuum chambers in general
several distinct static field solutions exist.

##  VIII Acknowledgments

We thank Hartmut Abele and Tobias Jenke for fruitful discussions. This article
was supported by the Austrian Science Fund (FWF): P 34240-N, and is based upon
work from COST Action COSMIC WISPers CA21106, supported by COST (European
Cooperation in Science and Technology).

##  Appendix A Numerical solutions of the differential equations

In order to solve

$\displaystyle\Delta\phi(\bm{x})=V_{\text{eff},\phi}(\phi(\bm{x});\rho(\bm{x}))\>,$
(24)

Newton’s method has been applied on the function space. This method has
already previously been used to simulate chameleon fields [72]. Starting from
an initial guess function $\phi^{(0)}(\bm{x})$ for the field, the differential
equation is expanded to first order around that guess function. The solution
serves as an improved guess. This results in a sequence of functions
$\phi^{(0)}$, …$\phi^{(n)}$, … with $\phi^{(n)}$ defined as the solution of
the $linear$ differential equation

$\displaystyle\Delta\phi^{(n)}(\bm{x})=V_{\text{eff},\phi}(\phi^{(n-1)}(\bm{x});\rho(\bm{x}))$
$\displaystyle+V_{\text{eff},\phi\phi}(\phi^{(n-1)}(\bm{x});\rho(\bm{x}))\left(\phi^{(n)}-\phi^{(n-1)}\right).$
(25)

These linear differential equations were solved repeatedly with the finite
element method (see e.g. [73]) until

$\displaystyle\frac{||\phi^{(n)}-\phi^{(n-1)}||_{2}}{||\phi^{(n-1)}||_{2}}<\epsilon\>,$
(26)

for some small value of $\epsilon$ (which was at least set to $10^{-10}$, but
much smaller for extreme parameters). To find the one node solution, the
effective potential was first expanded around its minimum

$\displaystyle\left[\frac{1}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial}{\partial r}\right)+\frac{\partial^{2}}{\partial
z^{2}}\right]\phi(r,z)=$ $\displaystyle\qquad
V_{\text{eff},\phi\phi}(\phi_{\rho(r,z)})\left(\phi(r,z)-\phi_{\rho(r,z)}\right),$
(27)

which resulted in a linear differential equation that can be solved
explicitly. Then, the solution was used as an initial guess for Newton’s
method to solve the full non-linear differential equation. The multiple node
solutions for the symmetron field were found by trying a large number of
physically plausible seeds.

##  Appendix B Derivation of the phase shift formula

The derivation below largely follows [74] (see also [75]). Starting from the
free Schrödinger equation for the neutron

$\displaystyle\hat{H}_{0}\phi_{0}$
$\displaystyle=-\frac{1}{2m_{n}}\,\Delta\phi_{0}$
$\displaystyle=E_{0}\phi_{0}\>,$ (28)

where $\phi_{0}\propto e^{i\bm{k}\cdot\bm{x}}$, the ansatz

$\displaystyle\phi$ $\displaystyle=\phi_{0}\chi\>,$ $\displaystyle E$
$\displaystyle=E_{0}\>,$ (29)

is performed. The full Schrödinger equation reads

$\displaystyle\big{(}\hat{H}_{0}+U\big{)}\phi=E\phi\>.$ (30)

A straight forward calculation gives

$\displaystyle\phi_{0}\Delta\chi+2i\phi_{0}\bm{k}_{0}\cdot\bm{\nabla}\chi=2m_{n}U\phi_{0}\chi\>.$
(31)

The semi-classical limit is given by $|\bm{\nabla}\chi|\ll k_{0}|\chi|$, which
allows to neglect the first term. Defining $k_{0}:=|\bm{k}_{0}|=2\pi/\lambda$
and the length parameter $s$ along the direction of $\bm{k}_{0}$ gives
$\bm{k}_{0}\cdot\bm{\nabla}=k_{0}\,\frac{d}{ds}$ and

$\displaystyle\chi$ $\displaystyle=e^{-i\frac{m_{n}}{k_{0}}\int ds\,U}\>.$
(32)

The resultant phase shift is therefore given by

$\displaystyle\delta\varphi=-\frac{m_{n}}{k_{0}}\int ds\,U\>,$ (33)

and for neutrons propagating through a cylindrical cavity with length $L$
along the symmetry $z$-axis at constant radius $r_{c}$

$\displaystyle\delta\varphi=-\frac{m_{n}}{k_{0}}\int_{0}^{L}dz\,U(r_{c},z)\>.$
(34)

The relative phase shift between vacuum chamber and air chamber takes the form

$\displaystyle\delta\varphi_{\text{relative}}=-\frac{m_{n}}{k_{0}}\int_{0}^{L}dz\left(U_{\text{Vac}}(r_{c_{1}},z)-U_{\text{Air}}(r_{c_{2}},z)\right).$
(35)

Next, the condition for the validity of the semiclassical limit will be
derived. From

$\displaystyle\bm{\nabla}\chi=-i\frac{m_{n}}{k_{0}}\,U(r_{c},z)\,\chi\,\bm{e}_{z}\>,$
(36)

follows with $|\chi|=1$

$\displaystyle|\bm{\nabla}\chi|=\frac{m_{n}}{k_{0}}\,|U(r_{c},z)|\>.$ (37)

Hence, the condition for the validity of the semiclassical limit becomes

$\displaystyle\frac{|\bm{\nabla}\chi|}{k_{0}}=\frac{m_{n}}{k_{0}^{2}}\,|U(r_{c},z)|\ll
1\>.$ (38)

Finally, the derivation assumes that outside both chambers, i.e. in air, the
perturbation potential $U=0$. However, scalar fields in air typically take on
a non-vanishing value. For all parameter values where constraints can be set,
however, the ranges of the fields in air are extremely small compared to the
dimensions of the experimental setup. Therefore, the field is to a very good
approximation constant outside the chambers, having its air expectation value.
In order to comply with the assumption $U=0$ outside the chambers, we have to
define the potentials for the dilaton (D), symmetron (S) and chameleon (C)
perturbation potentials in the following way

$\displaystyle U_{D}(\bm{x})$
$\displaystyle:=\frac{A_{2}m_{n}}{2m^{2}_{\text{Pl}}}\left(\phi^{2}(\bm{x})-\phi^{2}_{\text{Air}}\right),$
(39) $\displaystyle U_{S}(\bm{x})$
$\displaystyle:=\frac{m_{n}}{2M^{2}}\left(\phi^{2}(\bm{x})-\phi^{2}_{\text{Air}}\right),$
(40) $\displaystyle U_{C}(\bm{x})$
$\displaystyle:=\frac{m_{n}}{M_{c}}\left(\phi(\bm{x})-\phi_{\text{Air}}\right).$
(41)

The potentials above are only valid in the limit where the neutron can be
considered as a test particle in a scalar field background rather than as a
source of the field, i.e. if the screening of the neutron can be neglected. In
order to take its screening into account, we multiply the above potentials
with a screening charge

$\displaystyle U_{X}(\bm{x})\rightarrow\mathfrak{Q}_{X}U_{X}(\bm{x})\>,$ (42)

with $X\in\\{D,S,C\\}$. For additional information on the screening charge of
neutrons see Appendix A of Ref. [15]. Therefore, the resulting conditions for
the validity of the semiclassical limit are given by

$\displaystyle\mathfrak{Q}_{D}\frac{A_{2}m_{n}^{2}}{2m^{2}_{\text{Pl}}k_{0}^{2}}\left(\phi^{2}(\bm{x})-\phi^{2}_{\text{Air}}\right)$
$\displaystyle\ll 1\ (\text{dilaton})\>,$ (43)
$\displaystyle\mathfrak{Q}_{S}\frac{m_{n}^{2}}{2M^{2}k_{0}^{2}}\left(\phi^{2}(\bm{x})-\phi^{2}_{\text{Air}}\right)$
$\displaystyle\ll 1\ (\text{symmetron})\>,$ (44)
$\displaystyle\mathfrak{Q}_{C}\frac{m_{n}^{2}}{M_{c}k_{0}^{2}}\left(\phi(\bm{x})-\phi_{\text{Air}}\right)$
$\displaystyle\ll 1\ (\text{chameleon})\>.$ (45)

These conditions are always fulfilled on those edges of the constrained areas,
which do not come from imposed cut-offs.

##  Appendix C Screening charge for the chameleon field

The following derivation and physical reasoning closely follows [25] and hence
only key steps in the derivation of the screening charge for the chamleon
field are provided. In order to compute the chameleon field of a spherical
source, one has to solve

$\displaystyle\frac{d^{2}\phi}{dr^{2}}+\frac{2}{r}\frac{d\phi}{dr}=V_{\text{eff},\phi}(\phi;\rho)\>.$
(46)

Since the field outside the sphere approaches its asymptotic value $\phi_{V}$
and inside the sphere $\phi_{S}$, the effective potential is expanded around
these values

$\displaystyle\frac{d^{2}\phi}{dr^{2}}+\frac{2}{r}\frac{d\phi}{dr}$
$\displaystyle=\mu^{2}_{S}\left(\phi-\phi_{S}\right),$ $\displaystyle r<R\>,$
$\displaystyle\frac{d^{2}\phi}{dr^{2}}+\frac{2}{r}\frac{d\phi}{dr}$
$\displaystyle=\mu^{2}_{V}\left(\phi-\phi_{V}\right),$ $\displaystyle r\geq
R\>,$ (47)

where $R$ refers to the radius of the sphere, $\mu_{S}$ to the chameleon mass
for the field taking its minimum value inside the sphere and $\mu_{V}$ to the
chameleon mass for the field taking its VEV. The final solution outside the
sphere is given by

$\displaystyle\phi(r)=\phi_{V}-\mathfrak{Q}_{C}\frac{\mu^{2}_{S}R^{3}}{3}\left(\phi_{V}-\phi_{S}\right)\frac{e^{-\mu_{V}(r-R)}}{r}\>,$
(48)

where $\mathfrak{Q}_{C}$ is the screening charge of the chameleon field
defined by

$\displaystyle\mathfrak{Q}_{C}:=\frac{3}{\mu^{2}_{S}R^{2}}\frac{1-\frac{1}{\mu_{S}R}\,\text{tanh}(\mu_{S}R)}{1+\frac{\mu_{V}}{\mu_{S}}\,\text{tanh}(\mu_{S}R)}\>,$
(49)

which describes the screening of the sphere. The limits

$\displaystyle\mathfrak{Q}_{C}\rightarrow\begin{dcases}0\>,&\text{for
``screened" bodies with }\mu_{S}R\rightarrow\infty\>,\\\ 1\>,&\text{for
``unscreened" bodies with }\mu_{S}R\rightarrow 0\>,\end{dcases}$

are obeyed.

## References

  * [1] S. Perlmutter _et al._ (Supernova Cosmology Project).  _Nature_ 391, 51 (1998).  doi:10.1038/34124.  eprint astro-ph/9712212.
  * [2] A. G. Riess _et al._ (Supernova Search Team).  _Astron. J._ 116, 1009 (1998).  doi:10.1086/300499.  eprint astro-ph/9805201.
  * [3] B. P. Schmidt _et al._ (Supernova Search Team).  _Astrophys. J._ 507, 46 (1998).  doi:10.1086/306308.  eprint astro-ph/9805200.
  * [4] J. Sola.  _J. Phys. Conf. Ser._ 453, 012015 (2013).  doi:10.1088/1742-6596/453/1/012015.  eprint 1306.1527.
  * [5] A. Joyce _et al._ _Phys. Rept._ 568, 1 (2015).  doi:10.1016/j.physrep.2014.12.002.  eprint 1407.0059.
  * [6] J. Khoury and A. Weltman.  _Phys. Rev. D_ 69, 044026 (2004).  doi:10.1103/PhysRevD.69.044026.  eprint astro-ph/0309411.
  * [7] J. Khoury (2013).  eprint 1312.2006.
  * [8] P. Brax, C. Burrage and A.-C. Davis.  _JCAP_ 01, 020 (2013).  doi:10.1088/1475-7516/2013/01/020.  eprint 1209.1293.
  * [9] P. Brax and P. Valageas.  _Phys. Rev. D_ 90(2), 023507 (2014).  doi:10.1103/PhysRevD.90.023507.  eprint 1403.5420.
  * [10] A. I. Vainshtein.  _Phys. Lett. B_ 39, 393 (1972).  doi:10.1016/0370-2693(72)90147-5.
  * [11] T. Damour and A. M. Polyakov.  _Nucl. Phys. B_ 423, 532 (1994).  doi:10.1016/0550-3213(94)90143-0.  eprint hep-th/9401069.
  * [12] P. Brax _et al._ _Phys. Rev. D_ 82, 063519 (2010).  doi:10.1103/PhysRevD.82.063519.  eprint 1005.3735.
  * [13] J. Sakstein.  _Astrophysical Tests of Modified Gravity_.  Ph.D. thesis, Cambridge U., DAMTP (2014).  doi:10.17863/CAM.16133.  eprint 1502.04503.
  * [14] G. Cronenberg _et al._ _Nature Phys._ 14(10), 1022 (2018).  doi:10.1038/s41567-018-0205-x.  eprint 1902.08775.
  * [15] P. Brax and M. Pitschmann.  _Phys. Rev. D_ 97(6), 064015 (2018).  doi:10.1103/PhysRevD.97.064015.  eprint 1712.09852.
  * [16] M. Pitschmann.  _Phys. Rev. D_ 103(8), 084013 (2021).  doi:10.1103/PhysRevD.103.084013.  [Erratum: Phys.Rev.D 106, 109902 (2022)], eprint 2012.12752.
  * [17] C. Burrage and J. Sakstein.  _Living Rev. Relativ._ 21(1), 1 (2018).
  * [18] M. Gasperini, F. Piazza and G. Veneziano.  _Phys. Rev. D_ 65(2), 023508 (2001).
  * [19] T. Damour, F. Piazza and G. Veneziano.  _Physical Review D_ 66(4), 046007 (2002).
  * [20] T. Damour, F. Piazza and G. Veneziano.  _Physical Review Letters_ 89(8), 081601 (2002).
  * [21] J. Müller _et al._ _J. Geodesy_ 93, 2195 (2019).  doi:10.1007/s00190-019-01296-0.
  * [22] H. Abele _et al._ _Physical Review D_ 81(6), 065019 (2010).
  * [23] T. Jenke _et al._ _Nature Physics_ 7(6), 468 (2011).
  * [24] R. I. P. Sedmik and M. Pitschmann.  _Universe_ 7(7), 234 (2021).  doi:10.3390/universe7070234.  eprint 2107.07645.
  * [25] P. Brax _et al._ _Eur. Phys. J. C_ 82(10), 934 (2022).  doi:10.1140/epjc/s10052-022-10905-w.  eprint 2203.12512.
  * [26] H. Fischer _et al._ (2023).  eprint 2307.00243.
  * [27] C. Käding, M. Pitschmann and C. Voith.  _Eur. Phys. J. C_ 83(8), 767 (2023).  doi:10.1140/epjc/s10052-023-11939-4.  eprint 2306.10896.
  * [28] C. Burrage _et al._ _Phys. Rev. D_ 100(7), 076003 (2019).  doi:10.1103/PhysRevD.100.076003.  eprint 1812.08760.
  * [29] C. Burrage _et al._ _J. Phys. Conf. Ser._ 1275(1), 012041 (2019).  doi:10.1088/1742-6596/1275/1/012041.  eprint 1902.09607.
  * [30] C. Käding and M. Pitschmann.  _Phys. Rev. D_ 107(1), 016005 (2023).  doi:10.1103/PhysRevD.107.016005.  eprint 2204.08829.
  * [31] K. Hinterbichler and J. Khoury.  _Physical review letters_ 104(23), 231301 (2010).
  * [32] K. Hinterbichler _et al._ _Physical Review D_ 84(10), 103521 (2011).
  * [33] C. Burrage, E. J. Copeland and P. Millington.  _Phys. Rev. D_ 95(6), 064050 (2017).  doi:10.1103/PhysRevD.95.064050.  [Erratum: Phys.Rev.D 95, 129902 (2017)], eprint 1610.07529.
  * [34] C. A. J. O’Hare and C. Burrage.  _Phys. Rev. D_ 98(6), 064019 (2018).  doi:10.1103/PhysRevD.98.064019.  eprint 1805.05226.
  * [35] C. Burrage _et al._ _Phys. Rev. D_ 99(4), 043539 (2019).  doi:10.1103/PhysRevD.99.043539.  eprint 1811.12301.
  * [36] C. Käding.  _Astronomy_ 2(2), 128 (2023).  doi:10.3390/astronomy2020009.  eprint 2304.05875.
  * [37] C. Burrage, E. J. Copeland and E. Hinds.  _Journal of Cosmology and Astroparticle Physics_ 2015(03), 042 (2015).
  * [38] C. Burrage and E. J. Copeland.  _Contemporary Physics_ 57(2), 164 (2016).
  * [39] C. Burrage _et al._ _Journal of Cosmology and Astroparticle Physics_ 2016(12), 041 (2016).
  * [40] P. Brax and A.-C. Davis.  _Physical Review D_ 94(10), 104069 (2016).
  * [41] A. Upadhye.  _Physical review letters_ 110(3), 031301 (2013).
  * [42] G. Cronenberg _et al._ _Nature Physics_ 14(10), 1022 (2018).
  * [43] T. Jenke _et al._ _Eur. Phys. J. ST_ 230(4), 1131 (2021).  doi:10.1140/epjs/s11734-021-00088-y.  eprint 2012.07472.
  * [44] P. Brax, A.-C. Davis and B. Elder.  _Physical Review D_ 107(4), 044008 (2023).
  * [45] J. Khoury and A. Weltman.  _Physical review letters_ 93(17), 171104 (2004).
  * [46] A. N. Ivanov _et al._ _Phys. Rev. D_ 94(8), 085005 (2016).  doi:10.1103/PhysRevD.94.085005.  eprint 1606.06867.
  * [47] D. Hartley _et al._ _arXiv preprint arxiv:1909.02272_ (2019).
  * [48] D. Hartley _et al._ _Phys. Rev. D_ 99(10), 105002 (2019).  doi:10.1103/PhysRevD.99.105002.  eprint 1811.06927.
  * [49] H. Lemmel _et al._ _Phys. Lett. B_ 743, 310 (2015).
  * [50] K. Li _et al._ (The INDEX Collaboration).  _Phys. Rev. D_ 93, 062001 (2016).  doi:10.1103/PhysRevD.93.062001.  URL https://link.aps.org/doi/10.1103/PhysRevD.93.062001.
  * [51] H. Lemmel _et al._ _Journal of Applied Crystallography_ 55(4) (2022).
  * [52] H. Abele _et al._ _In_ _11th Patras Workshop on Axions, WIMPs and WISPs_ , 124–129 (2015).  doi:10.3204/DESY-PROC-2015-02/abele˙hartmut.  eprint 1510.03063.
  * [53] T. Jenke _et al._ _EPJ Web Conf._ 219, 05003 (2019).  doi:10.1051/epjconf/201921905003.
  * [54] M. Pitschmann and H. Abele.  _arXiv preprint arxiv:1912.12236_ (2019).
  * [55] R. I. P. Sedmik _et al._ _EPJ Web Conf._ 219, 05004 (2019).  doi:10.1051/epjconf/201921905004.  eprint 1908.09723.
  * [56] S. Sponar _et al._ _Nature Reviews Physics_ 3, 309 (2021).  doi:10.1038/s42254-021-00298-2.  URL https://doi.org/10.1038/s42254-021-00298-2.
  * [57] A. N. Ivanov, M. Wellenzohn and H. Abele.  _Phys. Lett. B_ 822, 136640 (2021).  doi:10.1016/j.physletb.2021.136640.  eprint 2109.09982.
  * [58] M. Suda _et al._ _Z. Naturforsch. A_ 77(9), 875 (2022).  doi:10.1515/zna-2022-0050.  eprint 2111.02769.
  * [59] N. Muto _et al._ _JINST_ 17(07), P07014 (2022).  doi:10.1088/1748-0221/17/07/P07014.  eprint 2201.04346.
  * [60] C. Käding, M. Pitschmann and H. Abele.  _Z. Naturforsch. A_ 78(7), 651 (2023).  doi:10.1515/zna-2023-0045.  eprint 2302.11429.
  * [61] P. Brax, C. Burrage and A.-C. Davis.  _Int. J. Mod. Phys. D_ 27(15), 1848009 (2018).  doi:10.1142/S0218271818480097.
  * [62] H. Rauch, W. Treimer and U. Bonse.  _Phys. Lett._ 47A(5), 369 (1974).  doi:DOI:10.1016/0375-9601(74)90132-7.
  * [63] H. Rauch and S. A. Werner.  _Neutron Interferometry_.  Clarendon Press, Oxford (2000).
  * [64] H. Rauch _et al._ _Phys. Lett._ 54A, 425 (1975).
  * [65] H. Rauch _et al._ _Physical Review A_ 53(2), 902 (1996).
  * [66] K. H. J. Khoury.  _arXiv preprint arxiv:1001.4525_ .
  * [67] A. Upadhye.  _Phys. Rev. Lett._ 110(3), 031301 (2013).  doi:10.1103/PhysRevLett.110.031301.  eprint 1210.7804.
  * [68] T. Vachaspati.  _Kinks and domain walls: An introduction to classical and quantum solitons_.  Cambridge University Press (2007).
  * [69] Ø. Christiansen _et al._ _Journal of Cosmology and Astroparticle Physics_ 2023(05), 009 (2023).
  * [70] C. Llinares and P. Brax.  _Physical Review Letters_ 122(9), 091102 (2019).
  * [71] K. Clements _et al._ _arXiv preprint arXiv:2308.01179_ (2023).
  * [72] C. Briddon _et al._ _. Cosmol. Astropart. Phys._ 2021(12), 043 (2021).
  * [73] H. P. Langtangen and K.-A. Mardal.  _Introduction to numerical methods for variational problems_ , vol. 21.  Springer Nature (2019).
  * [74] D. M. Greenberger and A. Overhauser.  _Reviews of Modern Physics_ 51(1), 43 (1979).
  * [75] M. Pitschmann.  _The High Precision Frontier: Search for New Physics with “Tabletop Experiments” & Beyond_.  habilitation, TU Wien (2023).