import mmf_setup;mmf_setup.nbinit()
This cell adds /home/docs/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src to your path, and contains some definitions for equations and some CSS for styling the notebook. If things look a bit strange, please try the following:
- Choose "Trust Notebook" from the "File" menu.
- Re-execute this cell.
- Reload the notebook.
Dark-Bright Soliton Trains#
This notebook is to replicated some of the work done in the paper “Generation of Dark-Bright Soliton Trains in Superfluid-Superfluid Counterflow” PhysRev 106 065302
The simulations some in the paper are 3D, here will use a 1D aproximation that asumes the radia directions have a guassian profile (called tube/tube2)
Experimental setup#
\(N_{tot} \approx 450000\)
States: \(\left|1,1 \right>\) and \(\left|2,2\right>\) with scatterling Lenghts \(a_{11} = 100.40 \text{ a.u.}\), \(a_{22} \approx a_{12} = 98.98 \text{ a.u.}\)
initial state 70% \(\left|1,1 \right>\) and 30% \(\left|2,2\right>\)
Initial states look to be around 600 microns long, \(R_{TF} = 300\)
Counterflow is ramped on linearly for 1 sec (1000ms) and held constant
Trap frequencies: \(2π·1.2\) Hz axially, \(2π·174\) Hz radially in the horizontal plane, and \(2π · 120 Hz\) radially in the vertical direction. \(w_s = 2\pi(1.2,174,120)\)
\(7 ms\) of free expansion before imaging
from gpe.soc_soliton import Experiment, u
from gpe import utils
class ExperimentST(Experiment):
states = ((1, -1), (2, -2))
Rx_TF = 600.0*u.micron
t_unit = u.ms
detuning_kHz = None
detuning_E_R = 0.
B_gradient_mG_cm = 7 # Magnetic field gradients are on the order of 5-10 mG/cm
recoil_frequency_Hz = 1843.0
#recoil_frequency_Hz = 1.
rabi_frequency_E_R = 0.
trapping_frequencies_Hz = (1.2, 174, 120)
initial_imbalance = 0.7 # Nothing RF transfered, all in state[...][0]
t__expand = 9000. # Time for expansion
t__counterflow_ramp = 1000.
@property
def t__expansion(self):
return (self.t__counterflow_ramp
+ self.t__expand)
def init(self):
Experiment.init(self)
self._B_gradient = utils.get_smooth_transition(
[0.0, self.B_gradient],
durations=[0.0],
transitions=[self.t__counterflow_ramp]
)
def B_gradient_t_(self, t_):
return self._B_gradient(t_)
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[2], line 1
----> 1 from gpe.soc_soliton import Experiment, u
2 from gpe import utils
3
4 class ExperimentST(Experiment):
ModuleNotFoundError: No module named 'gpe.soc_soliton'
Initial state#
%pylab inline --no-import-all
from gpe import utils, bec2, tube2, minimize, soc_soliton
from gpe.soc_soliton import u
reload(bec2)
reload(tube2)
reload(soc_soliton)
experiment = ExperimentST(
State=soc_soliton.State2tube,
Lxyz=(1000*u.micron,),
Nxyz=(2**7,),
N=4.5e5)
s0 = experiment.get_initial_state(cool_dt_t_scale=0.00114)
#s0.constraint = 'Nab'
#s0.plot()
# Don't fix N... this should respect mu and Rx_TF
m = minimize.MinimizeState(s0, fix_N=True)
s = m.minimize()
s.plot()
%pylab is deprecated, use %matplotlib inline and import the required libraries.
Populating the interactive namespace from numpy and matplotlib
[I 20:01:48 numexpr.utils] NumExpr defaulting to 2 threads.
[I 20:01:48 root] Patching zope.interface.document.asReStructuredText to format code
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[3], line 2
1 get_ipython().run_line_magic('pylab', 'inline --no-import-all')
----> 2 from gpe import utils, bec2, tube2, minimize, soc_soliton
3 from gpe.soc_soliton import u
4
5 reload(bec2)
File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/tube2.py:11
8 import numpy as np
10 from scipy.integrate import odeint
---> 11 from scipy.stats.mstats import gmean
13 from . import bec2
14 from .bec2 import u
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/scipy/stats/__init__.py:600
1 """
2 .. _statsrefmanual:
3
(...) 595
596 """ # noqa: E501
598 from ._warnings_errors import (ConstantInputWarning, NearConstantInputWarning,
599 DegenerateDataWarning, FitError)
--> 600 from ._stats_py import *
601 from ._variation import variation
602 from .distributions import *
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/scipy/stats/_stats_py.py:7306
7302
7303
7304
7305 @xp_capabilities(marray=True)
-> 7306 @_axis_nan_policy_factory(Power_divergenceResult, paired=True, n_samples=_pd_nsamples,
7307 too_small=-1)
7308 def chisquare(f_obs, f_exp=None, ddof=0, axis=0, *, sum_check=True):
7309 """Perform Pearson's chi-squared test.
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/scipy/stats/_axis_nan_policy.py:709, in _axis_nan_policy_factory.<locals>.axis_nan_policy_decorator(hypotest_fun_in)
707 doc['Parameters'].append(_keepdims_parameter_doc)
708 doc['Notes'] += _standard_note_addition
--> 709 doc = str(doc).split("\n", 1)[1].lstrip(" \n") # remove signature
710 axis_nan_policy_wrapper.__doc__ = str(doc)
712 sig = inspect.signature(axis_nan_policy_wrapper)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/scipy/_lib/_docscrape.py:611, in FunctionDoc.__str__(self)
608 print(f"Warning: invalid role {self._role}")
609 out += f".. {roles.get(self._role, '')}:: {func_name}\n \n\n"
--> 611 out += super().__str__(func_role=self._role)
612 return out
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/scipy/_lib/_docscrape.py:566, in NumpyDocString.__str__(self, func_role)
564 out += self._str_param_list(param_list)
565 out += self._str_section("Warnings")
--> 566 out += self._str_see_also(func_role)
567 for s in ("Notes", "References", "Examples"):
568 out += self._str_section(s)
File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/scipy/_lib/_docscrape.py:517, in NumpyDocString._str_see_also(self, func_role)
515 link = f"`{func}`_"
516 links.append(link)
--> 517 link = ", ".join(links)
518 out += [link]
519 if desc:
KeyboardInterrupt:
\[\]
def get_n(mus, gs):
mua, mub = mus
ga, gb, gab = gs
G = np.array([[ga,gab],
[gab,gb]])
print(np.linalg.inv(G))
na, nb = ns = np.asarray(np.linalg.solve(G,mus))
ia = np.where(na < 0.)
na[ia] = 0.0
nb[ia] = (np.maximum(mub, 0)/gb)[ia]
ib = np.where(nb < 0.)
nb[ib] = 0
na[ib] = (np.maximum(mua, 0)/gb)[ib]
return na, nb
#mus = np.asarray([-.1+s.mu,s.mu]).reshape(s[...].shape)
#mus = s.mu
na, nb = get_n(s.mu-s.get_VHO(), s.gs)
na.shape,nb.shape
plt.plot(s.xyz[0],na)
plt.plot(s.xyz[0],nb)
#plt.semilogy(s.xyz[0],na)
#plt.semilogy(s.xyz[0],nb)
s1 = s.copy()
new_state = np.ones(s1[...].shape)
s1[...] = np.asarray([np.sqrt(0.5017)*new_state[0],
np.sqrt(0.4981)*new_state[1]])
Hy = s1.get_Hy(subtract_mu=True)[...]
na, nb = (Hy.conj()*Hy).real
plt.semilogy(s.xyz[0],na)
plt.semilogy(s.xyz[0],nb)
va,vb,vab = s.get_V()
vea,veb,veab = s.get_Vext()
na,nb = s.get_density()
plt.figure(figsize=(20,10))
plt.subplot(311)
plt.plot(s.xyz[0], va)
plt.plot(s.xyz[0], vb)
plt.subplot(312)
plt.plot(s.xyz[0], vea)
plt.plot(s.xyz[0], veb)
plt.subplot(313)
plt.plot(s.xyz[0], (va-vea))
plt.plot(s.xyz[0], s.gs[1]*nb)
#print((va-vea),(vb-veb))
na,nb = s0.get_density()
x = s0.xyz[0]
plt.plot(x,na)
#plt.plot(x,nb)
%pylab inline --no-import-all
from gpe import utils, bec2, tube2, minimize, soc_soliton
from gpe.soc_soliton import u
reload(bec2)
reload(tube2)
reload(soc_soliton)
experiment = ExperimentST(
State=soc_soliton.State2tube,
Lxyz=(1000*u.micron,),
Nxyz=(2**7,),
N=4.5e5)
s0 = experiment.get_initial_state(cool_dt_t_scale=0.00114)
s0.constraint = 'Nab'
s0.plot()
# Don't fix N... this should respect mu and Rx_TF
m = minimize.MinimizeState(s0, fix_N=True)
s = m.minimize()
s.plot()
Evolve#
from IPython.display import clear_output, display
from pytimeode.evolvers import EvolverABM
from mmfutils.contexts import NoInterrupt
#s1 = s.copy()
s1.cooling_phase = 1+0.000001j
#s1.ws[0] = 0.
#s1[...][1] = 0.
e = EvolverABM(s1, dt=.0005*s.t_scale)
fig = None
with NoInterrupt() as interrupted:
while not interrupted:
e.evolve(200)
plt.clf()
fig = e.y.plot(fig)
display(fig)
clear_output(wait=True)
Magnetic Gradient Linear Ramp#
Currently get_smooth_transitions() is used to ramp on/off different parameters. This function is not the same as the linear ramp used in the experiments so this is an attempt to find a better linear ramp.
Slope in middle is:
\[
\left(\frac{\alpha\pi}{2}\right)^{1/(2p+1)}
\]
%pylab inline
from gpe.utils import step
x = np.linspace(0, 1, 10000)
#plt.plot(x, np.vectorize(step)(x, 1, alpha=0.4))
alpha = 1.0
def rt(p, x):
return np.sign(x)*abs(x)**(1./(2*p+1))
p = 1
alpha = 1.0
y = (1 + rt(p, np.tanh(alpha*np.tan(np.pi*(2*x-1)**(2*p+1)/2))))/2
plt.plot(x, y)
plt.plot(x, x)
%pylab inline
from gpe.utils import step
x = np.linspace(0, 1, 10000)
#plt.plot(x, np.vectorize(step)(x, 1, alpha=0.4))
alpha = 1.0
def rt(p, x):
return np.sign(x)*abs(x)**(1./(2*p+1))
alpha = 2./np.pi
y = (1 + np.tanh(alpha*np.tan(np.pi*(2*x-1)/2)))/2
plt.plot(x, y)
plt.plot(x, x)
np.