Visiting CERN for the Big Announcement

This morning I woke up rather early in order to go to the CERN campus to see the announcement of the new particle for myself. First of all:

Hell, yeah!

Higgs-like particle

Higgs-like particle

When Joe Incandela, the spokeperson for CMS, showed this plot everybody remained silent. He actually forgot to continue his speech for a few moments, sooo happy about being able to show such an insanely cool graph.

Obviously this is a great scientific discovery. Anyhow, it was not a surprise for most physicists, hence I will leave the discussion of the physics out of the way and focus on today’s event itself. (Just to be clear, this is indeed a new particle, however it is not know for sure whether it is the SM Higgs (I will be so happy if it is not (yes, click this link)))

Getting back to the event itself, the atmosphere was indeed great. The conference hall was too small for all the public to fit in, thus, in order to ensure that everybody is happy, CERN had the talks streamed to all meeting rooms on the campus. Regardless, there was a number of enthusiasts that slept in front of the doors to guarantee themselves a place in the main hall. Good for them as it seems that by already 5:30 am there were more than 150 people waiting in front of the hall (the seminar started at 9 o’clock).

After the presentations, the scientists that predicted the particle some 50 years ago were given the opportunity to say a few words. Gerald Guralnik’s comment was very appropriate:

It is wonderful to be at a physics event where there is applause like there is at at football game.

Peter Higgs had tears in his eyes when he said:

It is an incredible thing that it has happened in my lifetime.

In the spirit of these words, we should not forget all the people that were not as lucky and were not able to see the product of their work.

The CERN Director General Rolf Heuer closed the talks by saying:

As a laymen I would now say I think we have it. Do you agree?

And was greeted by ovations from the hall (and all the meeting rooms on campus). This person moderated the talks extremely well.  Indeed, it was a pleasure to watch him and the two collaborations’ spokepeople answering the journalists’ questions during the press conference that followed.

Lastly, a few places to check for more information (there is a number of articles on each blog, so check the prev/next buttons):

Gaussian optics in Sympy: studying a lens

In the last post I presented the small module for gaussian optics that I’ve created for sympy. Here I’ll try to show a real problem solved with it and with sympy.

[update]: The api is still changing. Also here I’m using ABCD matrices for something that is more easily done with conjugation relations. It’s that way so I can show how the matrices are used in the module.

The problem is finding how a laser beam is modified by passing trough a lens. What is the new waist? Where is it obtained?

Making the boilerplate

Let us first introduce the variables for the input beam and the lens. Notice that as input for the creation of the BeamParameter the Rayleigh length z_r is used and not the waist w_0 = \sqrt{\frac{z_r \lambda}{\pi}}. That’s the way for now just because it gives simpler expressions:

In [1]: from sympy.physics.gaussopt import *

In [2]:

In [3]: z_r = Symbol("z_r", positive=True)

In [4]: l = Symbol("lambda", positive=True)

In [5]: z = Symbol("z", real=True)

In [6]: f = Symbol("f", real=True)

In [7]: z_r
Out[7]: zᵣ

In [8]: l
Out[8]: λ

Defining the parameter and the lens and then calculating the parameter for the output beam:

In [9]: q_in = BeamParameter(z, z_r, l)

In [10]: tl = ThinLens(f)

In [11]:

In [12]: q_out= tl*q_in

Let’s see what does this give for the waist of the output parameter:

In [13]: w_0_out = q_out.w_0

In [14]: z_out = q_out.z

In [15]: w_0_out
Out[15]:
                     ⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽
  ⎽⎽⎽   ⎽⎽⎽⎽        ╱         1
╲╱ λ ⋅╲╱ zᵣ ⋅      ╱  ──────────────────
                  ╱              2     2
                 ╱        2⋅z   z    zᵣ
                ╱     1 - ─── + ── + ───
               ╱           f     2     2
             ╲╱                 f     f
─────────────────────────────────────────
                ⎽⎽⎽
              ╲╱ π

In [16]: z_out
Out[16]:
                        2                        2
z                      z                       zᵣ
────────────────── - ────────────────────── - ──────────────────────
           2     2     ⎛           2     2⎞     ⎛           2     2⎞
    2⋅z   z    zᵣ      ⎜    2⋅z   z    zᵣ ⎟     ⎜    2⋅z   z    zᵣ ⎟
