OEIS Twitterbot based on Raspberry Pi

I was actively finding fascinating afternoon-projects (as defibrillators) since the beginning of this semester, until I came across the @MoMARobot on twitter. This is a script robot spidering artworks from https://www.moma.org/  and then post them to twitter on a 6-hour basis. With the tasty Raspberry Pi 3 in my hand, I felt the urge to create my own twitter bot (on my own desk!) posting random entries from my all-time favorite website: https://oeis.org/.

The Online Encyclopedia of Integer Sequences

OEIS is an online database of integer sequences, e.g. A000796: {3, 1, 4, 1, 5, 9, 2,...}. The main purpose is to allow mathematicians or other scientists to find out if some sequence that turns up in their research has ever been seen before. Despite the fact that my love of this site is kind of an inner geek thing, it does generate a huge amount of interesting academical uses https://oeis.org/wiki/Works_Citing_OEIS.

Ever wondered the number of distinct n -carbon alkanes? No problem. Just use your high school chemistry to brutal force it for small n 's and dump it to OEIS, it will return you

A000602: Number of n-node unrooted quartic trees; number of n-carbon alkanes C(n)H(2n+2) ignoring stereoisomers.

1, 1, 1, 1, 2, 3, 5, 9, 18, 35, 75, 159, 355, 802, 1858, 4347, 10359, 24894, 60523, 148284, 366319, 910726, 2278658, 5731580, 14490245, 36797588, 93839412, 240215803, 617105614, 1590507121, 4111846763, 10660307791, 27711253769

Look at it. It's just beautiful! Cannot wait to have an hourly feed of it on twitter.

Another good news I forgot to mention is OEIS offers json data interchanging now! Do

url = "https://oeis.org/search?q="+ name +"&fmt=json"
j = requests.get(url).json()

you will get everything about a sequence!


As every major social network, twitter offers well designed APIs for status updating. Since I am a science student who cannot code properly :^, I use python package tweepy. This package contains various utilities you need to deploy a twitterbot. Simply authenticate by

auth = tweepy.OAuthHandler(consumer_key, consumer_secret)
auth.set_access_token(access_token, access_token_secret)
api = tweepy.API(auth)

and post using the api handle by

api.update_status(status = text)

Here the consumer_key, consumer_secret, access_token and access_token_secret can be obtained from the application page of you twitter account.

Schedule a Periodic Task on Raspberry Pi

Although we have got all the ingredients needed for tweeting one post, a while True loop relentlessly tweeting posts will blacklist your IP on the twitter side. However, on unix-like systems, one can schedule periodic tasks using cron. To add a cron task, execute

pi@raspberrypi:~ $ crontab -e

then use your favorite text editor to add one new line in the file

0 * * * * \abs\path\to\your\python \abs\path\to\your\twitterbot.py args

The 0 * * * * here is the cron schedule expression standing for "every hour at 0 minute". You can get assistance setting this up for your need at https://crontab.guru/.

Enjoy the Twitterbot

Yummy! The Raspberry Pi is tweeting. Follow my implementation here and get your hourly dose of integer sequences!

Notes on Vacuum Decay

This is the content I presented at the previous group meeting, and I save it here for the record. Please note that the content in this page is not fixed. The main part of this survey is based on V. Mukhanov's Physical Foundations of Cosmology.


True Vacuum, False Vacuum, Phase Transition

Vacuum Structure
Vacuum Structure

A vacuum is a minimum on a potential–state diagram. If the minimum is global, the vacuum is called a true vacuum(stable vacuum), else it is called a false vacuum(meta-stable vacuum).

A phase transition is a process that the system change from one vacuum state to another. If the state is discrete, the phase transition is called first order, else it is called continuous phase transition. First order phase transition have latent heat.

Non-Trivial Vacuum Structure

  1. If there exist both a true vacuum and a false vacuum in the potential–state diagram, the vacuum structure is called non-trivial.
  2. Not all model of electroweak interaction have the feature of non-trivial vacuum structure. Especially, the Standard Model do not have non-trivial vacuum structure. So the probe of vacuum structure is in a sense a probe of possible physics beyond SM. For instance

    V_\text{tree}(h)=\frac{1}{2}\mu^2 h^2+\frac{1}{4}\lambda h^4+\frac{\kappa}{8\Lambda^2}h^6

  3. In Quantum Field Theory considering temperature, the model’s coupling with temperature can also lead to the emergence of non-trivial vacuum structure. For instance

    V_\text{eff}(h,T)=V_\text{tree}(h)+V_1^{T=0}(h)+\Delta V_1^{T>0}(h,T)

