Interviewer: Audrey DeVault

Biography

Dr. Joseph Falson joined Caltech’s faculty in 2019 as a visiting associate with the Department of Applied Physics and Materials Science, and became an Assistant Professor of Materials Science in 2020. His research focuses on the synthesis and characterization of quantum materials that display emergent functionalities. Dr. Falson received his B.S from the University of New South Wales (2009), M.S. from Tohoku University (2012), and Ph.D. from the University of Tokyo (2015).

Trivia

• Favorite food: Very good sushi; I did my PhD in Japan and was exposed to very good sushi by my seniors. I go very seldom, but I appreciate it.

• Favorite country: Japan. It is one of the most beautiful countries on earth. It’s got everything, it’s really quite an amazing place.

• Favorite musical genre: I enjoy all types of music, but I have recently been enjoying LA’s KJazz 88.1 radio station. It’s quite famous, and you can stream it from around the world; it’s fantastic!

• Favorite theory/effect: The Fractional Quantum Hall Effect. It is a very beautiful, complex, layered, and emergent property of condensed matter systems. It deals with the concept of a system as a whole being more than the sum of its parts. The effect only emerges in many-body systems and is very difficult to predict a priori. It is nice to find properties of materials that are difficult to predict theoretically (and be ahead of the theorists). It’s very curiosity driven.

• A fun fact most people don’t know: Most people don’t know that I’m Australian. I haven’t lived in Australia in about 10 years now, so it’s very rare that people pick up on my accent. Of course, Australians would hear me and know, but otherwise it’s difficult for people to pinpoint.

What type of research does your lab focus on?

Probing the fundamental properties of nature. Specifically, we look at crystalline materials. We grow crystals in “thin film” form, on the scale of 10 or 100 nm thick. When crystals are made very thin, their properties begin to change due to dimensionality effects. In quantum mechanics you may have seen the example of a quantum well, where the energy levels in some spatial direction become quantized. So now you have a combination of very subtle material properties that are built into the crystal structure, and you can perturb it by making the crystals into the shape that you want.

The reason why we try to make really clean crystals is because if you have crystals that are highly defective, and electrons inside them, the electrons will begin to feel the effects of the lattice much more than the effects of the other electrons. So in some situations, the lattice becomes like a secondary property, and we can begin to see the subtleties of how the electrons talk to each other while dressed in some property of the lattice. Usually these properties are very delicate, and only appear in the cleanest materials, often at very low temperatures (about 1 K).

What motivates your research?

We are motivated by the search for new paradigms. Scientists are often asked about the real-world applications of their research, but for us, it isn’t about creating a new material for use in phones or cameras. We are focused on the exploration of fundamental properties of nature.

That’s not to say our research won’t have amazing applications in the future; everything around us started out as fundamental research. Semiconductor science used to be considered uncontrolled and unimportant in comparison to other more well-understood, controllable fields like chemistry and biology. But now I am looking at a billion semiconductors on my screen. Lasers too; when the laser was invented, it wasn’t clear what its applications would be. But now our internet service is based on fiber optic cables across the seas. So just because an application is opaque at the outset doesn’t mean that it’s going to remain that way going forward.

There aren’t that many places where you can pursue pure, fundamental science like this, but Caltech is one of the few places that you can.

How did you get into Materials Science?

100% by chance, I have had a very long and winding road to get here. I was dedicated to doing medicine when I entered university but found very quickly that it wasn’t for me. I graduated with an undergraduate degree in chemistry, but at that point I was already sure that I didn’t want to do traditional wet chemistry. But I didn’t have a great idea of what I wanted to do, I just knew I was interested in electronics and crystal growth.

My 3rd and 4th year of undergrad I explored outside of my chemistry coursework and took some physics and materials science courses and found that I enjoyed semiconductor physics. When I began applying to grad school, I applied for a Japanese government scholarship that I had been introduced to in a humanities course, and won it. So I cold emailed some professors in Japan. I met with a few chemistry professors, and one directed me towards another potential mentor he though I would click with, who studied crystal growth and materials physics.

It was like falling in love in a scientific sense.

I had a very good connection with the mentor he directed me to. So I applied, got the scholarship, got into the school, and then went there for 5 years and got my PhD. I had an absolutely fantastic time. A lot of stars aligned. From there I made another international connection that led to a postdoc position in Germany for four years, and then applied to Caltech, and got the job. I was just driven by my gut feeling.

What I try to tell my students is that in life, the key to success is finding something you are able to do for 10 hours a day and not grow to hate it. If you want to be top-class successful at anything in life, it doesn’t matter how smart you are if you don’t put the time in. Especially in experimental science, so many things go wrong, and you encounter such difficult procedures, so you have to just be tough. I never get up and go to work and want to go home, and I hope that the people around me find that as well.

What have been some of the challenges and rewards of building a new lab? (How has the pandemic affected things?)

Everything has been challenging because it’s been very difficult to get extra manpower, and everything gets delayed. A lot of the equipment is now set up, but I had to do that 100% with my one student. For the vacuum equipment, usually a team would have flown in for two weeks to set it up, so it was quite a monumental task, and it’s taken us about 3 months of work. The good thing that’s come out of that is that now I know that machine like it’s the back of my hand. That’s my baby now, I know all the intricate details. So as hard as it’s been, I’m pretty happy.

What do you think is the most exciting new development in materials science?

I would say the most exciting thing isn’t necessarily a development, but more of a tide change. Early on materials science was very much semi-conductor based, because they were novel. Then the field moved on to studying a select class of very high quality semiconductors very intensely. But now we know nearly everything about them, so recently 2-D materials have been a big thing, and we are entering a stage where there is a strong appetite for a large variety of materials. The field is maturing to be more interdisciplinary, and the tide is shifting towards a search for more exotic, highly crystalline materials. The momentum has grown a lot in recent years, and there are a lot of new computational tools that help predict new materials instead of just doing blind experiments. Theorists are using machine learning and large data sets to make these predications, and it enriches my half of the problem.

What does your typical day look like? What do you do in your spare time?

On a typical day I go to lab at around 8 am, I talk with my grad student and we discuss our goals for the day. Then we try to focus on setting up equipment. I teach twice a week, and spend a lot of time planning and thinking about my lectures. My days are sometimes pretty long and difficult, and physically challenging as well, but it’s nice to see something set up at the end of the day.

In my free time I like to garden. I like it because it takes me 100% away from thinking about anything academic or sophisticated. I also like the idea that it’s like multitasking. When I’m at work, the plants are also hard at work growing.

My wife and I recently bought a house, and there are a lot of things we need to fix up, so at the moment I spend all my weekends sanding or painting. But it’s fun to be free to build stuff how you want it.

What is it like to travel the world as a scientist?

Attending international conferences and talking to people has resulted in a lot of collaboration, friendship, and support. Science is a very natural way of bridging boundaries.

Obviously, it is very different from a vacation as you are traveling to work, but I do think that its incumbent on scientists to enjoy their time while traveling. I think it’s one of the major perks of being a scientist.

I try to travel internationally as much as possible, though I will admit I am enjoying this one year where I don’t have to travel. In the year before I came to Caltech, I would fly from Europe to Caltech for a meeting, and then fly home the next day. There were some hard days and lots of 15 hours flights, so I am enjoying not travelling this year.

What do you love most about Caltech?

I like that science is the number one priority. Caltech places a lot of trust in its students and staff, including young people like me, and I enjoy the independence. The weather, and the campus are also very nice; it’s a beautiful place. There’s something special about rolling up to Caltech and going to the Athenaeum; it’s a very unique experience, very dignified.