1 - ─── + ── + ───   f⋅⎜1 - ─── + ── + ───⎟   f⋅⎜1 - ─── + ── + ───⎟
     f     2     2     ⎜     f     2     2⎟     ⎜     f     2     2⎟
          f     f      ⎝          f     f ⎠     ⎝          f     f ⎠

In [17]: together(z_out) # to simplify it a bit
Out[17]:
  ⎛       2     2⎞
f⋅⎝f⋅z - z  - zᵣ ⎠
─────────────────────
 2            2     2
f  - 2⋅f⋅z + z  + zᵣ

Those are still mostly unreadable. To simplify the expressions further I’ll introduce the dimensionless parameters Z_i = \frac{z}{f} and Z_r = \frac{z_r}{f}. Also there is the obvious substitution for the input beam waist.

In [19]: Z_in = Symbol("Z_i", real=True)

In [20]: Z_r = Symbol("Z_r", real=True)

In [21]:

In [22]: Z_out = z_out/f

In [23]: Z_out = z_out.subs({z : Z_in*f, z_r : Z_r*f})

In [24]: Z_out = together(Z_out)

In [25]: Z_out
Out[25]:
  ⎛    2          2⎞
f⋅⎝- Zᵢ  + Zᵢ - Zᵣ ⎠
────────────────────
 2            2
Zᵢ  - 2⋅Zᵢ + Zᵣ  + 1

And same for the waist:

In [26]: w_0_in = Symbol("w_0^i", positive=True)

In [27]: w_0_in
Out[27]: wⁱ₀

In [28]: w_0_in_expression = q_in.w_0

In [29]: w_0_in_expression
Out[29]:
  ⎽⎽⎽   ⎽⎽⎽⎽
╲╱ λ ⋅╲╱ zᵣ
────────────
  ⎽⎽⎽
╲╱ π

In [30]:

In [31]: w_0_out = w_0_out.subs({w_0_in_expression : w_0_in, z : Z_in*f, z_r : Z_r*f})

In [32]: w_0_out
Out[32]:
          wⁱ₀
─────────────────────────
   ⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽
  ╱   2            2
╲╱  Zᵢ  - 2⋅Zᵢ + Zᵣ  + 1

Let’s study how the waist changes. First of all, there better be no solutions to the quadratic equations in the denominator. We want it always positive in order to have real waist. Let’s extract the expression.

In [41]: quad_equ = w_0_out.args[1].args[0]

In [42]: quad_equ
Out[42]:
 2            2
Zᵢ  - 2⋅Zᵢ + Zᵣ  + 1

Ok, obviously the discriminant is -Z_r^2 < 0 and we have no real solutions. But just as a crosscheck let’s pass it trough sympy’s solve.

In [44]: solve(quad_equ, Z_in)
 Out[44]: [-ⅈ⋅│Zᵣ│ + 1, ⅈ⋅│Zᵣ│ + 1]

Great. The waist is real and positive.

Some useful results

Until now we haven’t done anything really interesting from physical point of view. We have just derived some expressions. Here are some more interesting questions.

What is the maximal possible waist obtainable with a single lens?

To find it we minimize the denominator. I lost a few minutes searching for a maximize function in sympy but didn’t find one. But obviously the maximum is at Z_i = 1 and so w_{out}^{max} = \frac{w_{in}}{Z_r}=\frac{f \lambda}{\pi w_{in}}.

It’s quite interesting that the maximum is for Z_i = 1 (i.e. z_i = f). Let’s see what does that mean for a convergent lens (it’s mostly the same for a divergent one). The beam is converging to the focal point of the lens so after refraction the beam is parallel. And we know that in order to have parallel beam (i.e. minimize divergence) we must have big waist:

In [39]: q_in.divergence.subs({z_r : waist2rayleigh(w_0_in,l)})
Out[39]:
  λ
─────
π⋅wⁱ₀

What about the minimal waist?

The predicted waist goes to zero for very big Z_i. But the gauss beam is a solution to the paraxial equation which is invalid for great curvatures of the beam front. We have the following limit:

In [41]: q_out.waist_aproximation_limit
Out[41]:
2⋅λ
───
 π