Importance for Cosmology

  1. In the process of the Hot Big Bang Theory, the temperature keeps dropping after the Big Bang. This will lead to the emergence of non-trivial vacuum structure, therefore lead to the spontaneous symmetry breaking.

    Effective potential and Critical Temperature
  2. If the phase transition is first order, it will be quite strong and altering the energy-momentum tensor. This serves as a possible source for cosmological GW background.

What to Calculate about Vacuum Decay

  1. One of the most important parameter in many phase transitional gravitational wave model is the nucleation rate β, which represent the tunneling probability Γ from the false vacuum to the true vacuum.
  2. The main method we use to calculate the tunneling probability is Saddle Point Approximation or in a sense WKB approximation.
  3. The decay is obviously classically forbidden. So the nature of the tunneling is quantum fluctuation (instanton) & thermal fluctuation (sphaleron), so we will begin with a Quantum Mechanics situation and use analogy to go into QFT situation.

Quantum Mechanics Analogy

Potential of a 1D Particle
Potential of a 1D Particle

Decay via Instanton

The tunneling probability can be written using the propagator of quantum mechanics

P = |\langle t|f\rangle|^2 = \left|\int \mathcal{D}[x(t)]\exp(iS[x(t)]/\hbar)\right|^2

where S[x(t)] is the action of the particle defined by

S[x(t)]=\int_f^t\left(\frac{1}{2}m \left(\frac{\mathrm{d} x}{\mathrm{d} t}\right)^2 - V(x)\right) \mathrm{d} t

Doing the Wick Rotation  t \to \tau=it , we have

S[x(t)]=\int_f^t\left(-\frac{1}{2}m \left(\frac{\mathrm{d} x}{\mathrm{d} \tau}\right)^2 - V(x)\right) \mathrm{d} (i\tau)=iS_E[x(\tau)]

Where the S_E[x(\tau)] is the Euclidean action

S_E[x(\tau)]=\int_f^t\left(\frac{1}{2}m \left(\frac{\mathrm{d} x}{\mathrm{d} \tau}\right)^2 - [-V(x)]\right) \mathrm{d} \tau

while the tunneling probability can be written as

P = |\langle t|f\rangle|^2 = \left|\int \mathcal{D}[x(\tau)]\exp(-S_E[x(\tau)]/\hbar)\right|^2

Now introduce the saddle point approximation. Since \hbar is small, the main contribution of the integration comes from the extreme value of the Euclidean action \delta S_E[x(\tau)]=0 , this lead to the classical motion generated by S_E[x(\tau)] , i.e. particle moving x_\text{cl}(\tau) in a inverted potential -V(x) . This motion is called instanton, the corresponding action is called instanton action.

Thus we have

P \propto |\exp(-S_E[x_\text{cl}(\tau)]/\hbar)|^2=\exp(-S_E[x_\text{cl}^\text{cyc}(\tau)]/\hbar)

 Note that we absorb the power 2 into the instanton action and change the instanton motion into a cyclic loop, denoted by x_\text{cl}^\text{cyc}(\tau) .

Decay via Sphaleron

If the system is in a hot bath with temperature T , there is another way to decay from false vacuum to true vacuum beside quantum tunneling, that is thermal fluctuation. This way to decay adds another Boltzmann factor to the total tunnelling probability

P \propto \int \exp(-\frac{E}{kT}-S_E[x_{\text{cl},E}^\text{cyc}(\tau)]/\hbar)\ \omega(E)\mathrm{d} E

\omega(E) is the energy distribution in thermal equilibrium. The real interest of us, however, lies in the situation where the temperature $T$ is extremely high. In such case, we use saddle point approximation to calculate the dominant contribution.

The saddle point of the exponent gives

-\frac{1}{kT}-\frac{\partial}{\partial E}S_E[x_{\text{cl},E}^\text{cyc}(\tau)]/\hbar=0

 It can be easily shown (leave for Ex.) that