I am very passionate about undergraduate mentoring, and I like that undergrads are strongly encouraged to participate in research.

What advice would you give to your younger self/Caltech undergrads?

Reflecting on my undergrad time, I didn’t have a designated advisor. I was very naive and was very lucky that everything worked out. I wish that I had forged a relationship with a mentor when I was younger and had someone to guide me. Luckily, I did have a mentor in graduate school, and that’s when my career really started to take shape. I would recommend you take full advantage of having an advisor. Especially going into your 20s, you think you know what’s going on, but looking back now, I had no idea. You should hunt for a good mentor, it should be someone who can advise you in research and courses, who is also in tune with your personal needs. A good mentor is a good human.

I would suggest you get involved with a lab as heavily as you can, and that you should think about lab work not as something tangential to your studies, but as your life education. It’s what they’re not going to teach you in class. I think you should take the most difficult classes possible in the first three years, research whenever you can, and then do a senior thesis. I personally didn’t do a senior thesis, and I 100% regret it. I have plans for every undergraduate in my lab to do a senior thesis. I strongly encourage it, because that’s what is going to be the most useful when applying to grad school.

Above all, though, you’ve got to have fun and challenge yourself.

Author: Patryk Kozlowski
Mentor: Garnet K. Chan and Yang Gao
Editor: Hannah X. Chen

Introduction

Due to the imminent crisis of climate change, developing sustainable methods for fuel and chemical production has become more important than ever. A significant step towards progress comes from the electrochemistry processes which use renewable energy to convert molecules from the atmosphere into higher-valued products. Catalysts are instrumental in such processes, increasing the rate, efficiency, and selectivity of the chemical transformations involved¹.

One notable example is iron, which catalyzes the conversion of atmospheric nitrogen into ammonia (used in artificial fertilizers) in the Haber-Bosch process. In fact, the 2007 Nobel Prize in Chemistry was awarded to physicist Gerhard Ertl “for his studies of chemical processes on solid surfaces” in relation to the Haber-Bosch process². Clearly, understanding the chemistry occurring at catalyst surfaces will be crucial for the future development of enhanced catalysts. In particular, we need to be able to accurately compute surface energies – quantities that govern the chemistry at catalyst surfaces.

The Catalytic Mechanism

Compared to atoms that are part of the bulk material, atoms on a surface are unstable due to their relative lack of stabilizing intermolecular attractions (Fig. 1a). This inherent instability can cause surface atoms of a catalyst to interact with nearby gas molecules (Fig. 1b). Once adsorbed to the catalyst surface, reactants undergo a desired chemical transformation and are then desorbed from the surface, at which point the job of the catalyst is complete.

Theoretical Approaches

The surface energy measures the instability of a surface relative to the bulk material. Since surface energies thereby give insight into the adsorption process, we would really like to be able to compute them accurately using theory. Though density functional theory (DFT), the current workhorse computational method of surface science, has proven inadequate to this end, wavefunction-based methods that treat correlations between electrons have shown promise as a viable alternative⁵.

Coupled-cluster theory (CC), established as the “gold standard” of quantum chemistry at the molecular level, is one such framework which has been used to study materials^6–10. CCSD methods (like coupled cluster theory including perturbative singles and doubles) have shown increased accuracy with respect to DFT in computing the properties of various materials^6,11.

However, CCSD has only recently been used to study surface chemistry, which is relevant to catalysis. In this work, I compare the performance of the simpler MP2 and more robust CCSD perturbative frameworks with that of DFT for computing energies of the (111) surface of platinum, a common catalyst.

Computational Methodology

Achieving a Gapped Platinum System

Now, since platinum is a gapless metal, its highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) are degenerate in energy. This degeneracy was problematic for the perturbative framework I wanted to use, but one can numerically introduce a twisted periodic boundary condition to artificially “open a gap” in the platinum system. Thus, in my simulations I sampled the reciprocal space of platinum’s primitive cell to identify the ideal k-point mesh yielding the largest gap. I then used this k-point mesh in all proceeding calculations.

Lattice Spacing Optimization

We optimized the bulk platinum lattice spacings computationally, yielding fair agreement with the benchmark experimental value¹⁴ (Table 1). The minimal gth-szv and the more robust gth-dzvp basis sets overestimated and underestimated the bulk platinum lattice constant, respectively. A basis set is a set of functions representing the electronic wavefunction in DFT so that the partial differential equations of the model can be implemented efficiently on a computer.

Computing Surface Energies

We then calculated platinum (111) surface energies according to Equation 1, with E_bulk and E_slab, respectively the bulk platinum and platinum surface slab energies, computed using the Table 1 bulk-optimized lattice constants. N is the size of the simulated system.

Results and Discussion

Finite Size Effects

Material surfaces are effectively infinite in size. To compute their properties, we make two approximations: the size of the k-point mesh and the number of layers perpendicular to the surface slab sampled must be finite. We then consider the effects of this finite size on computed surface energies by independently varying the number of k-points and surface slab layers sampled in calculations of platinum (111) surface energies.

As expected, we observed that surface energies steadily approach a thermodynamic limit value as both more k-points (Fig. 4ab, 6) and more slab layers (Fig. 5ab) are sampled. We see this from the observed oscillatory convergence in the former case and monotonic convergence in the latter. Furthermore, as expected, a comparison of Figures 4a and 4b shows that surface energies approached in the thermodynamic limit by the more robust gth-dzvp basis are more accurate than those approached by the minimal gth-szv basis.

DFT Performance

It is worth mentioning that we corroborated a trend observed in the literature of DFT’s systemic overestimation of surface stability with respect to experiment⁵. That is, Fig. 4ab show that DFT consistently yields lower surface energies with respect to experiment when approaching the thermodynamic limit of more k-points sampled.

MP2 Performance

The perturbative MP2 method struggled with the still near energetic degeneracy of platinum’s HOMO and LUMO discussed previously. This is because the 2^nd order energetic perturbative correction, the component of the MP2 energy aiming to capture correlations between electrons, contains the difference between molecular orbital energies ε_i, ε_j, ε_a, and ε_b in its denominator (Equation 2).

The near energetic degeneracy of the HOMO and LUMO makes this denominator go to zero, causing the MP2 energy to go to infinity in a physically unrealistic manner.

We speculate that this degeneracy caused MP2’s poor surface energies relative to DFT, evidenced by a comparison of the results from Fig. 6 and 4b, respectively, at the robust basis set (gth-dzvp) and the largest comparable reciprocal space sampling (9 k-points).

CCSD Performance

The CCSD method performed poorly in computing platinum surface energies. Even an assessment of the accuracy of CCSD surface energies was inhibited by trouble with getting CCSD calculations, which rely on an iterative method, to converge numerically. CCSD, which treats correlations between electrons based on the perturbative CC ansatz, struggled to converge numerically due to the still near energetic degeneracy of platinum’s HOMO and LUMO.

We obtained CCSD convergence only using a minimal 1 k-point sampling, though results from such a minute k-point sampling are untrustworthy, customarily yielding surface energies over 10 J/m² off the experimental value (Fig. 5, 7). Additionally, the CCSD method has a steep scaling in computational cost of O(N⁶), where N is the size of the simulated system. Thus, computations performed at the robust basis set (gth-dzvp) and modest reciprocal space samplings (16 k-points or more) required for obtaining accurate platinum surface energies proved to be prohibitively expensive.

Conclusions