It would be interesting to know for what distance is this obtainable. The sad thing is that currently sympy can solve only very explicit equations. To rewrite it I do the following:

In [48]: waist_limit = q_out.waist_aproximation_limit

In [49]: equation = (w_0_out.args[0] / waist_limit)**2 - (1/w_0_out.args[1])**2

In [50]: equation
Out[50]:
                          2    2
    2            2       π ⋅wⁱ₀
- Zᵢ  + 2⋅Zᵢ - Zᵣ  - 1 + ───────
                             2
                          4⋅λ

In [51]:

In [52]: solve(equation, Z_in)
Out[52]:
⎡       ⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽         ⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎤
⎢      ╱       2  2    2    2         ╱       2  2    2    2 ⎥
⎢    ╲╱  - 4⋅Zᵣ ⋅λ  + π ⋅wⁱ₀        ╲╱  - 4⋅Zᵣ ⋅λ  + π ⋅wⁱ₀  ⎥
⎢1 - ─────────────────────────, 1 + ─────────────────────────⎥
⎣               2⋅λ                            2⋅λ           ⎦

There are many more interesting questions but that’s all for now.

Gaussian optics in Sympy

[update] The API shown here will most definitely change.

During my internship at LKB, Paris I’m expected to prepare a number of optical setups. I thought that it would be interesting to do the work in SymPy. That’s why I started working on the gaussopt module in sympy.physics (for the foreseeable future it won’t be upstream). I plan to make it small and basic (containing only the most basic functionality).

For the moment the only things on the blueprint are ray transfer matrices, geometric and gaussian beams.

To define a ray transfer matrix for free space or for a thin lens you just use the build in classes (a general matrix is also available):

In [1]: from sympy.physics.gaussopt import *

In [2]: distance = Symbol('distance', positive=True)

In [3]: focus = Symbol('focus', real=True)

In [4]: fs = FreeSpace(distance)

In [5]: tl = ThinLens(focus)

In [6]: fs
Out[6]:
⎡1  distance⎤
⎢           ⎥
⎣0     1    ⎦

In [7]: tl
Out[7]:
⎡  1    0⎤
⎢        ⎥
⎢  -1    ⎥
⎢─────  1⎥
⎣focus   ⎦

The module supports geometric beams but gaussian beams are more interesting, so I’ll show only them. The formalism that supports them is based on the complex beam parameter. Here is how you define a beam at its waist:

In [8]: waist = Symbol('waist', positive=True)

In [9]: wavelen = Symbol('wavelen', positive=True)

In [10]: z_rayleigh = waist**2 * pi / wavelen

In [11]: p = BeamParameter(0, z_rayleigh, wavelen)

In [12]: p.q
Out[12]:
         2
ⅈ⋅π⋅waist
──────────
wavelen

In [13]: p.w
Out[13]: waist

In [14]: p.w_0
Out[14]: waist

In [15]: p.R
Out[15]: 0

In [16]: p.z
Out[16]: 0

In [17]: p.z_r
Out[17]:
       2
π⋅waist
────────
wavelen

For the different properties shown here the names should be suggestive enough. The pretty_printer supports Greek letters for the wavelen but this is not so important here.

Now I would like to calculate how a telescopic system will act on the beam.  For example

In [18]: tl1 = ThinLens(focus/10)

In [19]: p_out = tl1*( fs*( tl*p))

The last line takes about 60 seconds and produces very UNsimplified expression. I suppose that I should do some implicit simplification. Also notice the parentheses. I should check if they are really required.

simplify(p_out) takes too much time so I’ll cheat a bit (from symbolics to numerics):

In [20]: p_out1 = p_out.subs({focus:1, distance:1*11/10, wavelen:532e-9, waist:5e-3})

In [21]: N(p_out1.q)
Out[21]: -0.110000000001061 + 1.47631233721549⋅ⅈ

In [22]: N(p_out1.w_0)
Out[22]: 0.000500000000000380

In [23]: N(p_out1.z)
Out[23]: -0.110000000001061

So with balanced telescopic system of power 10 the waist diminished 10 times and will be obtained at 11 cm after the second lens. Nice.

It would be nice to post some more realistic problems later on.