# Constructing Nano-Object Quantum Superpositions with a Stern-Gerlach
Interferometer

Ryan J. Marshman  Department of Physics and Astronomy, University College
London, Gower Street, WC1E 6BT London, United Kingdom.  Anupam Mazumdar  Van
Swinderen Institute, University of Groningen, 9747 AG Groningen, The
Netherlands.  Ron Folman  Department of Physics, Ben-Gurion University of the
Negev, Be’er Sheva 84105, Israel.  Sougato Bose  Department of Physics and
Astronomy, University College London, Gower Street, WC1E 6BT London, United
Kingdom.

###### Abstract

Probing quantum mechanics, quantum aspects of general relativity, along with
the sensing and the constraining of classical gravity can all be enabled by
unprecedented spatial sizes of superpositions of massive objects. In this
paper, we show that there is a feasible setup sourced by realisable magnetic
field gradients ${\cal O}(10-100)$ Tm-1 to construct a large spatial
superposition of ${\cal O}(10^{-4}-10^{-8})$ m for ${\cal
O}(10^{-17}-10^{-14})$ kg masses within a time of up to $0.1-10$ seconds. The
scale of superpositions is unrestricted as long as quantum coherence can be
maintained for a required amount of time.

##  I Introduction

At a microscopic level, three of the known forces of nature, electromagnetic
(EM), weak and strong, obey the principles of local quantum field theory (QFT)
[1]. However, there is no experimental proof yet of how the gravitational
interaction is being mediated. Recently, a tabletop experiment has been
suggested to explore the quantum origin of gravity in a tabletop experiment
[2, 3, 4]. The protocol crucially relies on the interaction of quantum matter
with quantum gravity, leading to the generation of entanglement between spins
embedded in two non-relativistic test masses [3]. The spin entanglement
witness will prove the graviton’s quantum-ness as a mediator of the force
between the two masses (quantumness of the linearized metric fluctuations
around a Minkowski background), and will further test the nature of the
gravitational interaction at microscopic distances [3, 5]. However, there
remains a demanding requirement: quantum spatial superpositions of distinct
localized states of neutral mesoscopic masses $m\sim 10^{-15}$ kg over spatial
separations of $\Delta x\sim 10$ $\mu$m [6], far beyond the scales achieved to
date (e.g., macromolecules $m\sim 10^{-22}$ kg over $\Delta x\sim 0.25$
$\mu$m, or atoms $m\sim 10^{-25}$ kg over $\Delta x\sim 0.5$ m) [7, 8]. Such
superpositions have also been shown to be of practical value in sensing of
weak forces, curvature [9], frame-dragging, and even a tabletop detection of
low-frequency gravitational waves [10].

Beyond the usage in sensing quantum and classical gravity, upgrading the mass
$m$ and the superposition $\Delta x$, naturally stretches the boundaries of
the validity of QM, which in itself is a worthy goal [11, 8]. However, there
is a gap in the literature at the moment as far as a realistic scheme for
achieving superposition sizes $\Delta x>1$ $\mu$m. While there are quite a few
schemes for a lower $\Delta x$, or $m$, investigated at various levels of
detail [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28,
29, 30, 31, 32], and while these may suffice to falsify various purported
modifications of QM [33, 34, 35], the only two predictions stemming from the
straightforward application of QM (local QFT to be more precise [3]) which
would give us non-trivial information (e.g. sensing the quantum nature of
gravity [2] and sensing extremely weak classical gravity [10]), necessarily
seem to require $\Delta x\gtrsim 1$ $\mu$m.

Rudimentary arguments on how to achieve such superpositions using a Stern-
Gerlach Interferometer (SGI) in high magnetic field gradients were presented
in [2], building on earlier ideas for smaller $m,\Delta x$ superpositions [20,
36]. A feasibility study building on atomic experiments was recently
performed, showing that SGI for massive objects is indeed possible [31].
However, beyond a simple scaling of mass, little work has been done [29] on
exploring what new issues will arise when these techniques are pushed beyond
the atomic level.

In this paper, we go further in several respects by incorporating some crucial
aspects missed in previous treatments while making it simultaneously less
demanding in certain aspects, and confirm the true potential for the spatial
splitting of massive objects (namely nanocrystals) achieved by SGI. The
proposal outlined here requires only moderate magnetic gradients,
significantly less than what was originally considered necessary[2] as shown
in Figure 3. We will be describing a 1D longitudinal interferometer, which
avoids the problems noted for 2D interferometers [37]. We take into account
gradients in other directions, as demanded by Maxwell equations. We take into
account the effects of the induced diamagnetism within the interferometric
mass as while this is not of concern with atomic Stern-Gerlach interferometry,
it becomes a dominant effect once larger masses are considered. We consider a
magnetic source which enables constant gradients over a large volume. In
contrast to the previous works [19, 20, 29], we have a scheme without a low
magnetic field in any region, so that to avoid the historically well known
phenomenon of Majorana spin flips [38, 39]. In this way we have combined well
studied problems in particle trapping (Majorana spin flips), and practical
realizations of Stern-Gerlach effect to show that a full-loop interferometer
[31], where wavepackets are spatially spilt in a spin dependent manner and
then brought back to overlap both in position and momentum to complete the
interferometry, is still possible with nanocrystals.

The above results are simultaneously achieved by a significant modification of
standard SGI by changing the initial conditions and incorporating a gradient-
free spatial region in which the diamagnetic force does not act, and yet, the
superposition continues to grow in view of a momentum difference. However, we
show what appears to be an unavoidable price to pay, that is a linear growth
of a relative phase between the superposed components in time. These features
allow us to present in detail a feasible interferometric accelerator for
micro-scale masses, and enable the creation of a massive Schrödinger’s cat.

Previous experimental configurations consider a magnetic field which
originates in a single current-carrying wire or a permanent magnet, whereby
the field goes as $\vec{B}\propto 1/r$, where $r$ is the distance from the
source. The magnetic field can then be linearised in a small region within
which it has an approximately constant gradient. However, a large splitting
will require long evolution times, and as the distance quickly increases
between the particle and the wire, significant magnetic gradients are no
longer available. We will therefore consider a configuration which enables a
constant gradient over a large region (e.g. quadrupole field from coils in an
anti-Helmholtz configuration). A typical resultant trajectory of such a
potential is shown in Figure 1.

Figure 1: Example interferometer path trajectories through the varying
magnetic fields used. The magnetic field transition time here is accentuated
for readability. Note that the internal spin states are reversed at time
$t\approx 0.9$ s to ensure the interferometer is closed.

##  II Model

We assume a host nano-crystal with a single spin embedded in it, and
concentrate, in this work, only on the translational motion of this crystal
under a magnetic field gradient. To explicitly model the spin and the host
crystal material, we will assume a nitrogen-vacancy (NV) centre spin in a
diamond nano/micro-crystal, although this analysis holds for all materials
with similar diamagnetic susceptibilities and generic electronic spins. We
consider the mass to be in a free fall along the $z$-axis. We therefore do not
treat its motion along the $z$ axis explicitly as here the $z$ motion commutes
with the motion along the other two axes. Moreover, $z$ motion would be absent
in a drop-tower experiment, for example, which may be necessary in such
experiments because of their long durations of $\sim 1$ s, and to mitigate
gravitational jitters [9]. The Hamiltonian of the system is given by [40, 29]

$H=\frac{\hat{P}^{2}_{x}+\hat{P}^{2}_{y}+\hat{P}^{2}_{z}}{2m}+\hbar
D\hat{S}^{2}_{x}-\frac{\chi_{m}m}{2\mu_{0}}\vec{B}^{2}-\mu_{B}\hat{s}\cdot\vec{B}\,,$
(1)

where $m$ is the mass of the diamond, $\chi_{m}=-6.2\times 10^{-9}$ m3kg-1 is
the mass magnetic susceptibility, $\mu_{B}$ is the Bohr magneton, $\hat{s}$ is
the spin operator, and $D=\left(2\pi\right)\times 2.8$ GHz is the NV zero-
field splitting. Note that in writing the above Hamiltonian, we are only
concentrating on the translational motion – we are assuming that we can
engineer a situation so that the torque and the rotational effects of the mass
is negligible/decoupled from translation (for possible mechanisms to cool
rotation see [41, 42, 43, 44]). The external magnetic field $\vec{B}$ varies
in the $x-y$ plane, and as we will see, with our choice of its profile,
couples predominantly to the $x$ motion via the $x$-component of the spin.
Under this setting, once the internal spin of the NV centre is initialised in
the state
$\frac{1}{\sqrt{2}}\left(\left|1\right\rangle+\left|-1\right\rangle\right)$,
we will realize a 1-dimensional SGI which, as time progresses, achieves first
a maximum wavefunction splitting $\Delta x$ in the $x$ direction, and then,
via appropriate spin flips/magnetic gradient changes, its subsequent
recombination to complete the interferometry. The second term of the
Hamiltonian drops out as a common phase and we will omit this term from our
analysis. This also means that our results will be completely general, and
hold for other (say spin$-1/2$) dopants in generic crystals.

There are known predicaments for generating superpositions by using a spin
coupled to a magnetic field: (1) The electronic spin will in general
experience Majorana spin-flips whenever the magnetic field magnitude becomes
small, thus no longer remaining in a well-defined eigenstate as required for
the coherent manipulation of the masses (this is particularly important for
our results to be inclusive of generic spin dopant atoms in a nanocrystal; it
has been perceived as an important problem for several atomic species [39]),
and (2) the effect of off-axis magnetic field gradients, which are generally
not considered [29], however, this _must_ be taken into account to satisfy
Maxwell’s equations.

First, we will consider the constraints on the magnetic field to address the
above problems. Noting that $\vec{\nabla}\cdot\vec{B}=0$ and
$\vec{\nabla}\times\vec{B}=0$, as we keep the mass away from the source of the
magnetic field. For simplicity, we assume the following profile in the $(x,y)$
plane: $\vec{B}=B_{x}(x,y)\hat{x}+B_{y}(x,y)\hat{y}$. We take the $x$ axis as
the desired superposition direction, and we require the magnetic field to be
linear along $x$ direction, with a constant magnetic field gradient along the
$x$ direction to take some constant value $\partial B$. To ensure that we can
take the x-axis as our quantisation axis, and that Majorana spin flips are
avoided, we require that
$\left|B_{x}(x,y)\right|\gg\left|B_{y}(x,y)\right|\forall x,y$ in the vicinity
of the controlled trajectories. The simplest general form for the magnetic
field which will satisfy all the above conditions can be given by:

$\vec{B}(x,y)=\left(B_{x}(0,0)-\partial Bx\right)\hat{x}+\partial
By\hat{y}\,,$ (2)

for any fixed value of the magnetic field at the origin $B_{x}(0,0)$ and
magnetic field gradient $\partial B$. Note that for positive $B_{x}(0,0)$ and
$\partial B$ the zero point of the magnetic field is always on the positive
$x$ axis. By ensuring that $y\approx 0$ during the entire interferometry, we
find a suitable definition for the zero magnetic field region which must be
avoided of

$x\notin\left[\frac{B_{x}(0,0)-\varepsilon}{\partial
B},\frac{B_{x}(0,0)+\varepsilon}{\partial B}\right]\,,$ (3)

where $\varepsilon$ is the minimum allowable magnetic field in the $x$
direction. For a sufficiently large $\varepsilon$, we can ascertain that the
spin states will always be approximately aligned along the quantisation ($x$)
axis and Majorana spin flips are avoided.

The spin state in the $y$ and $z$ basis will experience a rapid Larmor
precession, with a frequency set by $\omega_{L}=-\frac{ge}{2m_{e}}|B(x,y)|$,
where $g\approx 2$ is here the Landè g factor, $e$ is the electron charge and
$m_{e}$ is the electron mass. If we desire a minimum Larmor frequency of
$\omega_{L}^{min}$, we can define a minimum allowable magnetic field
magnitude:

$\varepsilon\sim-\frac{2m_{e}\omega_{L}^{min}}{ge}.$ (4)

Therefore, the particle must not enter the region given by Eq.(3). To achieve
this, we will use three linear magnetic field profiles sequentially in time,
which are specific solutions of the form of Eq.(2):

$\displaystyle\vec{B}\left(x,y\right)=$ $\displaystyle\left(B_{0}-\eta
x\right)\hat{x}+\eta y\hat{y}$ (5) $\displaystyle\vec{B}\left(x,y\right)=$
$\displaystyle B_{1}\hat{x}$ (6) $\displaystyle\vec{B}\left(x,y\right)=$
$\displaystyle-\left(B_{0}-\eta x\right)\hat{x}-\eta y\hat{y}$ (7)

with $B_{0},B_{1},\eta>0$. We then require $B_{1}\gg\varepsilon$, and the
timing of the switching between the magnetic fields will be done to ensure
that the particle never experiences a small magnetic field. That is, as the
particle approaches the disallowed region given by Eq. 3 (for
$B_{x}(0,0)=B_{0}$ and $\partial B=\eta$), the magnetic field is mapped to
that given by Eq. 6 (where now $B_{x}(0,0)=B_{1}$ and $\partial B=0$, thus the
region specified by Eq. 3 does not occur), and as it leaves that region, the
magnetic field is smoothly mapped to Eq. 7 (where $B_{x}(0,0)=-B_{0}$ and
$\partial B=-\eta$). The switching function will be modelled by
$\textrm{Sw}(t,t_{on},t_{off})=0.5\left(\tanh\left[\delta\left(t-t_{on}\right)\right]+1\right)\times
0.5\left(\tanh\left[\delta\left(t_{off}-t\right)\right]+1\right)$ where
$\delta$ is the switching frequency parameter, which we ensure to be
sufficiently slow such that the magnetic field change is always adiabatic.
This is to ensure that it does not complicate the spin dynamics. We consider
$\delta=10^{3}$ Hz, which is well below what is required to maintain the
adiabaticity conditions $\dot{\omega}_{L}/\omega_{L}^{2}\ll 1$ and
$\delta\ll\omega_{L}^{min}$. With this, since the magnetic field is along a
fixed direction, it effectively freezes the spin direction (with the other
spin components experiencing high-frequency Lamor precession) ensuring the
interferometer is effectively one dimensional.