This study confirmed the inadequacy of DFT for computing accurate surface energies as identified in the literature. Furthermore, we found that even for a maximally gapped platinum system, the perturbative MP2 and CCSD frameworks are unable to compute accurate surface energies. Platinum’s metallic character causes the former method to yield surface energies less accurate than DFT and plagues the latter method with numerical convergence issues. Ultimately, these results show the ineffectiveness of perturbative methods for computing surface energies of metals, like platinum.

A natural direction for future work will be to assess the performance of perturbative methods, including CCSD, for computing surface energies of non-metal catalysts, in which sufficiently large HOMO/LUMO gaps exist. We would also seek to compute other quantities that help elucidate catalytic mechanisms, such as the adsorption energy.

Further Reading

[1] A brief conceptual introduction to catalysis can be found here.

[2] An introduction to basic electronic structure theory of materials can be found in Chapters 1-4 of this reference.

[3] An introduction to coupled-cluster theory and a review of its recent application to materials science is found here.

Acknowledgements

I would like to thank Yang Gao for guiding me through my project on a daily basis in this difficult remote setting and Professor Chan for finding time to discuss my project each week throughout the summer. Lastly, I would like to sincerely thank the John Stauffer Charitable Trust for funding my SURF this summer, enabling me to get deeply involved in research I am passionate about.

References

(1)        Seh, Z. W.; Kibsgaard, J.; Dickens, C. F.; Chorkendorff, I.; Nørskov, J. K.; Jaramillo, T. F. Combining Theory and Experiment in Electrocatalysis: Insights into Materials Design. Science 2017, 355 (6321). https://doi.org/10.1126/science.aad4998.

(2)        The Nobel Prize in Chemistry 2007 https://www.nobelprize.org/prizes/chemistry/2007/press-release/ (accessed Feb 19, 2020).

(3)        What is Surface Energy? Calculation Models and More Explained https://www.ossila.com/pages/a-guide-to-surface-energy (accessed Sep 24, 2020).

(4)        Adsorption Energy & Surface Energy | Materials Square Tech Blog.

(5)        Schimka, L.; Harl, J.; Stroppa, A.; Grüneis, A.; Marsman, M.; Mittendorfer, F.; Kresse, G. Accurate Surface and Adsorption Energies from Many-Body Perturbation Theory. Nature Materials 2010, 9 (9), 741–744. https://doi.org/10.1038/nmat2806.

(6)        McClain, J.; Sun, Q.; Chan, G. K.-L.; Berkelbach, T. C. Gaussian-Based Coupled-Cluster Theory for the Ground-State and Band Structure of Solids. J. Chem. Theory Comput. 2017, 13 (3), 1209–1218. https://doi.org/10.1021/acs.jctc.7b00049.

(7)        Booth, G. H.; Grüneis, A.; Kresse, G.; Alavi, A. Towards an Exact Description of Electronic Wavefunctions in Real Solids. Nature 2013, 493 (7432), 365–370. https://doi.org/10.1038/nature11770.

(8)        Liao, K.; Grüneis, A. Communication: Finite Size Correction in Periodic Coupled Cluster Theory Calculations of Solids. J. Chem. Phys. 2016, 145 (14), 141102. https://doi.org/10.1063/1.4964307.

(9)        Gruber, T.; Liao, K.; Tsatsoulis, T.; Hummel, F.; Grüneis, A. Applying the Coupled-Cluster Ansatz to Solids and Surfaces in the Thermodynamic Limit. Phys. Rev. X 2018, 8 (2), 021043. https://doi.org/10.1103/PhysRevX.8.021043.

(10)      Wang, X.; Berkelbach, T. C. Excitons in Solids from Periodic Equation-of-Motion Coupled-Cluster Theory. J. Chem. Theory Comput. 2020, 16 (5), 3095–3103. https://doi.org/10.1021/acs.jctc.0c00101.

(11)      Gao, Y.; Sun, Q.; Yu, J. M.; Motta, M.; McClain, J.; White, A. F.; Minnich, A. J.; Chan, G. K.-L. Electronic Structure of Bulk Manganese Oxide and Nickel Oxide from Coupled Cluster Theory. Phys. Rev. B 2020, 101 (16), 165138. https://doi.org/10.1103/PhysRevB.101.165138.

(12)      Shepherd, J. J.; Grüneis, A. Many-Body Quantum Chemistry for the Electron Gas: Convergent Perturbative Theories. Phys. Rev. Lett. 2013, 110 (22), 226401. https://doi.org/10.1103/PhysRevLett.110.226401.

(13)      The first Brillouin zone of a face centered cubic lattice http://lampx.tugraz.at/~hadley/ss1/bzones/fcc.php (accessed Oct 23, 2020).

(14)      Properties: Platinum – Properties and Applications https://www.azom.com/properties.aspx?ArticleID=601 (accessed Jul 8, 2020).

(15)      Boer, F. R.; Boom, R.; Mattens, W. C. M.; Miedema, A. R.; Niessen, A. K. Cohesion in Metals: Transition Metal Alloys; North-Holland, 1988.

Author: Lucca S. de Mello^†
Mentor: Charles L. Steinhardt^‡
Editor: Laura Lewis

^†ldemello@caltech.edu, Caltech
^‡steinhardt@nbi.ku.dk, Cosmic Dawn Center

Abstract: In 2007, astronomers discovered highly energetic, millisecond-long radio waves of extragalactic origin. These Fast Radio Bursts (FRBs) are poorly understood, and there is no definitive explanation for their origin yet. Using a novel machine learning technique, we have made significant progress toward determining whether FRBs can be separated into discrete subclasses, or if all FRBs are instances of the same astrophysical phenomenon. The technique consists of modeling the entire intensity-versus-time curve of these bursts as vectors in a high-dimensional space and clustering these vectors in two dimensions. We expect that grouping FRBs by their physical properties will also group them by their physical origin.

Background

With advances in telescope technology, humanity has recently acquired the ability to observe the sky more
closely, noticing fascinating changes that occur within the span of a millisecond. In 2007, a new class of astronomical object was discovered: fast radio bursts, or FRBs (Fig. 1). FRBs are highly energetic¹ pulses in the radio spectrum that last for approximately two milliseconds at most. Due to their very recent discovery, little is known about them. Thus, the purpose of this project is to explore this field by further investigating the origin and properties of FRBs.

The high energy and short duration of FRBs make them an especially intriguing class of astronomical object. The special theory of relativity asserts that causality travels at the speed of light. So, in the span of a millisecond, causality travels approximately 300 km. This implies that the diameter of the source of any signal with a duration of one or two milliseconds – the typical duration of an FRB – is on the order of 300 km or smaller, otherwise a complex (and thus unlikely) synchronization mechanism would be required for the bursts to be so short. Relative to distances on an astronomical scale, the order of hundreds of kilometers is fairly small.

It is also known that this mysterious source, although small, must be associated with some extremely energetic astrophysical process¹. We know this because, in spite of their short duration and the cosmological distances they travel before reaching us, FRBs are as bright in the radio frequencies as a galaxy when viewed from Earth [1].

FRB science is a nascent but thriving field of astronomy. Although the dataset is still limited, the host of instruments scattered throughout the world listening for new FRBs will give us an accelerating influx of data, which can help us further understand the nature of FRBs. We have already observed that they originate from the entire sky, not just from directions coplanar with the Milky Way³, implying that most of them are of extragalactic origin⁴. Additionally, likely candidates for what could be the source of FRBs have been identified; these include neutron stars [2] and black holes. However, due to the very recent discovery of FRBs, there is currently no definitive explanation for what astronomical object is the source of FRBs, and we are currently far away from truly understanding what they are.