-\frac{\partial}{\partial E}S_E[x_{\text{cl},E}^\text{cyc}(\tau)]=\text{period of instanton with energy E}

 Thus the saddle point function gives

\frac{\hbar}{kT}=\text{period of instanton with energy E}

For extreme high temperature, the period of the instanton is approximately zero. That is the instanton motion localized around the minimum of -V(x) , i.e. the maximum of \max V(x)=V(x_m) . the system's motion x_\text{cl}(\tau)=x_m is called a sphaleron (Greek name for "ready to fall").

Thus the tunneling probability is given by substituting the action in the exponent part

P\propto \exp(-V(x_m)/T)


  • For instanton: The tunneling probability of a particle with energy $E$ is calculated by

    P \propto \exp(-S_E[x_{\text{cl},E}^\text{cyc}(\tau)]/\hbar)

     where x_{\text{cl},E}^\text{cyc}(\tau) is the periodic instanton motion with Euclidean energy -E generated by inverted potential -V(x) , which changes the potential barrier to a potential well.
  • For Sphaleron: At temperature T the tunneling probability is given by

    P \propto \int \exp(-\frac{E}{kT}-S_E[x_{\text{cl},E}^\text{cyc}(\tau)]/\hbar)\ \omega(E)\mathrm{d} E

     The dominant instanton energy E contribution to the tunneling probability at temperature T is given by

    \frac{1}{kT}=\text{period of instanton with energy E}

     At ultra-high temperature, the tunneling probability is given by

    P\propto \exp(-V(x_m)/T)

QFT Correspondence

QFT can be treated as Quantum Mechanics with infinite DOF, so the correspondence is

Screen Shot 2016-05-10 at 19.41.39

Thus we can get the corresponding EoM of Instanton and Sphaleron:

Screen Shot 2016-05-10 at 19.41.51

Equation of Motion: Bubble Solution

From the instanton and sphaleron EoM in variation form, one can obtain the EoM in its normal PDE form. If Spherical Symmetric solution is assumed (i.e. bubble solution), then h=h(r)

\begin{cases}\displaystyle\frac{\partial^2 h}{\partial r^2}-\frac{\alpha}{r}\frac{\partial h}{\partial r}=\frac{\partial V}{\partial h}\\h'(0)=0,\ h(\infty)=0\end{cases}


  • for instanton: \alpha=3 , r is the distance in Euclidean space \mathbb{E}^4
  • for instanton: \alpha=3 , r is the distance in Euclidean space \mathbb{E}^4

For different potential, one can obtain the numerical solution of the bubble profile and use the previous formula to calculate the tunneling probability. Note that the boundary condition for the ODE is strange: one condition at each side of the solution interval. This kind of ODE should be tackled by using overshoot and undershoot method.

One Numerical Solution Using scipy.optimize.odeint

One Numerical Solution Using h^6 Potential
One Numerical Solution Using h^6 Potential



I(f,A)=\int f(\vec{x})\exp\left(-\frac{1}{2}\sum_{i,j=1}^n A_{ij} x_i x_j\right)\mathrm{d}^n x.

其中 A 是半正定的实矩阵(来保证 -(1/2)\sum_{ij}A_{ij}x_i x_j 是半负定的二次型, 来让积分 I(f,A) 收敛).

但是在Wikipedia的这个页面上, 直接给出了下面这个公式:

I(f,A)=\sqrt{\frac{(2\pi)^n}{\det A}} \exp\left(\frac{1}{2}\sum_{i,j=1}^n (A^{-1})_{ij} \frac{\partial}{\partial x_i}\frac{\partial}{\partial x_j}\right)f(\vec{x})\bigg|_{\vec{x}=\vec{0}}.

前面还直接标着一个citation needed (╯‵□′)╯︵┻━┻. 这是什么鬼嘛, 给人一种钦定的感觉...

但是经过一点简单的思考, 我竟然发现可以给出一个简单的"证明"(打引号是因为这个证明显然是不严格的, 从数学系的角度来看). 下面简单的叙述一下.

首先, 矩阵 A 是实半正定的, 这意味着 A 可以正交相合对角化:

S^{T}AS=\mathrm{diag}(\lambda_1, \lambda_2, \cdots, \lambda_n).

其中 \lambda_i\geq 0 是矩阵 A 的本征值, S 是正交矩阵 SS^{T}=I .