The motion of the particle can be separated into different stages, largely
depending on the form of $\vec{B}$:

  1. 1.

$t<\tau_{1}$, $\vec{B}$ is given by Eq.5.

  2. 2.

$\tau_{1}\leq t<\tau_{2}$, $\vec{B}$ is switching adiabatically from Eq.5 to
Eq.6, via the switching function $\textrm{Sw}(t,t_{on},t_{off})$.

  3. 3.

$\tau_{2}\leq t<\tau_{3}$, $\vec{B}$ is given by Eq.6.

  4. 4.

$\tau_{3}\leq t<\tau_{4}$, $\vec{B}$ is switching adiabatically from Eq.6 to
Eq.7 via the switching function $\textrm{Sw}(t,t_{on},t_{off})$.

  5. 5.

$\tau_{4}\leq t$, $\vec{B}$ is given by Eq.7.

  6. 6.

$t=\tau_{5}$, when the spin states are reversed to close the superposition
while $\vec{B}$ is still given by Eq.7.

  7. 7.

$t=\tau_{6}$, when the two wavefunctions are brought to overlap in the
position and the momentum basis.

Away from the nearly zero-field region, where the magnetic field profile is
given by Eq.5 or 7, we can write the potential energy in a compact form as:

$U_{\pm}(\tilde{x}^{\prime})=-\frac{\chi_{m}m}{2\mu_{0}}\eta^{2}\tilde{x}^{\prime
2}+\mu_{\pm}\frac{\mu_{0}}{\chi_{m}m}$ (8)

where we have taken $y\approx 0$, while the second term is a constant energy,
and $\tilde{x}^{\prime}=x-C(s_{x})$, where

$C(s_{x}=\pm
1)={B_{0}}/{\eta}\pm\frac{ge\hbar}{2m_{e}}\frac{\mu_{0}}{\chi_{m}m\tilde{\eta}(t)}\,,$
(9)

where $\tilde{\eta}(t)=\eta$ when $t<\tau_{1}$ and $\tilde{\eta}(t)=-\eta$
when $t>\tau_{4}$. Thus the object sees a harmonic potential created by the
diamagnetic interaction whose centre is displaced in the $x$ direction by the
spin-magnetic field gradient interaction. It is clear from the above that the
object will roll in different potential wells corresponding to its spin state
and thereby develop a momentum difference as they approach the nearly zero
field region. It is in this region that we switch the gradient (the harmonic
potential) off, and let the object evolve in a magnetic field given by Eq.(6)
so that a continually increasing spatial splitting can develop due to the
momentum difference of the spin components, while Majorana spin flips are
avoided.

Figure 2: Figure (a) shows both the paths taken through the interferometer and
the corresponding velocity while figure (b) shows the superposition size with
time. These plots are the result of the combined numerical and analytical
analysis. The vertical grey lines represent times $\tau_{1}\approx\tau_{2}$,
$\tau_{3}\approx\tau_{4}$, $\tau_{5}$ and $\tau_{6}$. Both figures are for
$B_{0}=10^{-2}$ T, $\eta=46.483$ Tm-1 and $m=10^{-17}$ kg.

When the masses are in a coherent state of the harmonic potential, the
wavefunction will not spread, while when in the constant magnetic field (free
motion) it will hardly spread due to the largeness of the mass. Thus, it
suffices to model the centre of mass as following a classical trajectory for
each interferometric path. Thermal fluctuations in the initial state do not
limit the coherence in the final state as they factor out of the motion, see
the discussions in [12, 20]. To a good approximation, the two interferometric
paths can be modelled by $x(t)=A\cos\left(\omega t+\phi\right)+C(s_{x})$,
where the frequency of the diamagnetic trap is given by
$\omega=(-{\chi_{m}}/{\mu_{0}})^{1/2}\eta$, while note that $\chi_{m}<0$ as
the mass is diamagnetic. This leaves the amplitude $A$ and the phase $\phi$ to
be determined by requiring that the position and the momentum of each arm of
the interferometer are continuous throughout the trajectory. For the ease of
computation, we have solved the complete trajectories using an appropriate
combination of analytic and numerical solutions. For periods of evolution by a
time independent Hamiltonian, analytic solutions are used. However to allow
for a more realistic magnetic field switching to be modelled, numerics are
used, specifically during $\tau_{1}\leq t<\tau_{2}$ and $\tau_{3}\leq
t<\tau_{4}$. This is discussed in more detail in Appendix A. The resulting
motion can be seen in Fig.2. Note that at $\tau_{5}\approx 0.9$ s (marked by
the vertical lines in Fig.2), the spin reversal takes place by firing a rapid
microwave pulse to alter the internal spin state. We can find numerically the
final time $\tau_{6}$ to be

$\tau_{6}\approx 59\times({1{\rm~{}Tm^{-1}}}/{\eta})~{}{\rm~{}sec}\,,$ (10)

to ensure that the relative positions of the two paths $\Delta
x(t=\tau_{6})\approx 0$ m, and the relative velocity $\Delta
v(t=\tau_{6})\approx 0$ ms-1 (completion of interferometry).

Figure 3: Superposition size vs. time. This shows the maximal superposition
size achievable, as determined by the calculated paths. Note that for masses
$m=10^{-17}$ kg and heavier these results are well approximated by Eq. 11.

This results in the superposition sizes as seen in Fig.3, and the trajectories
through the interferometer can be seen as in Fig.2. Furthermore, the maximum
superposition which occurs between $\tau_{4}<t<\tau_{5}$ can be determined
numerically:

$\Delta x_{\textrm{max}}(m,\tau_{6})\approx\left(\frac{1.6\times
10^{-16}{\rm~{}Kg}}{m}\right)\left(\frac{\tau_{6}}{{\rm 1~{}sec}}\right)\times
10^{-6}\textrm{~{}m}$ (11)

for $B_{0}=10^{-2}$ T, $B_{1}=100\varepsilon$ and $m\gtrsim 10^{-17}$ kg. The
latter condition on the mass arises because smaller masses are subject to much
smaller induced diamagnetic potential. For a lighter diamond, if we do not
modify the times $\tau_{i}\left(\eta\right)$ and $B_{1}$, the masses would
inevitably move through the zero-field region of the magnetic field, and
therefore demands a different magnetic field setup, see [45, 31]. Therefore,
our analysis holds true for $m\geq 10^{-17}$ kg for $\chi_{m}\approx-6.2\times
10^{-9}$ kgm-3. Using Eq.11 we can estimate that achieving a superposition
size of $20$ $\mu$m with a $10^{-17}$ kg mass requires a total time of
$\tau_{6}\approx 1.25$ s, which corresponds to the moderate magnetic field
gradient $\eta\sim 46.8$ Tm-1, again using $B_{0}=10^{-2}$ T and
$B_{1}=100\varepsilon$, which can be achieved in a laboratory [45, 46].

##  III Spin Phase Evolution

Figure 4: Phase difference magnitude scaling with total runtime for
$m=10^{-15}$ kg, $m=10^{-16}$ kg, and $m=10^{-17}$ kg. This also shows the
phase difference as predicted by Equation 16.

To determine the output signal expected from such an interferometer it is
necessary to consider the phase evolution difference between the two paths
taken through the interferometer. The final output signal will be the result
of a combination of the phase accumulated by the phase due to the path it
takes through the interferometer, and external forces acting on the masses,
such as gravitational or electromagnetic interactions which differ across the
paths taken through the interferometer. For the moment however, we will
neglect these outside sources due to highly implementation dependent nature of
these. Given the entire Hamiltonian is used to determine the trajectories, it
is sufficient to calculate the action for a _free_ spin travelling along a
fixed trajectory given by the two trajectories through the interferometer.
This will then yield a phase, $\theta$ accumulated along a path $\gamma$ of

$\theta=\frac{1}{\hbar}\int_{\gamma}\frac{p(t)^{2}}{2m}\textrm{d}t$ (12)

where $p(t)$ is the time dependent momentum of the particle. To evaluate this
we can approximate the motion as being in three distinct segments:
$\tau_{0}=0\leq t\leq\tau_{1}$, before the particle has reached the zero field
region; $\tau_{2}<t\leq\tau_{3}$, after the particle has reached the zero
field region but before the spin has been reversed; and
$\tau_{3}<t\leq\tau_{4}$, the remainder of the evolution, ending when the two
wave packets are brought to overlap. Note that the time region
$t\in\left[\tau_{1},\tau_{2}\right]$ is extremely brief and does not
significantly change the end result. In each of these sections the general
particle momentum can be written as

$p^{+/-}(t)=-m\omega A_{i}^{+/-}\sin\left(\omega t+\phi_{i}^{+/-}\right)$ (13)

for the $+$ and $-$ arm of the interferometer and the value of $A^{+/-}$ and
$\phi^{+/-}$ differs from one segment to the next. Thus the phase difference
accumulated through any given stage of the interferometer will be

$\displaystyle\Delta\theta=$
$\displaystyle\sum_{i\in\left\\{1,3,4\right\\}}\Bigg{[}\frac{m}{2\hbar}\int_{\tau_{i-1}}^{\tau_{i}}\left(\omega
A_{i}^{+}\sin\left(\omega t+\phi_{i}^{+}\right)\right)^{2}$
$\displaystyle-\frac{m}{2\hbar}\int_{\tau_{i}}^{\tau_{i+1}}\left(\omega
A_{i}^{-}\sin\left(\omega t+\phi_{i}^{-}\right)\right)^{2}\Bigg{]}$
$\displaystyle=$
$\displaystyle\frac{m\omega}{8\hbar}\left[2\omega\tau_{1}-\sin\left(2\omega\tau_{1}\right)\right]\left[\left(A_{1}^{+}\right)^{2}-\left(A_{1}^{-}\right)^{2}\right]$
$\displaystyle+\sum_{i=3}^{4}\frac{m}{8\hbar}\Bigg{[}2\omega^{2}\left(T_{i}-T_{i-1}\right)\left(\left(A^{+}\right)^{2}-\left(A^{-}\right)^{2}\right)$
$\displaystyle+\sin\left(2\phi^{+}+2T_{i-1}\omega\right)-\sin\left(2\phi^{-}+2T_{i-1}\omega\right)$
$\displaystyle-\sin\left(2\phi^{+}+2T_{i}\omega\right)+\sin\left(2\phi^{-}+2T_{i}\omega\right)\Bigg{]}$
(14)

Now it is suitably insightful to consider in detail the phase difference
accrued when $t\in\left[\tau_{0},\tau_{1}\right]$, noting that
$\sin\left(2\omega\tau_{1}\right)\approx 0$:

$\displaystyle\Delta\theta\left(t=\tau_{1}\right)\approx$
$\displaystyle\frac{m}{8\hbar}\Bigg{[}2\omega^{2}\tau_{1}\left(\left(A_{1}^{+}\right)^{2}-\left(A_{1}^{-}\right)^{2}\right)$
$\displaystyle-\sin\left(2\tau_{1}\omega\right)+\sin\left(2\tau_{1}\omega\right)\Bigg{]}$
$\displaystyle=$
$\displaystyle\frac{m}{8\hbar}2\omega^{2}\tau_{1}\left[\left(\frac{B_{0}}{\eta}+\alpha\right)^{2}-\left(\frac{B_{0}}{\eta}-\alpha\right)^{2}\right]$
$\displaystyle=$ $\displaystyle\frac{geB_{0}}{2m_{e}}\tau_{1}$

Given that in each time segment considered here, the particles are simply
acting as harmonic oscillators, flipping between two different harmonic wells,
the entire path phase difference can be approximated as

$\Delta\theta\approx\frac{geB_{0}}{2m_{e}}\left(\frac{\tau_{1}}{\tau_{4}}-\frac{\tau_{3}-\tau_{1}}{\tau_{4}}+\frac{\tau_{4}-\tau_{3}}{\tau_{4}}\right)\tau_{4}$
(16)

where the extra terms in equation 14 are neglected and the negative time
scaling occurs when the internal spin direction is reversed relative to the
external field. Note that this occurs by reversing the field magnitude, not
the particle spin state. This entire term
$\left(\frac{\tau_{1}}{\tau_{4}}-\frac{\tau_{3}-\tau_{1}}{\tau_{4}}+\frac{\tau_{4}-\tau_{3}}{\tau_{4}}\right)\approx
0.5$ can be viewed as a comparison between the time spent with the masses
being accelerated away from one another and time spent being attracted towards
one another. This would typically be equal and as such this entire phase would
not usually appear. It is only due to the asymmetry in this interferometer set
up that it occurs. Thus, for the full trajectories, the path phase difference
will be $\Delta\theta\propto\tau_{4}$. This can be seen clearly in Figure 4,
which shows the phase scaling linearly with time and almost independently of
mass. This also shows the approximate phase difference given by Equation 16
and how it compares to the exact values. This also allows for the required
stability to be estimated, using equation 16:
$\Delta\theta\approx\frac{ge}{4m_{e}}B_{0}\tau_{4}\sim B_{0}\tau_{4}\times
10^{11}\textrm{ T}^{-1}\textrm{s}^{-1}$; to keep the final phase uncertainty
$\delta\left(\Delta\theta\right)<1$ will require

$\displaystyle 10^{-11}\textrm{ Ts}>$ $\displaystyle B_{0}\times\delta t$ (17)
$\displaystyle 10^{-11}\textrm{ Ts}>$ $\displaystyle\delta B_{0}\times t$ (18)