FRBs are not the only kind of short-duration, high-energy burst of radiation that has been observed. There are also the gamma-ray bursts (GRBs), which are longer and brighter than FRBs. GRBs have been studied for far longer – since the late 1960s – so they are better understood. Indeed, one important milestone in GRB study has been their separation into two discrete classes, “short” and “long”⁵. The realization that not all GRBs are the same was important because, by studying these classes separately, we were able to connect each of them with an astrophysical origin: we now know that “long” GRBs are generated from the collapse of massive stars, while “short” GRBs come from neutron stars⁶. This naturally raises the question of whether or not the situation for FRBs is similar to that of GRBs. Are all fast radio bursts the same?

Goals

Despite the important advancements that have been made in the process of understanding FRBs, the key step of classification has not yet been realized. As a result, it has not been possible to connect FRBs to an astrophysical origin – or origins. So the problem we chose to tackle in this research project was that of classifying FRBs.

To this end, we have applied a novel machine-learning technique known as “dimensionality reduction,” which has recently proven successful in solving the long-standing classification problem for gamma-ray bursts (GRBs)⁷. It was known that gamma-ray bursts appear to have two types – called “short” and “long”, as previously discussed – but it had not been possible to unambiguously separate these two types of objects. For example, some “short” GRBs are actually longer in duration than some “long” GRBs.

Earlier attempts at classifying GRBs took the light curve – the curve of intensity [3] versus time, integrated over all frequency channels – and make a small number of summary statistics (such as the duration of the GRBs) which lend themselves easily to physical interpretation. However, this always resulted in some ambiguous objects, which might have been either short or long, as Fig. 2 exemplifies.

To solve this problem, instead of only using two measurements, all of the data were use. The challenge with this is what is known as “the curse of dimensionality”⁹. Suppose we want to plot objects in N dimensions, where N is the number of measurements. As N increases, so does the sparsity of the data points, meaning that all objects appear to be dissimilar in many ways. The solution was to use what is known as a dimensionality-reduction algorithm. As the name suggests, dimensionality-reduction algorithms take in a set of high-dimensional vectors that model some dataset, then cluster those vectors in a lower-dimensional space – usually in two dimensions – grouping them in a way that preserves the structure of the high-dimensional data. This means that, in the low-dimensional output space, similar objects stay close in and dissimilar objects stay distant. This idea is illustrated in Fig. 3.

The result was Fig. 4, which shows two discrete categories of GRBs, with every object unambiguously classified as either short (orange) or long (purple). Since the problem of classifying GRBs appears to be similar to the problem of classifying FRBs, the specific goal of this project was to make progress on the latter by applying a promising technique that worked on the former. That is, our ultimate goal was to get something similar to Fig. 4, but with FRBs instead of GRBs [4]. Failure in doing so would indicate that either:

1. FRBs are, in essence, all instances of the same astrophysical phenomenon.
2. FRB subclasses do exist but may not be easily revealed by dimensionality reduction alone.
3. Unlike GRBs, FRBs can only be classified in a continuous spectrum.

Methods

As previously discussed, clustering the FRBs by representing each burst as a collection of summary statistics has not been a successful approach. So, instead, our technique consists of clustering the FRBs by the “similarity” of their light curves. For example, some FRBs have light curves that end just as abruptly as they start, while others have an exponential fallout tail after their maximum, like the burst “rings” for a period of time after peaking. Another component of this intuitive notion of similarity is that we are only concerned with the shape of the light curve; that is, if a light curve is simply a scaled-up version of another, we consider them to be similar. This is because the scaled-up (but otherwise similarly-shaped) light curve is simply a more intense version of the latter curve. Intensity should not be a factor because FRBs that appear to be more intense may simply be closer to us, so letting the bursts’ intensity affect their classification could lead to the FRBs being grouped by their sources’ proximity to Earth. Proximity is neither an intrinsic property of the FRB nor of its source, so such a clustering would be physically meaningless. For more details on the similarity of FRB light curves, see Fig. 5.

Quantitatively, we can encode the similarity of these light curves by processing the light curves into vectors in a way that:

1. Uses the raw light-curve data instead of summary statistics, making no assumptions on how strongly each physical measurement of the FRBs should influence the classification. Consequently, the clustering process is dictated by the raw data and nothing else.
2. Encodes the entire light curve. In practice, this means that it must be possible to reconstruct the light curve of an FRB from the vector representing it. In other words, the process of transforming FRBs into vectors must be reversible. As a result, our approach has perfect information on the light curves being clustered.
3. Ensures that the distance between two vectors is indeed a measure of dissimilarity in topology, or shape, of the corresponding light curves.

Since each vector stores all the information we need to reconstruct the corresponding light curve, the vectors are in a very high-dimensional space. Vectors are something we can naturally cluster, so it would stand to reason that the next step is to simply apply any clustering algorithm to the result. However, as discussed in the previous section (Goals), this would not give satisfactory results because of the “curse of dimensionality,” so we must apply a dimensionality-reduction algorithm.

The dimensionality-reduction algorithm we chose is the t-distributed stochastic neighbor embedding, or t-SNE. t-SNE assigns each high-dimensional data point to a temporary position in the low-dimensional space, called a “map point” [5], and gradually updates these positions in a way that attracts similar objects to each other and repels the objects that are dissimilar.

One physical analogy for how this is accomplished is as follows. First, the map points are connected to each other by springs. The smaller the distance between two high-dimensional data points, the stiffer the spring connecting their map points. As a result, map points are attracted if they are far apart but their corresponding objects are close. If two map points are close while their corresponding objects are dissimilar, they are repelled. Then, this system evolves until it reaches equilibrium. The equilibrium state associates each data point to its position in the low-dimensional space. For an illustration of what would result from this, refer again to Fig. 3 as the dimensionality-reduction algorithm used to compute that example was, in fact, t-SNE.

The example of classifying hand-drawn digits also illustrates how t-SNE preserves the macroscopic structure of the objects being clustered: fours look like nines, so those two clusters end up together; nines look like sevens, who
look like ones, etc. This structure-preserving property is essential for our goal of clustering FRB light curves by topological similarity, which is why we chose t-SNE for this project. Additionally, the previously-discussed success in unambiguously clustering GRBs (Fig. 4) was also accomplished using t-SNE, further motivating our choice of dimensionality-reduction algorithm.

Results

After fitting t-SNE to the processed, high-dimensional data and calibrating the algorithm’s parameters, we have obtained Fig. 6 as our program’s final output. As our figure illustrates, we have separated twenty-one FRBs into two groups: one group containing seven FRBs and the other containing fourteen. Recall that our goal was to obtain something similar to Fig. 4, but with FRBs instead of GRBs. In Fig. 4, the physical interpretation of the two subclasses is known: the orange cluster represents “short” GRBs, and the purple cluster represents “long” GRBs. Indeed, the short/long separation of GRBs was known before t-SNE was applied to GRBs (Fig. 2). Unlike GRBs, however, FRBs are still too shrouded in mystery for us to make a probable physical interpretation of the clusters in our result.

There is a significant variation in the measured intensity of FRBs. As previously discussed, while an FRB’s higher measured intensity may be an intrinsic property of the burst, it may also result from that FRB having originated closer to Earth. So, to ensure that our program is not biased by FRB intensity – that is, to show that we are not clustering FRBs by their proximity to Earth – we have also produced Fig. 7-8, in which FRBs are colored by peak intensity.

If our algorithm were biased by intensity, the resulting plot would have been, at least approximately, a spectrum from the “low-intensity” colors to the “high-intensity” colors. Instead, these figures show a mostly diverse distribution of FRBs, with both clusters containing bursts in the mid- and low- intensity ranges. The very few bursts that are in the high-intensity range are in the same cluster, but the previous observation – coupled with the fact that the light curves were normalized by peak intensity – leads us to believe that this is most likely not indicative of bias. Of course, a larger sample size is necessary for any definitive assertions on bias, or lack thereof, to be made.