做上面的正交变换是为了积分换元, 把积分测度变成标准的Gauss测度:

\mathrm{d} \mu_s(\vec{\rho})= \exp\left(-\frac{1}{2}\sum_{i=1}^n \lambda_i\rho_i^2\right)\mathrm{d}^n \rho.

为此, 把积分做变量代换:



|J|=\left| \frac{\partial \vec{x}}{\partial \vec{\rho}}\right|=|S|=1.

在这个变换下, 积分为:

\begin{align*} I(f,A) &=\int f(S\vec{\rho})\exp\left(-\frac{1}{2}\vec{\rho}^T S^T A S \vec{\rho}\right)|J|\mathrm{d}^n \rho \\ &=\int f(S\vec{\rho})\exp\left(-\frac{1}{2}\sum_{i=1}^n \lambda_i \rho_i^2\right)\mathrm{d}^n \rho.\end{align*}

里面用到了 S^T A S 的表达式.


I(f,A)=\int e^{-\frac{\lambda_n}{2}\rho_n^2}\mathrm{d}\rho_n \int e^{-\frac{\lambda_{n-1}}{2}\rho_{n-1}^2}\mathrm{d}\rho_{n-1} \cdots\int e^{-\frac{\lambda_1}{2}\rho_1^2}\mathrm{d}\rho_1 f(S\vec{\rho})

再假设 f 满足合适的解析条件, 就可以把 f 展开成其每个变量在0处的Taylor展开:

f( \cdots, \rho_i, \cdots)=\sum_{k=1}^\infty \frac{\rho_i^k}{k!}\frac{\partial^k f}{\partial \rho_i^k}( \cdots, \rho_i=0, \cdots)

显然带入并且交换积分和求和次序之后会出现多项式积分, 这里我们使用一个广为人知的多项式积分

\int_{-\infty} ^{+\infty} \rho^k e^{-\frac{\lambda}{2}\rho^2}\mathrm{d}\rho=\begin{cases}\sqrt{\frac{2\pi}{\lambda^{n+1}}}(k-1)!! &(k\text{ even})\\0 &(k\text{ odd})\end{cases}.

利用这个式子不停的对上面的式子进行逐次积分, 就可以得到最后的答案. 比如, 第1次积分是:

\begin{align*}I_1&=\int \mathrm{d}\rho_1 e^{-\frac{\lambda_1}{2}\rho_1^2}\sum_{k=0}^\infty\frac{\rho_1^k}{k!}\frac{\partial^k f}{\partial \rho_1^k}(0, \cdots) \\ &=\sum_{k=0}^\infty\frac{1}{k!}\frac{\partial^k f}{\partial \rho_1^k}(0, \cdots)\int \rho_1^k e^{-\frac{\lambda_1}{2}\rho_1^2}\mathrm{d}\rho_1 \\ &=\sum_{s=0}^\infty\frac{1}{(2s)!}\frac{\partial^{2s} f}{\partial \rho_1^{2s}}(0, \cdots)\sqrt{\frac{2\pi}{\lambda_1^{2s+1}}}(2s-1)!! \\ &=\sqrt{\frac{2\pi}{\lambda_1}} \sum_{s=0}^\infty \frac{1}{s!}\left[\left(\frac{1}{2\lambda_1}\frac{\partial^2}{\partial \rho_1^2}\right)^s f\right](0,\cdots)\\ &=\sqrt{\frac{2\pi}{\lambda_1}}\left[ \exp\left(\frac{1}{2\lambda_1}\frac{\partial^2}{\partial \rho_1^2}\right)f\right](0,\cdots).\end{align*}

上面的式子最后一行中的 \exp 是算子展开式 \sum_{s=0}^\infty \frac{1}{s!}\left[\left(\frac{1}{2\lambda_1}\frac{\partial^2}{\partial \rho_1^2}\right)^s f\right] 的幂级数记号. 同时注意到每次积分完成之后, 被积分过的变量在函数中置0, 没有积分过的变量仍然是自由的. 这个取值规则在下面的叙述中不再显示地写出了(不然式子就太长了, 摔!).