which places an strict requirement on both timing certainty, $\delta t$; and
bias field stability, $\delta B_{0}$. This will serve as a further challenge
that must be met before such an experiment can be fully realised.

##  IV Discussion

In this paper, we have provided a simple mechanism to accelerate heavy neutral
masses, with embedded spin, so as to create large spatial superpositions.
Using Eq.(11) we estimate that one can create, for example, the spatial
superposition of $\sim 20{\rm~{}\mu m}$ for masses, as heavy as $10^{-17}$ kg.
Our simple scheme fills the gap required for a realistic “wavefunction
splitting” of large masses to achieve a large spatial superposition, and
herewith opens up new vistas for testing and probing both the classical and
quantum nature of gravity while also giving access to unprecedented sensing
opportunities.

The time duration of coherence, and hence the experiment, is the only limiting
factor, but (a) spin coherence times are perpetually rising (approaching $1$s
[47, 48], even 30 s [49, 50])– adapting these to nanocrystals remains an open
problem, but is not fundamentally restricted [51], (b) spatial coherence times
can be made $\sim 100$ s [2, 6, 10, 52], by challenging, but achievable
pressures [53], temperatures [54], distances from other sources and
fluctuations. For example, a decoherence rate below $0.1$ Hz is achievable for
diamond spheres of masses $\sim 10^{-14}$ kg. This is expected [55] for
internal temperatures of $0.15$ K, an environmental temperature of $1$ K and a
environmental gas number density of $10^{8}$ m-3.

Effort will also have to be taken to initialise the particle with the spin
direction aligned with the external applied magnetic field to minimise
spurious torques arising from the spin-magnetic field coupling. This can be
accomplished already before being released from a trap, for example, by using
anisotropically shaped nanoparticles, which can be aligned with any given
direction in space by using linearly polarized lasers or electric fields [56]
or magnetic fields [57]. There will still be harmonic motion about this
orientation axis (called librations), which generally have much higher
frequencies [58] than typical centre of mass trapping frequencies in tweezers.
As the latter has already been cooled to the ground state by feedback cooling
[59, 60, 61, 62], one would expect the former (librations) to be possible to
be cooled to the ground state also. In fact, these librations have been
substantially cooled [63]. Once it is cooled close to the ground state, the
spin direction will be effectively exactly aligned to the external magnetic
field. If the magnetic field is maintained throughout the interference in the
same direction then the spin, and hence the nanoparticle, receives no torque
from it. We also note that it has already been shown that the internal degrees
of freedom (phonons) do not pose a practical problem [64]. We have shown that
there are also requirements on the magnetic field fluctuations and timing
certainty to ensure a stable interference signal, both at achievable levels,
however [65, 66, 67, 68]. To conclude, we have found an explicit scheme to
create large spatial superpositions, with the size increasing in proportion to
achievable coherence times, but otherwise not limited.

Acknowledgements: R. J. M. is supported by a EPSRC DTP departmental
studentship at UCL. AM’s research is funded by the Netherlands Organisation
for Science and Research (NWO) grant number 680-91-119. R. F. is supported by
the Israel Science Foundation. SB would like to acknowledge EPSRC Grant Nos.
EP/N031105/1 and EP/S000267/1.

## References

  * Weinberg [2005] S. Weinberg,  _The Quantum theory of fields. Vol. 1: Foundations_ (Cambridge University Press, 2005).
  * Bose  _et al._ [2017] S. Bose, A. Mazumdar, G. W. Morley, H. Ulbricht, M. Toro´s, M. Paternostro, A. Geraci, P. Barker, M. Kim, and G. Milburn, Spin Entanglement Witness for Quantum Gravity, Phys. Rev. Lett. 119, 240401 (2017), arXiv:1707.06050 [quant-ph] .
  * Marshman  _et al._ [2020a] R. J. Marshman, A. Mazumdar, and S. Bose, Locality and entanglement in table-top testing of the quantum nature of linearized gravity, Phys. Rev. A 101, 052110 (2020a).
  * Marletto and Vedral [2017] C. Marletto and V. Vedral, Gravitationally induced entanglement between two massive particles is sufficient evidence of quantum effects in gravity, Phys. Rev. Lett. 119, 240402 (2017).
  * Biswas  _et al._ [2012] T. Biswas, E. Gerwick, T. Koivisto, and A. Mazumdar, Towards singularity and ghost free theories of gravity, Phys. Rev. Lett. 108, 031101 (2012), arXiv:1110.5249 [gr-qc] .
  * van de Kamp  _et al._ [2020] T. W. van de Kamp, R. J. Marshman, S. Bose, and A. Mazumdar, Quantum Gravity Witness via Entanglement of Masses: Casimir Screening, Phys. Rev. A 102, 062807 (2020), arXiv:2006.06931 [quant-ph] .
  * Arndt  _et al._ [1999] M. Arndt, O. Nairz, J. Vos-Andreae, C. Keller, G. van der Zouw, and A. Zeilinger, Wave?particle duality of C60 molecules, Nature 401, 680 (1999).
  * Arndt and Hornberger [2014] M. Arndt and K. Hornberger, Testing the limits of quantum mechanical superpositions, Nature Physics 10, 271 (2014).
  * Toroš  _et al._ [2021] M. Toroš, T. W. van de Kamp, R. J. Marshman, M. S. Kim, A. Mazumdar, and S. Bose, Relative acceleration noise mitigation for nanocrystal matter-wave interferometry: Applications to entangling masses via quantum gravity, Phys. Rev. Research 3, 023178 (2021).
  * Marshman  _et al._ [2020b] R. J. Marshman, A. Mazumdar, G. W. Morley, P. F. Barker, S. Hoekstra, and S. Bose, Mesoscopic Interference for Metric and Curvature (MIMAC) $\&$ Gravitational Wave Detection, New J. Phys. 22, 083012 (2020b), arXiv:1807.10830 [gr-qc] .
  * Leggett [2002] A. J. Leggett, Testing the limits of quantum mechanics: motivation, state of play, prospects, Journal of Physics: Condensed Matter 14, R415 (2002).
  * Bose  _et al._ [1999] S. Bose, K. Jacobs, and P. L. Knight, Scheme to probe the decoherence of a macroscopic object, Physical Review A 59, 3204 (1999).
  * Armour  _et al._ [2002] A. D. Armour, M. P. Blencowe, and K. C. Schwab, Entanglement and decoherence of a micromechanical resonator via coupling to a cooper-pair box, Physical Review Letters 88, 148301 (2002).
  * Marshall  _et al._ [2003] W. Marshall, C. Simon, R. Penrose, and D. Bouwmeester, Towards quantum superpositions of a mirror, Physical Review Letters 91, 130401 (2003).
  * Sekatski  _et al._ [2014] P. Sekatski, M. Aspelmeyer, and N. Sangouard, Macroscopic optomechanics from displaced single-photon entanglement, Physical Review Letters 112, 080502 (2014).
  * Romero-Isart  _et al._ [2010] O. Romero-Isart, M. L. Juan, R. Quidant, and J. I. Cirac, Toward quantum superposition of living organisms, New Journal of Physics 12, 033015 (2010).
  * Romero-Isart  _et al._ [2011] O. Romero-Isart, A. C. Pflanzer, F. Blaser, R. Kaltenbaek, N. Kiesel, M. Aspelmeyer, and J. I. Cirac, Large quantum superpositions and interference of massive nanometer-sized objects, Physical Review Letters 107, 020405 (2011).
  * Khalili  _et al._ [2010] F. Khalili, S. Danilishin, H. Miao, H. Müller-Ebhardt, H. Yang, and Y. Chen, Preparing a mechanical oscillator in non-gaussian quantum states, Physical Review Letters 105, 070403 (2010).
  * Scala  _et al._ [2013] M. Scala, M. S. Kim, G. W. Morley, P. F. Barker, and S. Bose, Matter-wave interferometry of a levitated thermal nano-oscillator induced and probed by a spin, Physical Review Letters 111, 180403 (2013).
  * Wan  _et al._ [2016] C. Wan, M. Scala, G. W. Morley, A. A. Rahman, H. Ulbricht, J. Bateman, P. F. Barker, S. Bose, and M. S. Kim, Free nano-object ramsey interferometry for large quantum superpositions, Phys. Rev. Lett. 117, 143003 (2016).
  * Bateman  _et al._ [2014] J. Bateman, S. Nimmrichter, K. Hornberger, and H. Ulbricht, Near-field interferometry of a free-falling nanoparticle from a point-like source (2014).
  * Yin  _et al._ [2013] Z. Q. Yin, T. Li, X. Zhang, and L. M. Duan, Large quantum superpositions of a levitated nanodiamond through spin-optomechanical coupling, Physical Review A 88, 033614 (2013).
  * Pino  _et al._ [2018] H. Pino, J. Prat-Camps, K. Sinha, B. P. Venkatesh, and O. Romero-Isart, On-chip quantum interference of a superconducting microsphere, Quantum Science and Technology 3, 025001 (2018).
  * Clarke and Vanner [2019] J. Clarke and M. R. Vanner, Growing macroscopic superposition states via cavity quantum optomechanics, Quantum Science and Technology 4, 014003 (2019).
  * Ringbauer  _et al._ [2018] M. Ringbauer, T. J. Weinhold, L. A. Howard, A. G. White, and M. R. Vanner, Generation of mechanical interference fringes by multi-photon counting, New Journal of Physics 20, 053042 (2018).
  * Khosla  _et al._ [2018] K. E. Khosla, M. R. Vanner, N. Ares, and E. A. Laird, Displacemon electromechanics: how to detect quantum interference in a nanomechanical resonator, Physical Review X 8, 021052 (2018).
  * Kaltenbaek  _et al._ [2012] R. Kaltenbaek, G. Hechenblaikner, N. Kiesel, O. Romero-Isart, K. C. Schwab, U. Johann, and M. Aspelmeyer, Macroscopic quantum resonators (maqro), Experimental Astronomy 34, 123 (2012).
  * Romero-Isart [2017] O. Romero-Isart, Coherent inflation for large quantum superpositions of levitated microspheres, New Journal of Physics 19, 123029 (2017).
  * Pedernales  _et al._ [2020] J. S. Pedernales, G. W. Morley, and M. B. Plenio, Motional dynamical decoupling for interferometry with macroscopic particles, Phys. Rev. Lett. 125, 023602 (2020).
  * Hogan  _et al._ [2008] J. M. Hogan, D. Johnson, and M. A. Kasevich, Light-pulse atom interferometry, arXiv preprint arXiv:0806.3261 (2008).
  * Margalit  _et al._ [2021] Y. Margalit, O. Dobkowski, Z. Zhou, O. Amit, Y. Japha, S. Moukouri, D. Rohrlich, A. Mazumdar, S. Bose, C. Henkel, and R. Folman, Realization of a complete stern-gerlach interferometer: Toward a test of quantum gravity, Science Advances 7, eabg2879 (2021), https://www.science.org/doi/pdf/10.1126/sciadv.abg2879 .
  * Al Balushi  _et al._ [2018] A. Al Balushi, W. Cong, and R. B. Mann, Optomechanical quantum cavendish experiment, Physical Review A 98, 043811 (2018).
  * Bassi  _et al._ [2013] A. Bassi, K. Lochan, S. Satin, T. P. Singh, and H. Ulbricht, Models of wave-function collapse, underlying theories, and experimental tests, Reviews of Modern Physics 85, 471 (2013).
  * Kovachy  _et al._ [2015] T. Kovachy, P. Asenbaum, C. Overstreet, C. Donnelly, S. Dickerson, A. Sugarbaker, J. Hogan, and M. Kasevich, Quantum superposition at the half-metre scale, Nature 528, 530 (2015).
  * Murata and Tanaka [2015] J. Murata and S. Tanaka, A review of short-range gravity experiments in the lhc era, Classical and Quantum Gravity 32, 033001 (2015).
  * Wan [2017] C. Wan,  _Quantum superposition on nano-mechanical oscillator_ (Imperial College London (PhD Thesis), 2017).
  * Paraniak and Englert [2021] M. M. Paraniak and B.-G. Englert, Quantum dynamical simulation of a transversal stern?gerlach interferometer, Symmetry 13, 10.3390/sym13091660 (2021).
  * Majorana [1932] E. Majorana, Atomi orientati in campo magnetico variabile, Il Nuovo Cimento (1924-1942) 9, 43 (1932).
  * Inguscio [2007] M. Inguscio, Majorana "spin-flip" and ultra-low temperature atomic physics, PoS EMC2006, 008 (2007).
  * Loubser and van Wyk [1978] J. H. N. Loubser and J. A. van Wyk, Electron spin resonance in the study of diamond, Reports on Progress in Physics 41, 1201 (1978).
  * Kuhn  _et al._ [2017] S. Kuhn, A. Kosloff, B. A. Stickler, F. Patolsky, K. Hornberger, M. Arndt, and J. Millen, Full rotational control of levitated silicon nanorods, Optica 4, 356 (2017).
  * Schäfer  _et al._ [2021] J. Schäfer, H. Rudolph, K. Hornberger, and B. A. Stickler, Cooling nanorotors by elliptic coherent scattering, Phys. Rev. Lett. 126, 163603 (2021).
  * Rudolph  _et al._ [2021] H. Rudolph, J. Schäfer, B. A. Stickler, and K. Hornberger, Theory of nanoparticle cooling by elliptic coherent scattering, Phys. Rev. A 103, 043514 (2021).
  * van der Laan  _et al._ [2021] F. van der Laan, F. Tebbenjohanns, R. Reimann, J. Vijayan, L. Novotny, and M. Frimmer, Sub-kelvin feedback cooling and heating dynamics of an optically levitated librator, Phys. Rev. Lett. 127, 123605 (2021).
  * Machluf  _et al._ [2013] S. Machluf, Y. Japha, and R. Folman, Coherent stern–gerlach momentum splitting on an atom chip, Nature Communications 4, 2424 (2013).
  * Modena  _et al._ [2012] M. Modena, J. Garcia Perez, P. Thonet, C. Petrone, A. Vorozhtsov, O. Dunkel, E. Solodko, and D. Tommasini, Design, assembly and first measurements of a short model for clic final focus hybrid quadrupole qd0, in  _Conf. Proc._ , Vol. 1205201 (2012) p. THPPD010.
  * Bar-Gill  _et al._ [2013] N. Bar-Gill, L. M. Pham, A. Jarmola, D. Budker, and R. L. Walsworth, Solid-state electronic spin coherence time approaching one second, Nature communications 4, 1 (2013).
  * Abobeih  _et al._ [2018] M. H. Abobeih, J. Cramer, M. A. Bakker, N. Kalb, M. Markham, D. J. Twitchen, and T. H. Taminiau, One-second coherence for a single electron spin coupled to a multi-qubit nuclear-spin environment, Nature communications 9, 1 (2018).
  * Muhonen  _et al._ [2014] J. T. Muhonen, J. P. Dehollain, A. Laucht, F. E. Hudson, R. Kalra, T. Sekiguchi, K. M. Itoh, D. N. Jamieson, J. C. McCallum, A. S. Dzurak, _et al._ , Storing quantum information for 30 seconds in a nanoelectronic device, Nature nanotechnology 9, 986 (2014).
  * Farfurnik  _et al._ [2015] D. Farfurnik, A. Jarmola, L. M. Pham, Z. H. Wang, V. V. Dobrovitski, R. L. Walsworth, D. Budker, and N. Bar-Gill, Optimizing a dynamical decoupling protocol for solid-state electronic spin ensembles in diamond, Phys. Rev. B 92, 060301 (2015).
  * Knowles  _et al._ [2014] H. S. Knowles, D. M. Kara, and M. Atatüre, Observing bulk diamond spin coherence in high-purity nanodiamonds, Nature materials 13, 21 (2014).
  * Toroš  _et al._ [2020] M. Toroš, T. van de Kamp, R. Marshman, M. Kim, A. Mazumdar, and S. Bose, Relative Acceleration Noise Mitigation for Entangling Masses via Quantum Gravity, (2020), arXiv:2007.15029 [gr-qc] .
  * Fitzakerley  _et al._ [2016] D. Fitzakerley, M. George, E. Hessels, T. Skinner, C. Storry, M. Weel, G. Gabrielse, C. Hamley, N. Jones, K. Marable, _et al._ , Electron-cooled accumulation of 4$\times$ e9 positrons for production and storage of antihydrogen atoms, Journal of Physics B: Atomic, Molecular and Optical Physics 49, 064001 (2016).
  * Marx  _et al._ [2014] A. Marx, J. Hoess, and K. Uhlig, Dry dilution refrigerator for experiments on quantum effects in the microwave regime, arXiv preprint arXiv:1412.3619 (2014).
  * Schut  _et al._ [2022] M. Schut, J. Tilly, R. J. Marshman, S. Bose, and A. Mazumdar, Improving resilience of quantum-gravity-induced entanglement of masses to decoherence using three superpositions, Physical Review A 105, 032411 (2022).
  * Stickler  _et al._ [2021] B. A. Stickler, K. Hornberger, and M. Kim, Quantum rotations of nanoparticles, Nature Reviews Physics 3, 589 (2021).
  * Wood  _et al._ [2021] A. A. Wood, R. M. Goldblatt, R. E. Scholten, and A. M. Martin, Quantum control of nuclear-spin qubits in a rapidly rotating diamond, Phys. Rev. Research 3, 043174 (2021).
  * Liu  _et al._ [2017] S. Liu, T. Li, and Z. qi Yin, Coupling librational and translational motion of a levitated nanoparticle in an optical cavity, J. Opt. Soc. Am. B 34, C8 (2017).
  * Tebbenjohanns  _et al._ [2021] F. Tebbenjohanns, M. L. Mattana, M. Rossi, M. Frimmer, and L. Novotny, Quantum control of a nanoparticle optically levitated in cryogenic free space, Nature 595, 378 (2021).
  * Delić  _et al._ [2020] U. Delić, M. Reisenbauer, K. Dare, D. Grass, V. Vuletić, N. Kiesel, and M. Aspelmeyer, Cooling of a levitated nanoparticle to the motional quantum ground state, Science 367, 892 (2020).
  * Tebbenjohanns  _et al._ [2020] F. Tebbenjohanns, M. Frimmer, V. Jain, D. Windey, and L. Novotny, Motional sideband asymmetry of a nanoparticle optically levitated in free space, Phys. Rev. Lett. 124, 013603 (2020).
  * Magrini  _et al._ [2021] L. Magrini, P. Rosenzweig, C. Bach, A. Deutschmann-Olek, S. G. Hofer, S. Hong, N. Kiesel, A. Kugi, and M. Aspelmeyer, Real-time optimal quantum control of mechanical motion at room temperature, Nature 595, 373 (2021).
  * Delord  _et al._ [2020] T. Delord, P. Huillery, L. Nicolas, and G. Hétet, Spin-cooling of the motion of a trapped diamond, Nature 580, 56 (2020).
  * Henkel and Folman [2021] C. Henkel and R. Folman, Internal decoherence in nano-object interferometry due to phonons, arXiv preprint arXiv:2112.01263 (2021).
  * Ruster  _et al._ [2016] T. Ruster, C. T. Schmiegelow, H. Kaufmann, C. Warschburger, F. Schmidt-Kaler, and U. G. Poschinger, A long-lived zeeman trapped-ion qubit, Applied Physics B 122, 1 (2016).
  * Gohil  _et al._ [2020] C. Gohil, P. Burrows, N. B. Kraljevic, D. Schulte, and B. Heilig, Measurements of sub-nt dynamic magnetic field shielding with soft iron and mu-metal for use in linear colliders, Journal of Instrumentation 15 (12), P12030.
  * Sochnikov  _et al._ [2020] I. Sochnikov, D. Davino, and B. Kalisky, dc squid design with femtotesla sensitivity for quantum-ready readouts, Physical Review Applied 14, 014020 (2020).
  * Mourou  _et al._ [1981] G. Mourou, C. Stancampiano, A. Antonetti, and A. Orszag, Picosecond microwave pulses generated with a subpicosecond laser-driven semiconductor switch, Applied Physics Letters 39, 295 (1981).