Our program and the clustering we have obtained from it are both preliminary, but they are the first of their kind. Crucially, the quality of the classifications it produces will only increase with the amount of data it is given. With the global community of radio astronomers currently listening for FRBs, the amount of available FRB data is certain to grow in an accelerating fashion. We believe that, if our program is adapted to make full use of this progressively larger volume of data, it will be a significantly useful tool in solving the puzzle of FRBs.

Conclusion

We have developed a program that, given a standardized set of raw FRB data, clusters the FRBs by the similarity of their light-curves. Our program is dictated by the raw data alone, making no assumptions on how strongly each physical measurement of the FRBs should influence the classification. It has perfect information on each light curve being clustered, and everything indicates that it is not biased by each source’s proximity to Earth. Our program employs t-SNE, an algorithm that was successful in a problem very similar to that of classifying FRBs. The motivation for doing so is that we expect that grouping FRBs by their physical properties will also group them by their physical origin.

Our preliminary results imply that our program for processing and grouping FRB light-curves is a promising candidate as a solution to the challenge of classifying FRBs. Further, it implies that more definitive classification results obtained by running our program on a larger dataset will likely be unbiased with respect to proximity.

While our result does represent a necessary first step in the challenge of classifying FRBs, the problem has not
been solved. Due to several technical factors related to the very recent discovery of FRBs, we only had access to
a small sample size. One of the most difficult parts of studying FRBs is that the raw data is not easily accessible, as
there is no publicly-available catalog of raw FRB data yet [6]. Additionally, raw FRB data is stored in a large variety of file formats; if there were a standard for storage of raw FRB data, classification studies would probably advance much faster. For these two reasons, conclusive further research in this area will likely required the creation of a catalog of raw FRB data that uses a standardized file format [7].

Moreover, the algorithm developed in this project can only produce meaningful results if the dataset uses the same sampling period across all FRB recordings. This means that most datasets are not compatible with each other; some datasets, such as the one used in this project, are not even consistent within themselves. And, as previously discussed, FRB data is still limited. So, to obtain more definitive classification results through our algorithm, further research may have to wait until there is a larger collection of FRB observations with practically identical time resolutions.

Are all fast radio bursts the same? We do not know, and it may still take time until it is possible to answer this
question. But we have developed a tool that can be very useful in doing so; all we need is more data. We hope that,
with a large, standardized, and easily accessible catalog of raw FRB data, the use of machine learning in the classification of FRBs will prove to be a powerful way in which the nature of these mysterious objects can be better understood.

Acknowledgements

This research project is a part of, and was made possible by, the SURF program. I thank the Caltech Student-Faculty Office for their support, and I thank my mentor, Charles L. Steinhardt, for his help in executing this project. I think Duncan Lorimer for helping us obtain raw FRB data.

Appendix

The data and code used in this research can be found in this Github repository.

For this project, twenty-seven plaintext files of raw FRB light-curve data were acquired with the help of Dr. Duncan Lorimer. So that the light curves may be meaningfully compared, the twenty-one files (of those twenty-seven) whose data used the mode sampling period (0.97 ms) were selected. These 21 FRBs are in the “data” directory of the linked repository. The other six are in the “bad data” directory, which is not used by the program. The table below shows which 27 FRBs were used in this project, ordered by when they occurred [8]. The first 21 FRBs in the table were recorded with practically the same sampling period, and therefore can be meaningfully compared using our program. Accordingly, these twenty-one FRBs are the ones in Fig. 6-8.

FRB	Sampling Period
FRB010125	9.70000E-04
FRB010621	9.70000E-04
FRB090625	9.70000E-04
FRB110214	9.76562E-04
FRB110220	9.76562E-04
FRB110626	9.76562E-04
FRB110703	9.76562E-04
FRB120127	9.76562E-04
FRB121002	9.76562E-04
FRB130626	9.76562E-04
FRB130628	9.76562E-04
FRB130729	9.76562E-04
FRB131104	9.76562E-04
FRB140514	9.76562E-04
FRB150215	9.76562E-04
FRB150418	9.76562E-04
FRB150610	9.76562E-04
FRB150807	9.76562E-04
FRB151206	9.76562E-04
FRB151230	9.76562E-04
FRB160102	9.76562E-04
FRB171209	7.81221E-05
FRB180301	1.17000E-03
FRB180309	1.25016E-04
FRB180311	1.51286E-03
FRB180714	1.39082E-03
FRB180923	4.88281E-03

References

1. L. Billings. A Brilliant Flash, Then Nothing: New Fast Radio Bursts Mystify Astronomers (2013).
2. D. Lorimer. A bright millisecond radio burst of extragalactic origin (2007).
3. M. Caleb et al. The first interferometric detections of Fast Radio Bursts (2017).
4. M. Kiyoshi et al. Dense magnetized plasma associated with a fast radio burst (2015).
5. C. Kouveliotou et al. Identification of Two Classes of Gamma-Ray Bursts (1993).
6. E. Nakar. Short-Hard Gamma-Ray Bursts (2007).
7. C. Jespersen et al. An Unambiguous Separation of Gamma-Ray Bursts into Two Classes from Prompt Emission Alone (2020).
8. R. M. Shannon et al. The dispersion-brightness relation for fast radio bursts from a wide-field survey (2018).
9. R. E. Bellman. Adaptive control processes: a guided tour (1961).
10. D. Michilli et al. An extreme magneto-ionic environment associated with the fast radio burst source FRB 121102 (2018).
11. C. D. Bochenek et al. STARE2: Detecting Fast Radio Bursts in the Milky Way (2020).
12. J. P. U. Fynbo et al. The optical afterglow of the short gamma-ray burst GRB 050709 (2005).
13. A. Shahmoradi. Hardness as a spectral peak estimator for gamma-ray bursts (2010).
14. Pedregosa et al. Scikit-learn: Machine Learning in Python (2011).

Footnotes

[1] Fig. 5 of On the Normalised FRB Luminosity Function shows the luminosity function of FRBs, i.e. the number of FRBs per luminosity interval. Comparing it to the luminosity function of galaxies shows that, in general, the luminosity of FRBs is indeed similar to that of galaxies.

[2] Some FRBs are polarized, which indicates that at least those FRBs were emitted from a specific kind of neutron star called a magnetar¹⁰. Indeed, the first FRB to have been linked to a known source, FRB200428, was found to originate from a magnetar¹¹.

[3] By “intensity,” we mean Spectral flux density, measured in Janskys (Jy).

[4] Though we did not assume that the number of FRB subclasses, if they exist, is only two.

[5] Map points are computed using probabilistic methods beyond the scope of this paper.

[6] There is frbcat.org, a publicly-available catalog of FRB data, but it does not store raw FRB data; it is limited to summary statistics. While these are useful in general, using summary statistics has not been a fruitful strategy in the challenge of classifying FRBs.

[7] Or, alternatively, the extension of frbcat.org into one such catalog.

[8] Not necessarily when they were discovered; many FRBs, including the first one ever found, are discovered from old archival data.

Author: Gabrielle Dituri
Mentors: Maria Spiropulu, Christina Wang, Cristian Pena, and Si Xie
Editor: Stephanie Chen

The God Particle and Why it is Troublesome