很显然一般情况下, f 作为 \rho 的函数是无法变量分离的, 所以在做第2次积分的时候, 上式的最后一行的 \left[ \exp\left(\frac{1}{2\lambda_1}\frac{\partial^2}{\partial \rho_1^2}\right)f\right](0,\cdots) 必须一直跟着运算. 但是如果 f 有足够好的解析性质, 那么对于不同变量的偏微分是对易的. 由Baker–Campbell–Hausdorff公式的对易情形, 可以知道

[A,B]=0\to \exp(A)\exp(B)=\exp(A+B).

所以每次积分结果一直具有很好看的形式. 比如第 r+1 次的积分结果是

\begin{align*}I_{r+1}&=\int \mathrm{d}\rho_{r+1} e^{-\frac{\lambda_{r+1}}{2}\rho_{r+1}^2}I_r \\ &=\int \mathrm{d}\rho_{r+1} e^{-\frac{\lambda_{r+1}}{2}\rho_{r+1}^2} \left(\prod_{m=1}^r\sqrt{\frac{2\pi}{\lambda_m}}\right)\left[ \exp\left(\sum_{m=1}^r\frac{1}{2\lambda_m}\frac{\partial^2}{\partial \rho_m^2}\right)f\right]\\ &=\left(\prod_{m=1}^{r+1}\sqrt{\frac{2\pi}{\lambda_m}}\right)\exp\left(\frac{1}{2\lambda_{r+1}}\frac{\partial^2}{\partial \rho_{r+1}^2}\right)\left[ \exp\left(\sum_{m=1}^r\frac{1}{2\lambda_m}\frac{\partial^2}{\partial \rho_m^2}\right)f\right]\\&=\left(\prod_{m=1}^{r+1}\sqrt{\frac{2\pi}{\lambda_m}}\right)\left[ \exp\left(\sum_{m=1}^{r+1}\frac{1}{2\lambda_m}\frac{\partial^2}{\partial \rho_m^2}\right)f\right].\end{align*}

直到, 最终:

I(f,A)=I_n=\left(\prod_{m=1}^{n}\sqrt{\frac{2\pi}{\lambda_m}}\right)\left[ \exp\left(\sum_{m=1}^{n}\frac{1}{2\lambda_m}\frac{\partial^2}{\partial \rho_m^2}\right)f\right]\bigg|_{\vec{\rho}=\vec{0}}

为了方便, 一般把微分变量从 \rho 换回 x , 这只需要计算偏微分

\frac{\partial}{\partial \rho_m}=\sum_{j=0}^n \frac{\partial x_n}{\partial \rho_m}\frac{\partial }{\partial x_n} = \sum_{j=0}^n S_{nm} \frac{\partial }{\partial x_n}.

也就是说如果把微分算子排成一列, 看成列向量, 那么

\vec{\partial_\rho}=S^T \vec{\partial_x}.


\sum_{m=1}^{n}\frac{1}{\lambda_m}\frac{\partial^2}{\partial \rho_m^2}=\vec{\partial_\rho}^T \mathrm{diag}(1/\lambda_1, 1/\lambda_2, \cdots, 1/\lambda_n)\vec{\partial_\rho}.

这时把 A 的正交相合变换表达式的逆带入, 有

\sum_{m=1}^{n}\frac{1}{\lambda_m}\frac{\partial^2}{\partial \rho_m^2}=\vec{\partial_x}^T A^{-1} \vec{\partial_x}.


\det A=\lambda_1\lambda_2\cdots\lambda_n,


I(f,A)=\sqrt{\frac{(2\pi)^n}{\det A}} \exp\left(\frac{1}{2}\sum_{i,j=1}^n (A^{-1})_{ij} \frac{\partial}{\partial x_i}\frac{\partial}{\partial x_j}\right)f(\vec{x})\bigg|_{\vec{x}=\vec{0}}.

这个等式当然可以用来方便计算函数的Gaussian积分, 但是注意到如果函数的偏导数不在某阶之后截断, 那么最后要计算的是一个无穷级数的和. 所以本质上来说计算函数的Gaussian积分还是一个很复杂的问题. 只是对于多项式情形和少数容易计算无穷级数和的情形会比较简单.

数学系的同学请不要吐槽, 我知道中间用到了很多解析性质的条件. 我写这个证明只是帮助自己理解这个等式. 诸位数学系的同学要是能帮我指出整个推导过程之中 f 要满足的充分条件是什么的话, 灰常感激~.~