##  Appendix A Solving the particle dynamics

In this supplementary material, we provide our analytical and numerical
treatment of the equations of motion. The magnetic fields used to create the
superposition are

$\displaystyle\vec{B}\left(x,y\right)=$ $\displaystyle\left(B_{0}-\eta
x\right)\hat{x}+\eta y\hat{y}$ (19) $\displaystyle\vec{B}\left(x,y\right)=$
$\displaystyle B_{1}\hat{x}$ (20) $\displaystyle\vec{B}\left(x,y\right)=$
$\displaystyle-\left(B_{0}-\eta x\right)\hat{x}-\eta y\hat{y}$ (21)

with the magnetic field being mapped between Eq. 19 to Eq. 20 when
$t\in\left[\tau_{1},\tau_{2}\right]$, and between Eq. 20 to Eq. 21 when
$t\in\left[\tau_{3},\tau_{4}\right]$. Figure 5 shows the detail of the
experienced magnetic field for each arm of the interferometer during the
period in which the magnetic field is modified. This shows that provided the
value chosen for $B_{1}$ is sufficiently large that the minimum allowable
experienced magnetic field $\varepsilon$ can be avoided while ensuring a
smooth and adiabatic magnetic field transition. When not in this phase
($t\notin\left[\tau_{1},\tau_{4}\right]$) the motion is simply governed by the
harmonic oscillator potential, whose solution is given by:

$x(t)=A\cos\left(\omega t+\phi\right)+C(s_{x})$ (22)

where $A$ and $\phi$ are determined by the initial conditions,
$\omega=\sqrt{-\frac{\chi_{m}}{\mu_{0}}}\eta$ is the frequency of the
diamagnetic trap and $C(s_{x})$ is determined by the spin states and magnetic
field (Eqs. 19 and 21). Note that we are always within the adiabatic limit
throughout these times, such that $\dot{\omega}/\omega^{2}<1$. We also have

$C^{s_{x}=\pm 1}={B_{0}}/{\eta}\pm\alpha\,,$ (23)

where
$\alpha=-\mu_{\pm}\frac{\mu_{0}}{\chi_{m}m\tilde{\eta}(t)}=-\frac{ge\hbar}{2m_{e}}\frac{\mu_{0}}{\chi_{m}m\tilde{\eta}(t)}$.

Figure 5:  Magnetic field experienced by each arm of the interferometer with
time. The vertical grey lines represent times $\tau_{i}$ for $i={1,2,3,4}$
using $B_{0}=10^{-2}$ T, $\eta=40$ Tm-1 and $B_{1}=100\varepsilon$ T.

The trajectory through these times, $t\notin\left[\tau_{1},\tau_{4}\right]$,
can be constructed by simply assembling the solutions piecewise. To do this,
the values of the constants are given by the function of the initial
conditions at each time segment. Specifically, if $x(t_{0})=x_{0}$ and
$\dot{x}(t_{0})=v_{0}$ serve as the initial conditions, then

$\displaystyle A=$
$\displaystyle-\left(C-x_{0}\right)\left(1+\frac{v_{0}^{2}}{\omega^{2}}\left(C-x_{0}\right)^{-2}\right)^{1/2}$
(24) $\displaystyle\phi=$
$\displaystyle\tan^{-1}\left(\frac{v_{0}}{\omega\left(C-x_{0}\right)}\right)-\omega
t_{0}.$ (25)

For example, if we consider the particle to be initialised in a superposition
of spin states $\left|+1\right\rangle$ and $\left|-1\right\rangle$ at the
origin with zero initial velocity, then we can define the initial motion of
the two arms as:

$\displaystyle x_{1}^{+}(t)=$
$\displaystyle-\left(\frac{B_{0}}{\eta}+\alpha\right)\cos\left(\omega
t\right)+\left(\frac{B_{0}}{\eta}+\alpha\right)$ (26) $\displaystyle
x_{1}^{-}(t)=$
$\displaystyle-\left(\frac{B_{0}}{\eta}-\alpha\right)\cos\left(\omega
t\right)+\left(\frac{B_{0}}{\eta}-\alpha\right)\,.$ (27)

When the magnetic field is given by Eq. 20 there will be no spin dependant
acceleration, that is, there is no force acting to either create or destroy
the spatial superposition. Thus we want to minimise $\tau_{4}-\tau_{1}$. To do
this, the initial evolution, $x_{1}(t)$, should be maintained for as long as
possible while still ensuring that

$x\notin\left[\frac{B_{x}(0,0)-\varepsilon}{\partial
B},\frac{B_{x}(0,0)+\varepsilon}{\partial B}\right]\,,$ (28)

is satisfied. Only as the particle approaches the boundary set by Eq.28 is the
magnetic field is modified (this marks the time $\tau_{1}$). Specifically,
this is done as the magnetic field experienced by the forward most trajectory
approaches $B_{1}$:

$x_{1}^{+}(t\approx\tau_{1})=\frac{B_{0}-B_{1}}{\eta}\,,$ (29)

where we write
$x_{1}^{-}(\tau_{1})=A^{-}_{1}\cos\left(\omega\tau_{1}\right)+c^{-}_{1}$.

Similarly, the non-zero magnetic field gradient should be returned as soon as
possible while ensuring Eq. 28 holds. As such, the magnetic field gradient
should begin being restored as the magnetic field experienced by the rearmost
trajectory ($x_{2}^{-}(t)$) approximately the magnetic field as given by Eq.
21, that is:

$\displaystyle B_{1}=$ $\displaystyle-B_{0}+\eta
x^{-}_{2}\left(\tau_{3}\right)$ $\displaystyle
x^{-}_{2}\left(t\approx\tau_{3}\right)=$
$\displaystyle\frac{B_{1}-B_{0}}{\eta}$ (30)

where $x^{+}_{2}$ and $x^{-}_{2}$ are the trajectories when
$t\in\left[\tau_{1},\tau_{4}\right]$, and are found by numerically integrating
the equations of motion.

The final two stages of the trajectories are given by:

$\displaystyle x_{3}^{+}(t)=$ $\displaystyle A^{+}_{3}\cos\left(\omega
t+\phi^{+}_{3}\right)+C^{+}_{3}\,,$ (31) $\displaystyle x_{3}^{-}(t)=$
$\displaystyle A^{-}_{3}\cos\left(\omega t+\phi^{-}_{3}\right)+C^{-}_{3}\,,$
(32) $\displaystyle x_{4}^{+}(t)=$ $\displaystyle A^{+}_{4}\cos\left(\omega
t+\phi^{+}_{4}\right)+C^{+}_{4}\,,$ (33) $\displaystyle x_{4}^{-}(t)=$
$\displaystyle A^{-}_{4}\cos\left(\omega t+\phi^{-}_{4}\right)+C^{-}_{4}\,,$
(34)

where