The Higgs boson (nicknamed the “God Particle”) was discovered in 2012 at the Large Hadron Collider (LHC) by ATLAS and CMS¹. This nickname arrived due to a hilarious publishing decision in which the author of a book about the Higgs particle decided to shorten the original title (“the Goddamn Particle”) into a title that would appeal to a wider range of audiences (“the God Particle”). The original title was named to comedically highlight the difficulty in finding proof of the particle’s existence, which took about 50 years and an expensive particle accelerator. Despite this decision to shorten the name, the God Particle nickname accurately emphasizes the importance of the Higgs discovery to modern physics. The discovery of the particle provided credence to the Standard Model (SM), which is a model describing all the elementary particles and their interactions. The Higgs is an elementary particle previously predicted by the SM, but never before seen until 2012. The discovery of the particle confirmed that the Standard Model was able to correctly predict the existence of a particle before any data was collected. This is extremely important to particle physics as this shows that unknown particles can be predicted before we have the technology to find physical proof of them. According to the Standard Model, the Higgs can decay into fermions and massive gauge bosons (Fig. 1).

The SM dictates that the Higgs boson decays into SM particles including fermions, which are elementary particles with half odd integer spin (1/2, 3/2, etc.) that obey the Pauli Exclusion Principle.  However, we are looking to provide credence to an alternative theory³ which  suggests that the Higgs decays first into a pair of long lived particles (LLPs). These particles then decay into SM particles which are measured in the Cathode Strip Chamber (CSC), the endcap Muon System in CMS. The Higgs has a 20% probability of decaying into invisible particles, such as neutrinos or beyond the standard model (BSM) particles. Therefore, there is a possibility that the Higgs can decay to LLPs, which then decay back to decay modes that can be detected, like fermions⁴. LLPs are BSM particles, meaning they are extremely important to study as they could help solve many current problems in physics, such as the matter/antimatter asymmetry problem⁵. My research aims to verify the theory that the Higgs decays into LLPs by looking for this signature of the LLPs decaying back to SM particles in the CSC system. The decayed LLPs generally create a large cluster of RecHits (reconstructed hits in the detector) in the CSC thus allowing us to accurately study them.  Furthermore, we would like to understand if the number of RecHits are related to mass and the decay mode of the LLPs, and how the distribution is different from that of background. The background mostly consists of particles such as pions and kaons that are produced at the interaction point (the point where the protons collide) instead of particles produced from the displaced decays. 

Methods and Results

For this study, we utilized 10 simulated samples (Fig. 3, col. 1). Here, k is a kaon, pi is a pion, d is down quark, τ is tau lepton, e is an electron, and b is bottom quark. These particles were used as they represent a broad spectrum of particles, including some charged/neutral, quarks/leptons, etc.

Acceptance vs Lifetime

We were able to plot a graph of acceptance (the number of events that have at least 1 LLP decay in CSC) vs proper lifetime (Fig. 2). The larger the number of signal events, the more sensitive the detector is to a decay. From this graph, we could clearly see that the most sensitive lifetime is highly dependent on the LLP mass and that the Muon System is sensitive to low mass LLP with shorter lifetime.

Acceptance and Clustering Efficiency

We then calculated the acceptance for events pertaining to LLPs with differing lifetimes (Fig. 3). We also looked at how many of those events had a cluster matched to LLP, where a cluster has 50 or more RecHits. The cluster matched simply means that a cluster was produced from a specific LLP.  By analyzing these two percentages, we were able to conclude the Clustering Efficiency in RecHitCluster. The Clustering Efficiency is defined as the number of events with a cluster matched to LLP/ number of events that have at least 1 LLP decay in CSC. This is important as it tells us if the formation of clusters is related to the LLP mass and the type of particles that they decay to.

Decay Particles, ms (GeV), cτ (mm)	% Acceptance	% Events with a Cluster Matched to LLP	Clustering Efficiency (%)
ee, 0.5, 500	3.04	0.90	29.50
ee, 1, 100	15.65	5.00	31.94
ee, 1, 500	6.59	1.95	29.67
pi0pi0, 1, 100	15.65	5.00	31.98
pi0pi0, 1, 500	6.58	1.99	30.18
pi+pi–, 1, 500	6.55	4.05	61.90
k+k–, 1.5, 500	8.80	5.70	64.83
dd, 1, 1000	3.74	1.96	52.43
ττ, 7, 1000	14.15	6.45	45.58
bb, 15, 1000	16.35	11.94	73.04

Relationship Between Number of Events and Mass

We wanted to find the relationship between the number of events/hits and the lifetime and mass of the LLP the particles decayed from. This gave us insight into how sensitive the detector is concerning these factors. One of the ways we did this was by analyzing the relationship between events and X Spread (how much the clusters are spread in the X direction of the detector) , Y Spread (how much the clusters are spread in the Y direction of the detector), N Station (station in the detector where the hits are detected), Cluster Size (number of RecHits in the clusters), pseudorapidity of the cluster (represents the angle of the particle relative to the z axis), and Avg Station (average of the station number of the RecHits, 4 stations for the CSC).

• 

• 

• 

• 

• 

• 

Based on the table (Fig. 3) and the graphs above, the initial conclusion was that the detector is sensitive to the branching ratio (the ratio of Higgs decaying to LLP) of similar order of magnitude for LLPs with different masses. We concluded this because the Clustering Efficiency and the RecHits distribution are similar for the different signal samples. Additionally, we can see that the detector is sensitive to decay particles since pi⁺pi^– and pi⁰pi⁰ have extremely different Clustering Efficiencies despite being produced from LLPs with the same mass and lifetime.

Clustering Efficiency and Geometry

We plotted the Clustering Efficiency with respect to several factors, including the LLP decay position and the energy of the LLP. This shows us how the Clustering Efficiency varies in the detector and which showers penetrate the stations (Fig. 10-12).

• 

• 

• 

The LLP decay vertex Z/R plots clearly show the shape of the detector: when the LLP decays in the sensitive detector volume the Clustering Efficiency increases. Fig. 13 directly demonstrates why the Clustering Efficiency fluctuates – the iron/steel in the detector prevents a high Clustering Efficiency while the sensitive detector region encourages a high Efficiency. There also seems to be a small increase of Clustering Efficiency with LLP energy, but it is not highly dependent. As all these graphs exhibit the same general shape, we concluded that the detector is sensitive to different particles. More specifically, we can see from Fig. 10-12 that Clustering Efficiency is dependent on the type of shower – a particle shower is a cascade of particles produced from a particle collision. Hadronic shower particles like bb tend to have a higher Clustering Efficiency, and electromagnetic (EM) shower particles like ee tend to have a lower Clustering Efficiency. Because the Clustering Efficiency is higher for hadronic shower particles, we know that the level of sensitivity is different depending on the type of shower.

2D Clustering Efficiency

Finally, we wanted to analyze the 2D Clustering Efficiency plots, which would help to further show the difference between hadronic and EM showers. We knew from the previous graphs that the Clustering Efficiency changed according to the type of shower due to the geometry of the detector (Fig. 13). However, we wanted to visualize the dependence on detector geometry in both the R and Z direction of the detector. Below are the 2D Clustering Efficiency graphs, with the red boxes representing the detector volume.

• 

• 

• 

• 

In a Nutshell (more like Particle Detector)

In summary, we saw that the Muon System is sensitive to all LLP decay modes – ee, bb, dd, ττ, pi⁰pi⁰, pi⁺pi^–, k⁺k^–, but that the level of sensitivity is different between these decay modes. Also, our preliminary study shows that Clustering Efficiency and Cluster Size are correlated with the LLP energy and momentum. Furthermore, we found that hadronic showers have more penetration than EM showers. Most importantly, we found that the Muon System is sensitive to LLPs with mass as low as 0.4 GeV; this is the first study at the CMS to set limits on the sensitivity of low mass LLPs.