$\displaystyle A^{+}_{3}=-$
$\displaystyle\left(C^{+}_{3}-x^{+}_{2}(\tau_{2})\right)\left(1+\frac{\left(v^{+}_{2}(\tau_{2})\right)^{2}}{\omega^{2}}\left(C^{+}_{3}-x^{+}_{2}(\tau_{2})\right)^{-2}\right)^{1/2}\,,$
$\displaystyle A^{-}_{3}=-$
$\displaystyle\left(C^{-}_{3}-x^{-}_{2}(\tau_{2})\right)\left(1+\frac{\left(v^{-}_{2}(\tau_{2})\right)^{2}}{\omega^{2}}\left(C^{-}_{3}-x^{-}_{2}(\tau_{2})\right)^{-2}\right)^{1/2}\,,$
$\displaystyle\phi^{+}_{3}=$
$\displaystyle\tan^{-1}\left(\frac{v^{+}_{2}(\tau_{2})}{\omega\left(C^{+}_{3}-x^{+}_{2}(\tau_{2})\right)}\right)-\omega\tau_{2}\,,$
$\displaystyle\phi^{-}_{3}=$
$\displaystyle\tan^{-1}\left(\frac{v^{-}_{2}(\tau_{2})}{\omega\left(C^{+}_{3}-x^{-}_{2}(\tau_{2})\right)}\right)-\omega\tau_{2}\,,~{}~{}$
$\displaystyle C^{+}_{3}=$
$\displaystyle-\left(\frac{B_{0}}{\eta}-\alpha\right)\,,$ $\displaystyle
C^{-}_{3}=$ $\displaystyle-\left(\frac{B_{0}}{\eta}+\alpha\right)\,,$ (35)

and

$\displaystyle A^{+}_{4}=$
$\displaystyle-\left(C^{+}_{4}-x^{+}_{3}(\tau_{3})\right)\left(1+\frac{\left(v^{+}_{3}(\tau_{3})\right)^{2}}{\omega^{2}}\left(C^{+}_{4}-x^{+}_{3}(\tau_{3})\right)^{-2}\right)^{1/2}\,,$
$\displaystyle A^{-}_{4}=$
$\displaystyle-\left(C^{-}_{4}-x^{-}_{3}(\tau_{3})\right)\left(1+\frac{\left(v^{-}_{3}(\tau_{3})\right)^{2}}{\omega^{2}}\left(C^{-}_{3}-x^{-}_{3}(\tau_{3})\right)^{-2}\right)^{1/2}\,,$
$\displaystyle\phi^{+}_{4}=$
$\displaystyle\tan^{-1}\left(\frac{v^{+}_{3}(\tau_{3})}{\omega\left(C^{+}_{4}-x^{+}_{3}(\tau_{3})\right)}\right)-\omega\tau_{3}\,,$
$\displaystyle\phi^{-}_{4}=$
$\displaystyle\tan^{-1}\left(\frac{v^{-}_{3}(\tau_{3})}{\omega\left(C^{+}_{4}-x^{-}_{3}(\tau_{3})\right)}\right)-\omega\tau_{3}\,,$
$\displaystyle C^{+}_{4}=$
$\displaystyle-\left(\frac{B_{0}}{\eta}+\alpha\right)\,,$ $\displaystyle
C^{-}_{4}=$ $\displaystyle-\left(\frac{B_{0}}{\eta}-\alpha\right).$ (36)

The values for the times $\tau_{5}$ and $\tau_{6}$ are fixed by the following
conditions:

  1. 1.

$\Delta x(\tau_{6})=x_{4}^{+}(\tau_{6})-x_{4}^{-}(\tau_{6})=0$ and $\Delta
v(\tau_{6})=v_{4}^{+}(\tau_{6})-v_{4}^{-}(\tau_{6})=0$, such that the two arms
of the interferometers are brought together to overlap in both the position
and the momentum space, respectively.

  2. 2.

$x_{4}^{+}(\tau_{6})>\frac{B_{0}+\varepsilon}{\eta}$, such that the Majorana
spin flip region of the magnetic field is again avoided.

Note that the appropriate time scales which met the above conditions
(conditions 1. and 2.) were found solely numerically. This has lead to a set
of times and the corresponding magnetic field gradients applied during the
time period, which we have tabulated in the table 1.

Table 1: $\tau_{i}$ values used to calculate wave-packet trajectories for example magnetic field gradients ($\eta$)  |  $\eta=4$ Tm-1 |  $\eta=40$ Tm-1 |  $\eta=400$ Tm-1
---|---|---|---
$\tau_{1}$ | 5.39 s | 0.534 s | 0.0493 s
$\tau_{2}$ | 5.39 s | 0.539 s | 0.0539 s
$\tau_{3}$ | 5.80 s | 0.580 s | 0.0580 s
$\tau_{4}$ | 5.80 s | 0.584 s | 0.0626 s
$\tau_{5}$ | 9.01 s | 0.902 s | 0.0913 s
$\tau_{6}$ | 14.8 s | 1.48 s | 0.148 s

There are a couple of points to note here; the time $\tau_{6}$, set such that
$\Delta v(\tau_{6})=0$ will automatically minimise $\Delta x(\tau_{6})$. Also,
the value of $\Delta x$ evaluated at time $\tau_{6}$ is continuous in
$\tau_{5}$, and there exists times $t_{1}$ and $t_{2}\in\mathbb{R}$ such that,
when $\tau_{5}=t_{1}~{}\Delta x(\tau_{6})>0$ and when $\tau_{5}=t_{2}~{}\Delta
x(\tau_{6})<0$. In a nut-shell the procedure to find the correct values of
$\tau_{5}$ and $\tau_{6}$ are:

  1. 1.

To make an initial guess for the value of $\tau_{5}$ which was used to
calculate a complete trajectory, typically this is $\tau_{5}=2\tau_{4}$.

  2. 2.

From this we have evaluated the value of $\tau_{6}$, which we have determined
using the relation $\Delta v(\tau_{6})=0$, and

  3. 3.

the corresponding value of $\Delta x(\tau_{6})$ was then evaluated. This
allows the the assumed value for $\tau_{5}$ to be optimised accordingly.
Specifically noting that increasing $\tau_{5}$ will lead to decrease the value
of $\Delta x(\tau_{6})$.


Supplementary Materials for

A Variational Approach to Trap Macromolecules

Luofu Liu1,2, Chao Duan1, and Rui Wang1,2,∗

1\. Department of Chemical and Biomolecular Engineering, University of
California Berkeley, CA 94720, USA

2\. Materials Sciences Division, Lawrence Berkeley National Lab, Berkeley,
California 94720, USA

* correspondence to<EMAIL_ADDRESS>
##  I I. Derivation of the Variational Approach for Trapping the Center-of-
Mass of Polymers

In this section, we provide a detailed derivation of the variational approach
for trapping the center of mass (c.m.) of polymers. The formulation for a
single neutral polymer is in given Sec. 1.1. The generalization of the
approach to a charged polymer is provided in Sec. 1.2, and the generalization
to the simultaneous trap of two polymers is provided in Sec. 1.3.

1.1. Variational Approach for a Single Neutral Polymer

We consider a system of one polymer immersed in $n_{\rm s}$ solvent molecules
which are connected with a pure solvent reservoir to maintain the chemical
potential of solvent $\mu_{\rm s}$. The polymer is modeled by a continuous
Gaussian chain of $N$ Kuhn segments with Kuhn length $b$. To trap this
polymer, we fix its c.m. at a specific position $\bm{\xi}$, giving rise to the
semicanonical partition function

$\displaystyle Z=\sum_{n_{\rm s}}\frac{e^{\beta\mu_{\rm s}n_{\rm
s}}}{n_{s}!v_{\rm p}^{N}v_{\rm s}^{n_{\rm
s}}}\int{\rm\hat{D}}\\{\bm{R}\\}\prod_{j=1}^{n_{\rm s}}\int{\rm
d}\bm{r}_{j}\exp(-\beta H)\prod_{\bm{r}}\delta[1-\hat{\phi}_{\rm
p}(\bm{r})-\hat{\phi}_{\rm s}(\bm{r})]\delta\left[\frac{1}{N}\int_{0}^{N}{\rm
d}s\bm{R}(s)-\bm{\xi}\right]$ (S1)

where $\int{\rm\hat{D}}\\{\bm{R}\\}=\int{\rm D}\\{\bm{R}\\}{\rm
exp}[{-(3/2b^{2})\int^{N}_{0}ds(\partial{\bm{R}}}/{\partial s})^{2}]$ denotes
the integration over all configurations weighted by Gaussian-chain statistics.
The Hamiltonian

$\displaystyle\beta H=\frac{\chi}{v}\int{\rm d}\bm{r}\hat{\phi}_{\rm
p}(\bm{r})\hat{\phi}_{\rm s}(\bm{r})$ (S2)

includes the polymer-solvent interaction in terms of the Flory $\chi$
parameter, with $\hat{\phi}_{\rm p}$ and $\hat{\phi}_{\rm s}$ the
instantaneous volume fractions of the polymer and solvent, respectively. For
simplicity, the volumes of polymer segments and solvent molecules are taken to
be the same, $v_{\rm p}=v_{\rm s}=v$. The first $\delta$-functional accounts
for the incompressibility whereas the second one enforces the constraint on
the c.m. of the polymer.

The transformation from the particle-based to the field-based representation
is achieved by inserting the following identities into the partition function:

$\displaystyle 1=$ $\displaystyle\int{\rm
D}\phi_{\kappa}\prod_{\bm{r}}\delta[\phi_{\kappa}({\bm{r}})-\hat{\phi}_{\kappa}({\bm{r}})]=\int{\rm
D}\phi_{\kappa}{\rm D}W_{\kappa}\exp\left\\{i\int
d{\bm{r}}W_{\kappa}({\bm{r}})[\phi_{\kappa}({\bm{r}})-\hat{\phi}_{\kappa}({\bm{r}})]\right\\}\leavevmode\nobreak\
\leavevmode\nobreak\ (\kappa={\rm p,s})$ (S3)

where the right-hand side of the equation arises from the Fourier
representation of the $\delta$-function with $W_{\kappa}({\bm{r}})$ being the
Fourier conjugate field to $\phi_{\kappa}({\bm{r}})$. Similarly, the Fourier
representations of the incompressibility condition and constraint of the
polymer c.m. are

$\displaystyle\delta[1-\hat{\phi}_{\rm p}(\bm{r})-\hat{\phi}_{\rm
s}(\bm{r})]=\int{\rm D}\Gamma\exp\left\\{i\int{\rm
d}{\bm{r}}\Gamma({\bm{r}})\left[1-\hat{\phi}_{\rm p}(\bm{r})-\hat{\phi}_{\rm
s}(\bm{r})\right]\right\\}$ (S4)
$\displaystyle\delta\left[\frac{1}{N}\int_{0}^{N}{\rm
d}s{\bm{R}}(s)-{\bm{\xi}}\right]=\int{\rm
d}{\bm{F}}\exp\left\\{i{\bm{F}}\cdot\left[\frac{1}{N}\int_{0}^{N}{\rm
d}s{\bm{R}}(s)-{\bm{\xi}}\right]\right\\}$ (S5)

By performing the identity transformations, we obtain the field-based
partition function as

$\displaystyle Z=\int{\rm D}\phi_{\rm p}{\rm D}\phi_{\rm s}{\rm D}W_{\rm
p}{\rm D}W_{\rm s}{\rm D}\Gamma{\rm d}{\bm{F}}\sum_{n_{\rm
s}}\frac{e^{\beta\mu_{\rm s}n_{\rm s}}}{n_{s}!v^{N+n_{\rm
s}}}\int{\rm\hat{D}}\\{\bm{R}\\}\prod_{j=1}^{n_{\rm s}}\int{\rm d}\bm{r}_{j}$
$\displaystyle\cdot\exp\left\\{\frac{1}{v}\int{\rm
d}{\bm{r}}\left[-\chi\phi_{\rm p}({\bm{r}})\phi_{\rm s}({\bm{r}})+iW_{\rm
p}({\bm{r}})\left(\phi_{\rm p}({\bm{r}})-\hat{\phi}_{\rm
p}({\bm{r}})\right)+iW_{\rm s}({\bm{r}})\left(\phi_{\rm
s}({\bm{r}})-\hat{\phi}_{\rm
s}({\bm{r}})\right)+i\Gamma({\bm{r}})\left(\phi_{\rm p}({\bm{r}})+\phi_{\rm
s}({\bm{r}})-1\right)\right]\right\\}$
$\displaystyle\cdot\exp\left\\{i{\bm{F}}\cdot\left[\frac{1}{N}\int_{0}^{N}{\rm
d}s{\bm{R}}(s)-{\bm{\xi}}\right]\right\\}$ (S6)

Here, Eq. I can be decomposed as $Z=Z_{0}Z_{F}$, where $Z_{0}$ denotes the
contribution irrelevant to the trap and also exists in other constraint-free
systems. $Z_{F}$ comes from the constraint of the c.m. They are given by

$\displaystyle Z_{0}=\int{\rm D}\phi_{\rm p}{\rm D}\phi_{\rm s}{\rm D}W_{\rm
p}{\rm D}W_{\rm s}{\rm D}\Gamma\exp\left\\{e^{\beta\mu_{\rm s}}Q_{\rm
s}+\frac{1}{v}\int{\rm d}{\bm{r}}\left[-\chi\phi_{\rm p}\phi_{\rm s}+iW_{\rm
p}\phi_{\rm p}+iW_{\rm s}\phi_{\rm s}+i\Gamma\left(\phi_{\rm p}+\phi_{\rm
s}-1\right)\right]\right\\}$ (S7a) $\displaystyle Z_{F}=\int{\rm
d}{\bm{F}}\exp(-L_{F})$ (S7b)

where $Q_{\rm s}=\int{\rm d}{\bm{r}}e^{-iW_{\rm s}({\bm{r}})}$ is the
partition function of solvents. The “action” $L_{F}=-\ln Q_{\rm P}$, with
$Q_{\rm P}$ the single-chain partition function in the auxiliary fields
$W_{\rm p}$ and $\bm{F}$:

$Q_{\rm
p}=\frac{1}{v^{N}}\int{\rm\hat{D}}\\{{\bm{R}}\\}\exp\bigg{\\{}-\int_{0}^{N}{\rm
d}s\bigg{[}iW_{\rm
p}({\bm{R}}(s))-i\frac{\bm{F}}{N}\cdot({\bm{R}}(s)-{\bm{\xi}})\bigg{]}\bigg{\\}}$
(S8)

To focus on the effect of the c.m. fluctuation, we perform a variational
treatment on $Z_{F}$ while taking the saddle-point approximation to evaluate
the functional integrals of fields included in $Z_{0}$ [1]. For the latter
case, $W_{\rm p}$, $W_{\rm s}$, and $\Gamma$ will be replaced by their saddle-
point values $-iw_{\rm p}$, $-iw_{\rm s}$, and $-i\gamma$. The free energy
corresponding to $Z_{0}$ is given by

$\displaystyle W_{0}$ $\displaystyle=-\ln Z_{0}$
$\displaystyle\approx-e^{\beta\mu_{\rm s}}Q_{\rm s}+\frac{1}{v}\int{\rm
d}{\bm{r}}\left\\{\chi\phi_{\rm p}(\bm{r})\phi_{\rm s}(\bm{r})-w_{\rm
p}(\bm{r})\phi_{\rm p}(\bm{r})-w_{\rm s}(\bm{r})\phi_{\rm
s}(\bm{r})+\gamma(\bm{r})\left[\phi_{\rm p}(\bm{r})+\phi_{\rm
s}(\bm{r})-1\right]\right\\}$ (S9)

For $Z_{F}$, the saddle-point contribution of ${\bm{F}}$, as used in previous
work [2, 3, 4], fails to capture the fluctuation of the c.m. which cannot be
ignored as the polymer takes a coil state in $\theta$ and good solvents. To
self-consistently account for this fluctuation effect, the Gibbs-Feynman-
Bogoliubov variational approach is used to estimate the integral of $\bm{F}$
for evaluating $Z_{F}$ [5]:

$Z_{F}=Z_{\rm ref}\langle\exp[-(L_{F}-L_{\rm ref})]\rangle_{\rm ref}\approx
Z_{\rm ref}\exp[-\langle L_{F}-L_{\rm ref}\rangle_{\rm ref}]$ (S10)

where the average $\langle...\rangle_{\rm ref}$ is taken in the reference
ensemble with action $L_{\rm ref}$. We take the reference action to be the
$2n$-power modified Gaussian form centered around the average force
$-i{\bm{f}}$:

$e^{-L_{\rm
ref}[{\bm{F}},{\bm{f}},{\bm{A}}]}=\prod_{\kappa=x,y,z}(F_{\kappa}+if_{\kappa})^{2n}\exp\left[-\frac{(F_{\kappa}+if_{\kappa})^{2}}{2A_{\kappa}}\right]$
(S11)

$Z_{\rm ref}=\int{\rm d}\bm{F}\exp(-L_{\rm
ref})=\prod_{\kappa}(2\pi)^{1/2}(2n-1)!!A^{n+1/2}_{\kappa}$ is the
normalization factor. The average force $\bm{f}$ and coefficients
${A}_{\kappa}$ are taken to be the variational parameters. This construction
of $L_{\rm ref}$ assumes the independence of the c.m. fluctuations in the
three directions. Note that a general $2n$-power modified Gaussian is adopted
here instead of a standard ($0$th-power) Gaussian to reinforce the confinement
of the c.m. Although the variational approach using the standard Gaussian [6,
7, 8] can also lead to the correct scaling (see Fig. S2), this reinforcement
is necessary for an accurate prediction of the coil size in $\theta$ and good
solvents (see Fig. S1). Evaluation of $e^{-\langle L_{F}\rangle_{\rm ref}}$
gives an approximation of the variational single-chain partition function
$Q_{\rm var}$ as follows:

$\displaystyle Q_{\rm var}$ $\displaystyle=\frac{1}{v^{N}}\int{\rm
D}\\{{\bm{R}}\\}\exp\left\\{-\frac{3}{2b^{2}}\int_{0}^{N}{\rm
d}s\left[\frac{\partial{\bm{R}}(s)}{\partial s}\right]^{2}-\int_{0}^{N}{\rm
d}sw_{\rm
p}[{\bm{R}}(s)]+{\bm{f}}\cdot\overline{\Delta{\bm{R}}}\right\\}\langle\exp\left(i{\bm{g}}\cdot\overline{\Delta{\bm{R}}}\right)\rangle_{\rm
ref}$ (S12)

where ${\bm{g}}={\bm{F}}+i{\bm{f}}$, and
$\overline{\Delta{\bm{R}}}=N^{-1}\int_{0}^{N}{\rm d}s[{\bm{R}}(s)-{\bm{\xi}}]$
is the difference of the instantaneous c.m. from the targeted position
$\bm{\xi}$. The average in the last term can be evaluated via the Hubbard-
Stratonovich transform [1] as

$\displaystyle\langle\exp\left(i{\bm{g}}\cdot\overline{\Delta{\bm{R}}}\right)\rangle_{\rm
ref}$ $\displaystyle=\frac{1}{Z_{\rm ref}}\prod_{\kappa=x,y,z}\int{\rm
d}g_{\kappa}g^{2n}_{\kappa}\exp\left(ig_{\kappa}\overline{\Delta
R}_{\kappa}-\frac{g^{2}_{\kappa}}{2A_{\kappa}}\right)$
$\displaystyle=\prod_{\kappa=x,y,z}\exp\left(-\frac{A_{\kappa}\overline{\Delta
R}^{2}_{\kappa}}{2}\right)\sum^{n}_{m=0}C_{n,m}(A_{\kappa}\overline{\Delta
R}^{2}_{\kappa})^{m}$ (S13)

with the coefficient $C_{n,m}=(-1)^{m}[(2n)!(2n-2m-1)!!]/[(2m)!(2n-1)!!]$. The
variational single-chain partition function $Q_{\rm var}$ in Eq. S12 can then
be written as

$\displaystyle Q_{\rm var}$ $\displaystyle=\frac{1}{v^{N}}\int{\rm
D}\\{{\bm{R}}\\}e^{-H[\bm{R}]}\left[\prod_{\kappa=x,y,z}\left(\sum^{n}_{m=0}C_{n,m}A^{m}_{\kappa}\overline{\Delta
R}^{2m}_{\kappa}\right)\right]$ (S14)

with the Hamiltonian $H[\bm{R}]$ defined by

$\displaystyle H[{\bm{R}}]$ $\displaystyle=\frac{3}{2b^{2}}\int_{0}^{N}{\rm
d}s\left[\frac{\partial{\bm{R}}(s)}{\partial s}\right]^{2}+\int_{0}^{N}{\rm
d}sw_{\rm
p}[{\bm{R}}(s)]-{\bm{f}}\cdot\overline{\Delta{\bm{R}}}+\sum_{\kappa=x,y,z}\frac{A_{\kappa}}{2}\overline{\Delta
R}^{2}_{\kappa}$ (S15)

Note that the complete variational single-chain Hamiltonian in $Q_{\rm var}$
should be
$H[{\bm{R}}]-\sum_{\kappa}\ln(\sum^{n}_{m=0}C_{n,m}A^{m}_{\kappa}\overline{\Delta
R}^{2m}_{\kappa})$. Here, we separate $H[\bm{R}]$ from the complete
Hamiltonian to facilitate the evaluation of $Q_{\rm var}$. In this sense,
$Q_{\rm var}$ can be re-expressed as an ensemble average based on $H[\bm{R}]$:

$\displaystyle Q_{\rm var}$ $\displaystyle=\frac{v^{-N}\int{\rm
D}\\{{\bm{R}}\\}e^{-H[\bm{R}]}\left[\prod_{\kappa=x,y,z}\left(\sum^{n}_{m=0}C_{n,m}A^{m}_{\kappa}\overline{\Delta
R}^{2m}_{\kappa}\right)\right]}{v^{-N}\int{\rm
D}\\{{\bm{R}}\\}e^{-H[\bm{R}]}}\frac{1}{v^{N}}\int{\rm
D}\\{{\bm{R}}\\}e^{-H[\bm{R}]}$
$\displaystyle=\left[\prod_{\kappa=x,y,z}\left(\sum^{n}_{m=0}C_{n,m}A^{m}_{\kappa}\langle\overline{\Delta
R}^{2m}_{\kappa}\rangle_{H}\right)\right]Q_{H}$ (S16)

where $Q_{H}$ is the single-chain partition function with the Hamiltonian
$H[\bm{R}]$

$\displaystyle Q_{H}=\frac{1}{v^{N}}\int{\rm D}\\{{\bm{R}}\\}e^{-H[\bm{R}]}$
(S17)

A common way to evaluate $Q_{H}$ requires replacing the complicated functional
integral with the product of chain propagators. The chain propagators satisfy
the modified diffusion equation [1]. However, this becomes difficult here due
to the nonlinear term appeared in $A_{\kappa}\overline{\Delta R}^{2}_{\kappa}$
in $H[\bm{R}]$ (Eq. S15) which contains 2-body contributions. To circumvent
this difficulty, we decompose the Hamiltonian $H[\bm{R}]$ into a 1-body term
$H^{(1)}[{\bm{R}}(s)]$ and a 2-body term $H^{(2)}[{\bm{R}}(s),{\bm{R}}(t)]$:

$\displaystyle H^{(1)}[{\bm{R}}(s)]$ $\displaystyle=\int_{0}^{N}{\rm
d}s\left\\{\frac{3}{2b^{2}}\left[\frac{\partial{\bm{R}}(s)}{\partial
s}\right]^{2}+w_{\rm
p}[{\bm{R}}(s)]-\frac{\bm{f}}{N}\cdot[{\bm{R}}(s)-{\bm{\xi}}]+\sum_{\kappa=x,y,z}\frac{A_{\kappa}}{2}[R_{\kappa}(s)-\xi_{\kappa}]^{2}\right\\}$
(S18) $\displaystyle H^{(2)}[{\bm{R}}(s),{\bm{R}}(t)]$
$\displaystyle=\int_{0}^{N}{\rm d}s\int_{0}^{N}{\rm
d}^{\prime}t\left\\{\sum_{\kappa=x,y,z}\frac{A_{\kappa}}{2}[R_{\kappa}(s)-\xi_{\kappa}][R_{\kappa}(t)-\xi_{\kappa}]\right\\}$
(S19)

where the prime in $\int^{N}_{0}{\rm d}^{\prime}t$ indicates the case of $t=s$
is excluded in the integral. Then we can approximate $Q_{H}$ by evaluating the
2-body part $H^{(2)}[{\bm{R}}(s),{\bm{R}}(t)]$ via an averaging under the
reference ensemble with 1-body part $H^{(1)}[{\bm{R}}(s)]$:

$\displaystyle Q_{H}$ $\displaystyle=\int{\rm
D}{\bm{R}}e^{-H^{(1)}[{\bm{R}}(s)]-H^{(2)}[{\bm{R}}(s),{\bm{R}}(t)]}=Q^{(1)}_{H}{\langle
e^{-H^{(2)}[{\bm{R}}(s),{\bm{R}}(t)]}\rangle}_{H^{(1)}}$ (S20)

with $Q^{(1)}_{H}=\int D{\bm{R}}e^{-H^{(1)}[{\bm{R}}(s)]}$. To evaluate the
two-body Hamiltonian in the one-body reference in the right-hand side of Eq.
S20, we first fix the position of the $s$-th monomer ${\bm{R}}(s)$ and
calculate the average with respect to the position of the $t$-th monomer
${\bm{R}}(t)$ to obtain ${\langle
e^{-H^{(2)}[{\bm{R}}(s),\overline{{\bm{R}}(t)}]}\rangle}_{H^{(1)}}$. We then
take the average with respect to ${\bm{R}}(s)$ as:

$\displaystyle Q_{H}$ $\displaystyle=\int{\rm
D}{\bm{R}}e^{-H^{(1)}[{\bm{R}}(s)]}{\langle
e^{-H^{(2)}[{\bm{R}}(s),\overline{{\bm{R}}(t)}]}\rangle}_{H^{(1)}}\approx\int{\rm
D}{\bm{R}}e^{-H^{(1)}[{\bm{R}}(s)]-{\langle
H^{(2)}[{\bm{R}}(s),\overline{{\bm{R}}(t)}]\rangle}_{H^{(1)}}}$ (S21)

where the second equality assumes that higher order $2n$-body correlations
($n\geq 2$) are neglected. Then $Q_{H}$ can be evaluated via the forward and
backward chain propagators $q({\bm{r}};s)$, $q^{\dagger}({\bm{r}};s)$

$\displaystyle Q_{H}$ $\displaystyle=\frac{1}{v}\int
d{\bm{r}}q^{\dagger}({\bm{r}};s)q({\bm{r}};s)$ (S22)
$\displaystyle\frac{\partial q({\bm{r}};s)}{\partial
s}=\frac{b^{2}}{6}\nabla^{2}q({\bm{r}};s)-U_{\rm p}({\bm{r}};s)q({\bm{r}};s)$
(S23a) $\displaystyle-\frac{\partial q^{\dagger}({\bm{r}};s)}{\partial
s}=\frac{b^{2}}{6}\nabla^{2}q^{\dagger}({\bm{r}};s)-U_{\rm
p}({\bm{r}};s)q^{\dagger}({\bm{r}};s)$ (S23b)

where the total interaction field experienced by the $s$-th monomer $U_{\rm
p}({\bm{r}};s)=U^{(1)}_{\rm p}({\bm{r}})+U^{(2)}_{\rm p}({\bm{r}};s)$ with the
1-body interaction

$\displaystyle U^{(1)}_{\rm p}({\bm{r}})=w_{\rm
p}({\bm{r}})-\frac{\bm{f}}{N}\cdot({\bm{r}}-{\bm{\xi}})+\sum_{\kappa=x,y,z}\frac{A_{\kappa}}{2N}(\kappa-\xi_{\kappa})^{2}$
(S24)

and the 2-body interaction

$\displaystyle U^{(2)}_{\rm p}({\bm{r}};s)$
$\displaystyle=\sum_{\kappa=x,y,z}\frac{A_{\kappa}}{2N^{2}}(\kappa-\xi_{\kappa})\Biggl{\\{}\int^{s}_{0}{\rm
d}t\frac{1}{q_{(1)}({\bm{r}};s)}\int{\rm
d}{\bm{r}^{\prime}}(\kappa^{\prime}-\kappa)q_{(1)}({\bm{r}}^{\prime};t)g_{(1)}({\bm{r}}^{\prime},{\bm{r}};t,s)$
$\displaystyle+\int^{N}_{s}{\rm
d}t\frac{1}{q^{\dagger}_{(1)}({\bm{r}};s)}\int{\rm
d}{\bm{r}^{\prime}}(\kappa^{\prime}-\kappa)g_{(1)}({\bm{r}},{\bm{r}}^{\prime};s,t)q^{\dagger}_{(1)}({\bm{r}}^{\prime};t)\Biggl{\\}}$
(S25)

where the chain propagator $q_{(1)}({\bm{r}};s)$
($q^{\dagger}_{(1)}({\bm{r}};s)$) and intra-chain correlation function
$g_{(1)}({\bm{r}},{\bm{r}}^{\prime};s,t)$
($g_{(1)}({\bm{r}}^{\prime},{\bm{r}};t,s)$) are subject to the same recurrence
relation (Eq. S23a (Eq. S23b)) as $q({\bm{r}};s)$ ($q^{\dagger}({\bm{r}};s)$)
but with interaction field $U^{(1)}_{\rm p}$.

Then the variational free energy corresponding to $Z_{F}$ is then given by

$\displaystyle W_{F}$ $\displaystyle=-\ln Z_{F}$ $\displaystyle\approx-\ln
Q_{H}-\frac{1}{2}\sum_{\kappa=x,y,z}\ln
A_{\kappa}-\sum_{\kappa=x,y,z}\ln\left(\sum^{n}_{m=0}C_{n,m}A^{m}_{\kappa}\langle\overline{\Delta
R}^{2m}_{\kappa}\rangle_{H}\right)$ (S26)

where the second term in the second equality comes from $-\ln Z_{\rm
ref}=-(n+1/2)\sum_{\kappa}\ln A_{\kappa}$ and $-\langle L_{\rm
ref}\rangle_{\rm ref}=n\sum_{\kappa}\ln A_{\kappa}$. Note that we ignore the
constants independent of the variational parameters. The last term is from the
$2n$-power in the modified Gaussian reference. Using $\langle
x^{2m}\rangle=\int{\rm d}xx^{2m}e^{-ax^{2}/2+bx}/\int{\rm
d}xe^{-ax^{2}/2+bx}=\sum^{m}_{l}(-1)^{l}(2m-1)!!C_{m,l}b^{2l}a^{-(m+l)}$, the
moments $\langle\overline{\Delta R}^{2m}_{\kappa}\rangle_{H}$ can be estimated
by neglecting the elastic energy and interaction field $w_{\rm p}$ in Eq. S15.
Thereafter, the last term in Eq. I becomes a logarithm of linear superposition
of power functions for the variational parameter $A_{\kappa}$, in the form of
$-\sum_{\kappa}\ln(\sum^{n}_{p=0}d_{\kappa,p}A^{-p}_{\kappa})$. $d_{\kappa,p}$
are prefactors independent of $A_{\kappa}^{-p}$.

Minimizing $W_{F}$ with respect to $A_{\kappa}$, we obtain

$\displaystyle-\frac{1}{Q_{H}}\frac{\partial Q_{H}}{\partial
A_{\kappa}}=\frac{1}{A_{\kappa}}\left[\frac{1}{2}+\frac{\sum^{n}_{p=0}(-p)d_{\kappa,p}A^{-p}_{\kappa}}{\sum^{n}_{p=0}d_{\kappa,p}A^{-p}_{\kappa}}\right]\leavevmode\nobreak\
\leavevmode\nobreak\ (\kappa=x,y,z)$ (S27)

with the help of Eq. S15, the left-hand side yields

$\displaystyle-\frac{1}{Q_{H}}\frac{\partial Q_{H}}{\partial
A_{\kappa}}=\frac{1}{2}\left[R^{2}_{g,\kappa}+\delta^{(2)}_{\kappa}\right]$
(S28a) $\displaystyle R^{2}_{g,\kappa}=\frac{\int{\rm
d}{\bm{r}}(\kappa-\xi_{\kappa})^{2}\phi_{\rm p}(\bm{r})}{\int{\rm
d}{\bm{r}}\phi_{\rm p}(\bm{r})}$ (S28b)
$\displaystyle\delta^{(2)}_{\kappa}=\frac{1}{\int{\rm d}{\bm{r}}\phi_{\rm
p}(\bm{r})}\int{\rm d}{\bm{r}}(\kappa-\xi_{\kappa})\phi_{\rm
p}(\bm{r})\Biggl{\\{}\int^{s}_{0}{\rm d}t\frac{1}{q_{(1)}({\bm{r}};s)}\int{\rm
d}{\bm{r}^{\prime}}(\kappa^{\prime}-\kappa)q_{(1)}({\bm{r}}^{\prime};t)g_{(1)}({\bm{r}}^{\prime},{\bm{r}};t,s)$
$\displaystyle+\int^{N}_{s}{\rm
d}t\frac{1}{q^{\dagger}_{(1)}({\bm{r}};s)}\int{\rm
d}{\bm{r}^{\prime}}(\kappa^{\prime}-\kappa)g_{(1)}({\bm{r}},{\bm{r}}^{\prime};s,t)q^{\dagger}_{(1)}({\bm{r}}^{\prime};t)\Biggl{\\}}$
(S28c)

where $\delta^{(2)}_{\kappa}$ comes from the 2-body interaction $U^{(2)}_{\rm
p}({\bm{r}};s)$ in Eq. I. For the right-hand side of Eq. S27, we neglect the
$A_{\kappa}$-dependence of the ratio
$r=\sum^{n}_{p=0}(-p)d_{\kappa,p}A^{-p}_{\kappa}/\sum^{n}_{p=0}d_{\kappa,p}A^{-p}_{\kappa}$,
and treat $r$ as a general fitting parameter for simplicity. Then Eq. S27
becomes

$\displaystyle
A_{\kappa}=\frac{\lambda}{R^{2}_{g,\kappa}+\delta^{(2)}_{\kappa}}\leavevmode\nobreak\
\leavevmode\nobreak\ (\kappa=x,y,z)$ (S29)

with $\lambda=1+2r$. Similarly, $W_{F}\approx-\ln
Q_{H}-(\lambda/2)\sum_{\kappa}\ln A_{\kappa}$. Especially, if we take the
standard (0-th power) Gaussian reference, we have $\lambda=1$. While the ratio
$r$ arising from the $2n$-power modified Gaussian makes $\lambda$ adjustable,
so that we can calibrate its value based on the criterion of the known size of
the infinitely long ideal chain: $R_{\rm g}^{2}=Nb^{2}/6$ at $\chi=0.5$ as
$N\to\infty$.

Taking the saddle-point approximation for the density of polymer and solvents
$\phi_{\rm p}(\bm{r})$, $\phi_{\rm s}(\bm{r})$, conjugate fields $w_{\rm
p}(\bm{r})$, $w_{\rm s}(\bm{r})$, incompressibility field $\gamma(\bm{r})$,
and optimizing $W_{F}$ with respect to the average force ${\bm{f}}$, we obtain

$\displaystyle w_{\rm p}(\bm{r})=\chi\phi_{\rm s}(\bm{r})+\gamma(\bm{r})$
(S30a) $\displaystyle w_{\rm s}(\bm{r})=\chi\phi_{\rm
p}(\bm{r})+\gamma(\bm{r})$ (S30b) $\displaystyle\phi_{\rm
p}(\bm{r})=\frac{1}{Q_{H}}\int{\rm d}sq^{\dagger}(\bm{r};s)q(\bm{r};s)$ (S30c)
$\displaystyle\phi_{\rm s}(\bm{r})=\exp[\beta\mu_{\rm s}-w_{\rm s}(\bm{r})]$
(S30d) $\displaystyle\phi_{\rm p}(\bm{r})+\phi_{\rm s}(\bm{r})-1=0$ (S30e)
$\displaystyle\int{\rm d}\bm{r}(\bm{r}-\bm{\xi})\phi_{\rm p}(\bm{r})=\bm{0}$
(S30f)

The resulting free energy of the system is

$\displaystyle W$ $\displaystyle=W_{0}+W_{F}=-\ln Q_{H}-e^{\beta\mu_{\rm
s}}Q_{\rm s}-\frac{\lambda}{2}\sum_{\kappa=x,y,z}\ln A_{\kappa}$
$\displaystyle+\frac{1}{v}\int{\rm d}\bm{r}\left\\{\chi\phi_{\rm
p}(\bm{r})\phi_{\rm s}(\bm{r})-w_{\rm p}(\bm{r})\phi_{\rm p}(\bm{r})-w_{\rm
s}(\bm{r})\phi_{\rm s}(\bm{r})+\gamma(\bm{r})\left[\phi_{\rm
p}(\bm{r})+\phi_{\rm s}(\bm{r})-1\right]\right\\}$ (S31)

In most cases, we can safely assume that the two-body contribution is
negligible compared to the one-body term such that $U^{(2)}_{\rm
p}({\bm{r}};s)\ll U^{(1)}_{\rm p}({\bm{r}})$ and $\delta_{\kappa}^{(2)}\ll
R_{\rm g,\kappa}^{2}$. Under this assumption, Eqs. S29, S23a, S24, and I
become Eqs. 6e-9 in the main text.

1.2. Generalization to a Charged Polymer

We consider a semicanonical ensemble consisting of one polyelectrolyte (PE)
chain connected to a reservoir of a salt solution with bulk ion concentration
$c^{b}_{\pm}$ that maintains the chemical potentials of solvent $\mu_{s}$ and
ions $\mu_{\pm}$. The PE is assumed to be a homogeneous Gaussian chain with
smeared backbone charge density $\alpha$ [9]. Mobile ions are taken as point
charges with valency $z_{\pm}$. The semicanonical partition function can be
written as

$\displaystyle Z=\frac{1}{v_{\rm
p}^{N}}\int{\rm\hat{D}}\\{\bm{R}\\}\prod_{\zeta={\rm s},\pm}\sum_{n_{\rm
s}}\frac{e^{\beta\mu_{\zeta}n_{\zeta}}}{n_{\zeta}!v_{\zeta}^{n_{\zeta}}}\prod_{j=1}^{n_{\zeta}}\int{\rm
d}\bm{r}_{j}\exp(-\beta H)\prod_{\bm{r}}\delta[1-\hat{\phi}_{\rm
p}(\bm{r})-\hat{\phi}_{\rm s}(\bm{r})]\delta\left[\frac{1}{N}\int_{0}^{N}{\rm
d}s\bm{R}(s)-\bm{\xi}\right]$ (S32)

with the Hamiltonian

$\displaystyle H=\frac{\chi}{v}\int{\rm d}\bm{r}\hat{\phi}_{\rm
p}(\bm{r})\hat{\phi}_{\rm s}(\bm{r})+\frac{1}{2}\int{\rm d}{\bm{r}}{\rm
d}{\bm{r}}^{\prime}\hat{\rho}_{e}({\bm{r}})C({\bm{r}},{\bm{r}}^{\prime})\hat{\rho}_{e}({\bm{r}}^{\prime})$
(S33)

where
$\hat{\rho}_{e}({\bm{r}})=z_{+}\hat{\rho}_{+}({\bm{r}})+z_{-}\hat{\rho}_{-}({\bm{r}})+\alpha\hat{\phi}_{\rm
p}({\bm{r}})/v$ is the net charge density. Note that $z_{\pm}$ and $\alpha$
contain the sign of their charges, and we have assumed $v_{P}=v_{S}=v$.
$C({\bm{r}},{\bm{r}}^{\prime})$ is the Coulomb operator which satisfies

$\displaystyle-\nabla\cdot[\epsilon({\bm{r}})\nabla
C({\bm{r}},{\bm{r}}^{\prime})]=\delta({\bm{r}}-{\bm{r}}^{\prime})$ (S34)

where $\epsilon({\bm{r}})=kT\epsilon_{0}\epsilon_{r}({\bm{r}})/e^{2}$ is the
scaled permittivity with $\epsilon_{0}$ the vacuum permittivity and
$\epsilon_{r}({\bm{r}})$ the local dielectric constant which depends on the
local composition of the system [9].

We follow a similar procedure as Sec. 1.1, which involves: (1) identity,
Fourier and Hubbard-Stratonovich transformations to decouple the interacting
system into noninteracting polymer $\phi_{\rm p}({\bm{r}})$, solvents
$\phi_{\rm s}({\bm{r}})$, and ions in the fluctuating fields $W_{\rm
p}({\bm{r}})$, $W_{\rm s}({\bm{r}})$, $\Gamma({\bm{r}})$, $\Psi({\bm{r}})$ and
force ${\bm{F}}$. $\Psi({\bm{r}})$ is the electrostatic potential field
conjugate to the net charge; (2) the saddle-point approximation to simplify
the evaluation of the functional integral over the fluctuating fields, i. e.,
$(W_{\rm p},W_{\rm s},\Gamma,\Psi)\rightarrow(-iw_{\rm p},-iw_{\rm
s},-i\gamma,-i\psi)$; and (3) a variational treatment to account for the
fluctuation of force that controls the c.m. position of polymer chain. Then
similar self-consistent field equations as those in Eqs. S30 can be obtained
with the electrostatics relevant modification to the conjugate field of
polymer