Why Does This (anti)Matter?

Dark matter may be comprised of these BSM particles which is why this field is important to study: dark matter makes up most of the universe, and we still do not know exactly what it is. In the future, I hope to continue my research on this topic and find the exact level of sensitivity to LLP mass and decay mode. I additionally hope to research the lowest mass that the detector is sensitive to, which would help in looking for BSM particles. In general, conducting studies like this one, not bound by the SM, brings us physicists closer than ever to finding a unifying theory to describe the world around us.

Acknowledgements

Thank you so much to my SURF supervisor, Maria Spiropulu. Without her mentorship, I never would have undertaken such an amazing project. Also thank you to Cristian Pena and Si Xie who constantly reviewed my work and gave me great feedback. Special thanks to Christina Wang who guided me through every aspect of the project and helped me realize my passion for this field!

References

1. Gray, Heather, and Bruno Mansoulié, “The Higgs Boson: the Hunt, the Discovery, the Study and Some Future Perspectives.” ATLAS Experiment at CERN, 4 July 2018.

2. Particles of the Standard Model of particle physics (Image: Daniel Dominguez/ CERN).

3. Craig, Nathaniel, et al. “Naturalness in the Dark at the LHC.” Journal of High Energy Physics, vol. 2, 23 Mar. 2015, doi:10.1007/jhep07(2015)105.

4. Blas, J. De, et al. “Higgs Boson Studies at Future Particle Colliders.” Journal of High Energy Physics, vol. 1, 9 May 2019, doi:10.1007/jhep01(2020)139.

5. The Matter-Antimatter Asymmetry Problem. CERN, home.cern/science/physics/matter-antimatter-asymmetry-problem.

Author: Annabel Gomez
Mentor: Dr. Xiaoqing Pi
Editor: Suchitra Dara

Introduction

Solar activities, such as solar flares and coronal mass ejections, produce perturbations in solar wind and the interplanetary magnetic field (IMF). These perturbations interact with the Earth’s magnetosphere, resulting in geomagnetic storms, charged particle precipitation, and plasma convection in the high-latitude ionosphere. These events can produce various ionospheric disturbances depending on location, geomagnetic conditions, presence of ionospheric currents, plasma convection, etc. The ionospheric disturbances, including the process of plasma instabilities, give rise to irregular structures in ionospheric density distribution, also known as ionospheric irregularities. 

The Wide Area Augmentation System (WAAS) was developed by the Federal Aviation Administration (FAA) of the United States to support aviation navigation [1]. WAAS provides a ground-based global positioning system (GPS) tracking network, master stations for data processing, and geostationary satellites to relay WAAS data. To increase navigation accuracy, WAAS generates ionospheric correction data and broadcasts it to its users. Ionospheric irregularities cause scintillation of GPS signal power and phase, which can interrupt radio signal reception and degrade the integrity and accuracy of technology applications that rely on the Global Navigation Satellite System (GNSS).

Ionospheric irregularities can be explained and measured as described by Figure 1. The progression of rate of total electron content (TEC) is measured as the rate of TEC (ROT). Both ROT and rate of TEC index (ROTI) are measurements of ionospheric irregularities that can cause L-band radio signal scintillations; the ROTI index is the standard deviation of ROT [2][3]. A value of ROTI at or near zero represents a low level of disturbances, while a high value represents the presence of ionospheric irregularities indicating a stormy period. 

The goal of this research is to investigate which space weather and geophysical measurements are useful in machine-learning (ML) techniques for a prediction system that uses historical GPS data to predict the location, time, and intensity of concerning irregularities. The system is targeted to reduce the impact of scintillation effects on GNSS data and applications. 

Correlation Analyses

The ionospheric data used for this analysis was ROTI data with 5-minute intervals. The ROTI data was obtained by processing GPS phase data collected at Fairbanks, Alaska (station name: FAIR, station coordinates: 64.98º geographic latitude and 212.5º geographic longitude). At each epoch, <ROTI> is obtained by averaging ROTI values over a number of satellites in view. 

The first main parameter used was auroral electrojet ((AE)) indices derived from longitudinally distributed magnetic field data in the polar region. This data includes four indices: 

• AU: The upper envelope expressing the strongest intensity of eastward auroral electrojets.
• AL: The lower envelope expressing the strongest intensity of westward auroral electrojets.
• AE:
• AO:

A correlation analysis was first conducted between the ROTI and (AE) data. On average, the <AE> index, averaged over 5-minute intervals, had a higher correlation with <ROTI>, though <AU> occasionally correlated more with <ROTI> such as the case for March 01, 2015 (Figure 2). In general, there were more complications with <AL> and <AO> indices, as is evident with the negative values in contrast to positive <ROTI> values. Additionally, while <AE> did correlate more with <ROTI>, the correlation coefficient was 0.63, indicating that there were periods of time where <AE> data was disturbed and <ROTI> was not. Furthermore, from Figure 3, it is apparent that the relationship between <AE> and ROTI is non-linear. Thus (AE) data is not ideal for forming <ROTI> predictions, since it is measured along many longitudes (twelve to be exact) while ROTI measures localized irregularities. 

Since (AE) did not correlate well with <ROTI>, new data were needed. Localized magnetometer data, provided by the SuperMag Organization, was the next main parameter used. This data was collected at the Alaskan station of College (station name: CMO, station coordinates: 64.87º geographic latitude and 212.14º geographic longitude). To obtain a better correlation with <ROTI>, a combination of the magnetic field’s horizontal and vertical components was calculated and used: First, Bh = , was computed where Bn and Be are the northern and eastern components of the magnetic field, respectively. Then, Bhz = Bh – Bz was determined to finally quantify <<Bhz>> = 5-minute average of Bhz for 24 hours, starting at 2.5 minutes Universal Time Coordinated (UTC). 

Additionally, several other measurements, such as solar radio flux, solar wind, IMF components, and geomagnetic field indices, were used for analysis, as shown below in Table 1.

Table 1: Parameters and Measurements

All the data used, including <<Bhz>> and <ROTI>, are measured in UTC. Tests were also conducted to determine if the correlation improves when data is compared at Local Solar Time using the following formula:

Using magnetometer data collected at the same longitude as the ROTI data gave an improved correlation. In fact, looking at the same day, March 01, 2015, the correlation between <<Bhz>> and <ROTI> is 0.82, an increase of 0.2 compared to using (AE) data. Referring to Figure 4, the periods of disturbed and undisturbed times match for both <<Bhz>> and <ROTI>. Additionally, Figure 5 shows a more linear relationship between <<Bhz>> and <ROTI> in comparison to that between (AE) and <ROTI>. 

Machine-Learning and Ridge Regressions

Three integral components of a ML algorithm are variance, bias (regularization), and overfitting, as seen in Figure 6 [4]. Variance is a measure of how accurately a model can change when presented with a different data set. High variance is associated with a model failing to fit the testing dataset. Bias is the inability of a method to capture the true relationship of data. Typically, bias is used to develop a model that can predict targets for data that follow a general pattern, rather than a specific one. Bias is associated with a model failing to fit the training set of data. Finally, overfitting occurs when a polynomial with a higher degree than what is needed is used to model data. Overfitting causes a model to perform well for training samples but to fail when it comes to generalizing. 