$\displaystyle w_{\rm p}(\bm{r})=\chi\phi_{\rm
s}(\bm{r})+\gamma(\bm{r})-\frac{v}{2}\frac{\partial\epsilon({\bm{r}})}{\partial\phi_{\rm
p}({\bm{r}})}[\nabla\psi({\bm{r}})]^{2}+\alpha\psi({\bm{r}})$ (S35)

where the electric potential $\psi$ satisfies the Poisson equation:

$\displaystyle-\nabla\cdot[\epsilon({\bm{r}})\nabla\psi({\bm{r}})]=z_{+}c_{+}({\bm{r}})+z_{-}c_{-}({\bm{r}})+\frac{\alpha}{v_{0}}\phi_{\rm
p}({\bm{r}})$ (S36)

with ion concentration

$\displaystyle c_{\pm}({\bm{r}})=\lambda_{\pm}e^{-z_{\pm}\psi({\bm{r}})}$
(S37)

where $\lambda_{\pm}=e^{\beta\mu_{\pm}}/v_{K}$ is the fugacity of the salt
ion.

Finally, the resulting free energy of the system becomes

$\displaystyle W$ $\displaystyle=-\ln Q_{H}-e^{\beta\mu_{\rm s}}Q_{\rm
s}-\frac{\lambda}{2}\sum_{\kappa=x,y,z}\ln A_{\kappa}$
$\displaystyle+\frac{1}{v}\int{\rm d}\bm{r}\left\\{\chi\phi_{\rm
p}(\bm{r})\phi_{\rm s}(\bm{r})-w_{\rm p}(\bm{r})\phi_{\rm p}(\bm{r})-w_{\rm
s}(\bm{r})\phi_{\rm s}(\bm{r})+\gamma(\bm{r})\left[\phi_{\rm
p}(\bm{r})+\phi_{\rm s}(\bm{r})-1\right]\right\\}$ $\displaystyle+\int{\rm
d}{\bm{r}}\left\\{\frac{\alpha}{v}\phi_{\rm
p}({\bm{r}})\psi({\bm{r}})-\frac{\epsilon({\bm{r}})}{2}[\nabla\psi({\bm{r}})]^{2}-c_{+}({\bm{r}})-c_{-}({\bm{r}})\right\\}$
(S38)

1.3. Generalization to the Simultaneous Trap of Two Polymers

Generalization of the variational approach for trapping two polymers is
straightforward by independently fixing the c.m. positions of the two
indivdidual polymers. For simplicity, without loss of generality, we consider
two identical neutral polymers with their c.m. fixed at position
${\bm{\xi}_{1}}$ and ${\bm{\xi}_{2}}$, respectively. The semicanonical
partition function becomes

$\displaystyle Z=\prod_{\zeta=1,2}\frac{1}{v_{\rm
p}^{N}}\int{\rm\hat{D}}\\{\bm{R}_{\zeta}\\}\delta\left[\frac{1}{N}\int_{0}^{N}{\rm
d}s\bm{R}_{\zeta}(s)-{\bm{\xi}}_{\zeta}\right]\prod_{j=1}^{n_{\rm
s}}\sum_{n_{\rm s}}\frac{e^{\beta\mu_{\rm s}n_{\rm s}}}{n_{s}!v_{\rm
s}^{n_{\rm s}}}\int{\rm d}\bm{r}_{j}\exp(-\beta
H)\prod_{\bm{r}}\delta[1-\hat{\phi}_{\rm p}(\bm{r})-\hat{\phi}_{\rm
s}(\bm{r})]$ (S39)

Following the same procedure as in Sec. 1.1, the polymer density becomes

$\displaystyle\phi_{\rm p}(\bm{r})$
$\displaystyle=\sum_{\zeta=1,2}\phi_{\zeta}(\bm{r})=\sum_{\zeta=1,2}\frac{1}{Q_{\zeta}}\int{\rm
d}sq^{\dagger}_{\zeta}(\bm{r};s)q_{\zeta}(\bm{r};s)$ (S40)

where propagators $q_{\zeta}({\bm{r}};s)$ and
$q^{\dagger}_{\zeta}({\bm{r}};s)$ are subject to the same recurrence relation
in Eq. S23a and Eq. S23b as $q({\bm{r}};s)$ ($q^{\dagger}({\bm{r}};s)$) with
the average force ${\bm{f}}_{\zeta}$, variational parameter ${\bm{A}}_{\zeta}$
and c.m. position ${\bm{\xi}}_{\zeta}$. Similarly, $Q_{\zeta}$ and
${\bm{A}}_{\zeta}$ can be evaluated the same as $Q_{H}$ in Eq. S22 and
${\bm{A}}$ in Eq. S29, respectively. ${\bm{f}}_{\zeta}$ can be determined from

$\displaystyle\int{\rm
d}{\bm{r}}({\bm{r}}-{\bm{\xi}}_{\zeta})\phi_{\zeta}(\bm{r})={\bm{0}}$ (S41)

Finally, the semi-canonical free energy is given by

$\displaystyle W$ $\displaystyle=-\sum_{\zeta=1,2}\left(\ln
Q_{\zeta}+\frac{\lambda}{2}\sum_{\kappa=x,y,z}\ln
A_{\zeta,\kappa}\right)-e^{\beta\mu_{\rm s}}Q_{\rm s}$
$\displaystyle+\frac{1}{v}\int{\rm d}\bm{r}\left\\{\chi\phi_{\rm
p}(\bm{r})\phi_{\rm s}(\bm{r})-w_{\rm p}(\bm{r})\phi_{\rm p}(\bm{r})-w_{\rm
s}(\bm{r})\phi_{\rm s}(\bm{r})+\gamma(\bm{r})\left[\phi_{\rm
p}(\bm{r})+\phi_{\rm s}(\bm{r})-1\right]\right\\}$ (S42)

##  II II. The Calibration of $\lambda$

In this Section, we present the details of the calibration of $\lambda$.
Figure S1a plots $R_{\rm g}$ as a function of $\chi$ for different $\lambda$.
The reference action $L_{\rm ref}$ using a standard Gaussian corresponds to
$\lambda=1$, which largely overestimates $R_{\rm g}$ at $\chi=0.5$ as $R_{\rm
g}\approx 2.5\sqrt{Nb^{2}/6}$. This deviation indicates that it is necessary
to calibrate $\lambda$ for a quantitative prediction of $R_{\rm g}$. $R_{\rm
g}$ decreases as $\lambda$ increases. For the case of $N=200$ and $b^{3}/v=1$,
$R_{\rm g}$ exactly equals $\sqrt{Nb^{2}/6}$ at $\chi=0.5$ when
$\lambda=3.77$. Further increasing $\lambda$ leads to an underestimation of
$R_{\rm g}$. Therefore, the first step of the calibration is to determine the
$\lambda$ value at a finite chain length $N$, $\lambda(N)$, which satisfies
$R_{\rm g}=\sqrt{Nb^{2}/6}$ at $\chi=0.5$.

Figure S1: Calibration of $\lambda$. (a) $R_{\rm g}$ plotted as a function of
$\chi$ for different $\lambda$. $N=200$ and $b^{3}/v=1$. (b) The plot of
$\lambda(N)$ satisfying $R_{\rm g}=\sqrt{Nb^{2}/6}$ at $\chi=0.5$ for finite
$N$. The limiting value yields $\lambda(\infty)=3.61$.

The next step is to extrapolate $\lambda(N)$ to $N\to\infty$. As plotted in
Fig. S1b, $\lambda(N)$ has a linear relationship with $N^{-2/3}$:

$\lambda(N)=\lambda(\infty)+\alpha N^{-2/3}$ (S43)

The intercepts tells $\lambda(\infty)$ which satisfies the well-known
criterion $R_{\rm g}^{2}=Nb^{2}/6$ at $\chi=0.5$ as $N\to\infty$. $\lambda$
only depends on chain stiffness $b^{3}/v$, uniquely determined for a given
polymer regardless of the chain length. For example, $\lambda=3.61$ for
$b^{3}/v=1$ as used in our calculation.

It should be noted that, while taking $L_{\rm ref}$ to be a standard Gaussian
(i.e. $\lambda=1$) fails to accurately predict $R_{\rm g}$, it can still
predict the correct scaling relationship between $R_{\rm g}$ and $N$ in good,
$\theta$ and poor solvent regimes, as shown in Fig. S2.

Figure S2: Scaling relationship between $R_{\rm g}$ and $N$ in good, $\theta$
and poor solvents for $\lambda=1$.

##  III III. Comparison with the non-variational method using a fixed
reference of ideal coil

In this section, we compare our variational approach with Xu’s method of
constraining the middle segment via a harmonic spring [10]. As shown in Fig.
S3, our variational approach predicts the $N^{3/5}$ scaling, consistent with
the well-known result in the absence of density fluctuation. In contrast, the
application of Xu’s method to polymers in good solvent yields a wrong scaling
as $R_{\rm g}\sim N^{\sim 0.55}$. In Xu’s method, the spring constant is
derived from a fixed reference of the ideal coil conformation (where $R_{\rm
g}\sim N^{1/2}$), which fails to capture the feedback of the actual polymer
conformation. The application of this spring constant to polymers in the
swollen coil conformation leads to such deviation in the scaling exponent.

Figure S3: Comparison between our variational approach and Xu’s method of
constraining the middle segment via a harmonic spring [10]. The radius of
gyration $R_{\rm g}$ is plotted against $N$ for both methods. $b^{3}/v=1$.
$\lambda=3.61$ is used in our theory.

##  IV IV. Parameterization for Comparison in Fig. 1d

In this Section, we provide the details of the parameterization for the
comparison with experiments in Fig. 1d in the main text. Ahmed et al. studied
the kinetics of PNIPAM hydrophobic collapsing with a UV resonance Raman
spectroscopy [11]. Their results indicate that one Kuhn monomer of PNIPAM
contains $\approx 10$ chemical monomers, leading to the molecular weight of a
Kuhn monomer $M_{k}=1131.6{\rm g/mol}$. Using the single-chain molecular
weight of $M=1.3*10^{7}{\rm g/mol}$ in Wu and Wang’s work [12], the chain
length in terms of the number of Kuhn monomers is given by $N=M/M_{k}\approx
11488$. The volume of a Kuhn segment $v$ can be calculated from its molecular
weight and density ($\rho=1.1{\rm g/cm^{3}}$) as $v=M_{k}/\rho\approx
1.7\rm{nm^{3}}$.

We calculate the Kuhn length $b$ from the radius of gyration of PNIPAM at its
ideal coil state, $R_{\rm g,id}$. This state is estimated to be the position
where $R_{\rm g}$ has the most rapid change with temperature. $R_{\rm g,id}$
is estimated to be $\approx 80{\rm nm}$ from Wu and Wang’s work [12], giving
the Kuhn length $b=R_{\rm g,id}/\sqrt{N/6}\approx 1.8{\rm nm}$.

The temperature dependence of $\chi$ is obtained from fitting the experimental
data by Nakamoto et al. [13] The experimental data was measured from the
swelling equilibrium of PNIPAM gel. As shown in Fig. S4, $\chi$ is fitted to
two linear regimes as a function of inverse temperature $1/T$ ($T$ in Kelvin),
given by

$\chi(T)=\begin{cases}-5.5371*\frac{1000}{T}+18.6588&{\rm for}\
\frac{1000}{T}\leqslant 3.276\ \ (T\geqslant 32.1^{\rm o}{\rm C})\\\
-0.6528*\frac{1000}{T}+2.6566&{\rm for}\ \frac{1000}{T}>3.276\ \ (T<32.1^{\rm
o}{\rm C})\end{cases}$ (S44)

Figure S4: The fitting of $\chi$ as a function of inverse temperature into two
linear regimes. Temperature $T$ is in Kelvin. The experimental data is adopted
from the work by Nakamoto et al. [13]

## References

  * Fredrickson [2006] G. Fredrickson,  _The equilibrium theory of inhomogeneous polymers_ (Oxford University Press, 2006).
  * Grosberg and Kuznetsov [1992] A. Y. Grosberg and D. V. Kuznetsov, Macromolecules 25, 1991 (1992).
  * Liu  _et al._ [2022] L. Liu, C. Duan, and R. Wang, Polymer 258, 125312 (2022).
  * Liu  _et al._ [2024] L. Liu, C. Duan, and R. Wang, J. Chem. Phys. 160 (2024).
  * Kardar [2007] M. Kardar,  _Statistical physics of fields_ (Cambridge University Press, 2007).
  * Agrawal and Wang [2022] N. R. Agrawal and R. Wang, Phys. Rev. Lett. 129, 228001 (2022).
  * Wang [2010] Z.-G. Wang, Phys. Rev. E 81, 021501 (2010).
  * Wang and Wang [2015] R. Wang and Z.-G. Wang, J. Chem. Phys. 142, 104705 (2015).
  * Duan  _et al._ [2020] C. Duan, W. Li, and R. Wang, J. Chem. Phys. 153, 064901 (2020).
  * Xu and Wang [2021] Y. Xu and Z.-G. Wang, Macromolecules 54, 10984 (2021).
  * Ahmed  _et al._ [2009] Z. Ahmed, E. A. Gooding, K. V. Pimenov, L. Wang, and S. A. Asher, J. Phys. Chem. B 113, 4248 (2009).
  * Wu and Wang [1998] C. Wu and X. Wang, Phys. Rev. Lett. 80, 4092 (1998).
  * Nakamoto  _et al._ [1996] C. Nakamoto, T. Kitada, and E. Kato, Polym. Gels. Netw. 4, 17 (1996).