The algorithm used to produce <ROTI> predictions is a Ridge Regression. A Ridge Regression is a Linear Regression (sum of squares) with a small bias [4]. This results in a significant drop in variance. However, starting with a slightly deficient fit results in an improved long-term prediction since the same pattern that the model was trained on is not memorized. A Ridge Regression is ideal for working with data that shows multicollinearity tendencies, such as <ROTI> data. In comparison, an Ordinary Least Squares (OLS) model treats all variables equally and is unbiased. The OLS model therefore becomes more complex as the number of variables increases. As displayed in Figure 7, the optimal model complexity is at the middle of the plot, where the error is at its lowest. Unlike the OLS model, the regularization parameter of a Ridge Regression can be tuned to change the model’s coefficients and produce more accurate predictions. 

The Ridge Regression has several variables and parameters: 

• α = Hyperparameter that determines the regularization strength. The larger the alpha value, the stronger the regularization. A large alpha increases the bias of the model significantly. An alpha value equal to 1 is equivalent to a Linear Regression.
• w = Vector of coefficients
• B = Regression coefficients to be estimated
• e = Error residuals
• X = Matrix of independent variables (x)
• Y = Dependent variable (prediction)

The Ridge Regression equation, in matrix form, used to calculate desired predictions is shown in (1) [6]. 

                                           Y = XB + e                                  (1)

The loss function of the Ridge Regression is the linear least squares function, as shown in (2), and the regularization is given by the L2 norm, which calculates the distance of a vector coordinate from the origin of the vector in space. It is equivalent to the OLS model plus α multiplied by the summation of the squared coefficient values [4]. 

(2)

Overall, when calculating predictions, the Ridge Regression aims to minimize the loss function by finding a balance between minimizing the residual sum of squares as well as the coefficients (3) [4].

                                        (3)

Coefficients of the Ridge Regression are calculated by first taking the cost function, performing some algebra, taking the partial derivative of the cost function with respect to w, and setting it equal to 0 to solve (4) [4].

        (4)

The Ridge Regression model used to predict <ROTI> was programmed using Python and three main libraries: Pandas, Numpy, and Scikit-Learn. It was fit and tested on data from a total of 22 disturbed and stormy days with 18 days used for training and 4 days used for testing. The overall goal was to test which space weather and geophysical measurements would be useful for a ML prediction system. The test was performed with historical GPS and magnetometer data to predict the location, time, and intensity of concerned irregularities, under various geophysical and space weather conditions. 

Two measures of performance were used to determine the accuracy of the Ridge Regression model: Root Mean Square Error (RMSE) and R squared (R²). RMSE is the average magnitude of the residuals, or error, and is represented by (5). It is optimal to obtain an RMSE as close to 0 as possible [7].

      (5)

represents the proportion of variance for a target variable, explained by the independent variables, and is expressed as a percentage. It can be represented by (6), where is the sum of residual squares, , and is the sum of squares of the measurement deviation from the mean, . For our standards, an value between 30% and 50% is considered weak, an value between 50% and 70% is considered moderate, and an value greater than 70% is considered strong [7].           

                                                                (6)

Experimental Results 

As shown in Tables 2-4, many parameters and combinations of parameters were tested using a Ridge Regression model. The outputs of the Regression process include an ideal α value, the w vector, the Ridge Regression coefficients, the , RMSE value, and the correlation coefficient between the predicted and actual <ROTI> for the specific group of inputs.

Table 2: Ridge Regression Results (Part 1/3)

Table 3: Ridge Regression Results (Part 2/3)

Table 4: Ridge Regression Results (Part 3/3)

The combination of parameters with the highest correlation coefficient and value was <<Bhz>>, SYM/H, Kp, ap, and F10.7 (all in UT time) (Trial 32). The most impactful parameters were <<Bhz>> and Kp. The increase in from Trial 1, where <<Bhz>> was a single parameter, to Trial 12, where both <<Bhz>> and Kp were parameters, was ~9.7%. It is also evident in Figure 8 that the scattering of the predicted vs. actual ROTI plot became tighter than that of the other plots. It was observed that the parameters converted to LST had lower correlations with <ROTI> than those in UTC. For example, the value in Trial 7, where <<Bhz>> was in LST, decreased by ~8.6% in comparison to Trial 1 where <<Bhz>> was in UT. Additionally, <AE> was used as a parameter in Trials 6, 13, 16, 22, 23, 25, 30, 31, and 33, and it was clear that the regression with (AE) performed somewhat worse for predictions of irregularities. This may be because (AE) is a longitudinally averaged index while the irregularities are localized.   

Conclusions

Since the largest value obtained from the Ridge Regression predictions was 53.5%, the accuracy of the predictions is not ideal. There are several reasons why this may have occurred. To begin with, the data samples were relatively small, since only 18 days’ worth of data samples were used to train the model. However, this was intentional; the main objective was to experiment which parameters and combinations thereof correlated best with <ROTI>. A different model may have been more appropriate, but the Ridge Regression was selected since it allowed an easy transition between various inputs. It is important to note that the parameters selected were all space weather and geomagnetic measurements, while the irregularities may be driven more directly by other nonlinear processes, contrasting the Ridge Regression’s linear approach. 

Understanding how and why these different physical phenomena are correlated is crucial in making accurate predictions. As the development of the ML algorithm continues, its output will aid and contribute to the representation and understanding of the variations of the ionospheric state. The goal for the algorithm, once it is implemented, is to predict the location, time, and intensities of ionospheric irregularities under various geophysical and space weather conditions, reducing the impact of ionospheric scintillation on receptions of GNSS signals and data and improving navigation systems. 

Acknowledgements

The author would like to thank Dr. Xiaoqing Pi for providing research guidance, technical assistance, and for processing and providing the ROTI data; Caltech Northern California Associates for the fellowship; Peyman Tavallali, Jet Propulsion Laboratory, California Institute of Technology, for suggesting the use of the Ridge Regression Algorithm; IGS for the management of the GNSS Data Archives; and SuperMag Organization for the distribution of magnetometer data.

 

References 

[1] (2020, August 03). Retrieved from https://www.faa.gov/about/office_org/headquarters_offices/ato/service_units/techops/navservices/gnss/waas/ 

[2] Pi, X., A. J. Mannucci, U. J. Lindqwister, and C. M. Ho, Monitoring of global ionospheric irregularities using the worldwide GPS network, Geophys. Res. Lett., 24, 2283, 1997.

[3] Pi, X., Mannucci, A.J., Valant-Spaight, B., Bar-Sever, Y., Romans, L. J., Skone, S., Sparks, L., Hall, G. Martin, “Observations of Global and Regional Ionospheric Irregularities and Scintillation Using GNSS Tracking Networks,” Proceedings of the ION 2013 Pacific PNT Meeting, Honolulu, Hawaii, April 2013, pp. 752-761.

[4] Maklin, C. (2019, August 19). Ridge Regression Python Example. Retrieved from https://towardsdatascience.com/ridge-regression-python-example-f015345d936b 

[5] Qshick. (2019, January 03). Ridge Regression for Better Usage. Retrieved from https://towardsdatascience.com/ridge-regression-for-better-usage-2f19b3a202db

[6] Chapter 335 Ridge Regression – NCSS .Chapter 335 Ridge Regression Introduction Ridge Regression is – [PDF Document]. (n.d.). Retrieved from https://fdocuments.in/document/chapter-335-ridge-regression-ncss-chapter-335-ridge-regression-introduction.html 

[7] Singh, D. (2019, May 17). Deepika Singh. Retrieved from https://www.pluralsight.com/guides/linear-lasso-ridge-regression-scikit-learn

[8] Magnetometer Data. (n.d.). Retrieved from http://supermag.jhuapl.edu/mag/?fidelity=low&start=2001-01-01T00:00:00.000Z&interval=1:00:00