PyTransport: A Python package for the calculation of inflationary correlation functions

PyTransport constitutes a straightforward code written in C++ together with Python scripts which automatically edit, compile and run the C++ code as a Python module. It has been written for Unix-like systems (OS X and Linux). Primarily the module employs the transport approach to inflationary cosmology to calculate the tree-level power-spectrum and bispectrum of user specified models of multi-field inflation, accounting for all sub and super-horizon effects. The transport method we utilise means only coupled differential equations need to be solved, and the implementation presented here combines the speed of C++ with the functionality and convenience of Python. This document details the code and illustrates how to use it with a worked example. It has been updated to be a companion to the second version of the code, PyTransport 2.0, which includes functionality to deal with models of inflation with a curved field space metric.


I. INTRODUCTION
PyTransport is distributed under a GNU GPL licence.The most recent version can be obtained by visiting transportmethod.com.If you use PyTransport you are kindly asked to cite Ref. [1] as well as the archive version of this user guide in any resulting works.
The main purpose of this document is to teach those interested how to use, and if so desired adapt, the PyTransport package.It has now been updated to be a companion to the second version of the code, PyTransport 2.0, which includes functionality to deal with models of inflation with a curved field space metric.We have had to make some minimal syntax changes in this second version in order to support new functionality, as discussed below.Users of the original package will unfortunately not be able to switch to the new one without amending their scripts.The original user guide can still be found as the arXiv version 1 of this document.The philosophy behind the implementation is simplicity and ease of use.Python was selected as the language though which to interact with the code because it enables rapid scripting and provides a flexible and powerful platform.In particular, it has many readily available tools and packages for analysis and visualisation, and for tasks such as parallelisation (using for example Mpi4Py).
As an interpreted language, however, Python can be slow for some tasks.This is circumvented here by using C ++ code, which is compiled into a Python module, to perform numerically intensive tasks with the result that the speed of the package is nearly indistinguishable from pure C ++ .The C ++ code itself is kept as simple and clean as possible and can therefore easily be edited if required.PyTransport has been developed on OS X using Python 2.7.We have also performed limited testing on Linux systems, and attempted to ensure compatibility with versions of Python 3. It can also be adapted to Windows systems, but this functionality has not yet been incorporated into the released package 1The code is intended to be a reusable resource for inflationary cosmology.It enables users to quickly create a complied Python module(s) for any given model(s) of multi-field inflation.The primary function of the complied module is to calculate the power-spectrum and bi-spectrum of inflationary perturbations produced by multi-field inflation.To this end, the module contains a number functions that can be called from Python and that perform tasks such as calculating the background evolution of the cosmology, as well as the evolution of the two and three point functions.We also provide a number of further functions written in Python that perform common tasks such as calculating the power spectrum or bispectrum over a range of scales by utilising the compiled module.The true power of the approach, however, is that users can rapidly write their own scripts, or adapt ours, to suit their own needs.
The transport approach to inflationary perturbation theory that the code employs can be seen as the differential version of the integral expressions of the In-In formalism.It is helpful numerically because it provides a set of ordinary differential equations for the correlation functions of inflationary perturbations.The code solves these equations from deep inside the horizon until some desired time after horizon crossing using a standard variable step size ordinary differential equation (ODE) routine with error control.Such off the shelf routines are extremely well tested, and provide an easy way to change the required accuracy.This is helpful in order to check convergence of the numerical solutions, or to respond to needs of models with very fine features.Details of the transport method itself that the code is based on can be found in the recent papers [1] and [2], the second of which updates the method to allow for the analysis of models with a curved field space metric.We highly recommend reading this guide in combination with those papers.
In this guide, we first we give some brief background and motivation for the code, much more can be found in Refs.[1,2], before giving an overview of its structure and how it can be set up.In the appendices we give some more detail about the structure of the underlying C ++ code, give full details of all the functions the complied module provides, and all the functions provided by Python scripts which accomplish common tasks.The best way to learn how to use the package, however, is by example.We present an extended example below spread between the "Getting going" and "Examples" sections, complete with screen shots of the code in use.Other examples that come with the distribution are discussed in the "Examples" section.Throughout, familiarity with Python and to some extent C ++ is assumed, though in reality users can just probably get a long way by looking at the examples and modifying to their needs.
Finally, we would also like to refer readers to the complementary package developed in tandem with the work in Ref. [1] and with PyTransport: CppTransport [3].This is a platform for inflationary cosmology developed fully in C ++ and recently also updated to deal with curved field space metrics [4].In comparison with PyTransport it has more external dependancies (in the sense that the dependancies of PyTransport are mainly Python modules), but provides more sophisticated parallelisation and data management capabilities.In limited testing it is also found to be marginally faster.For users with modest aims in terms of CPU hours and data generation, however, it is likely to have a higher overhead in getting started, but may well be beneficial for intensive users.PyTransport is intended to be more lightweight with users encouraged to utilise the power of Python in combination with PyTransport to achieve their specific aims and data management needs.

II. BACKGROUND
Calculations of the correlation functions of perturbations produced by inflation are now extremely mature.In the single field context, the In-In formalism is routinely used to calculate the equal time correlation functions of the curvature perturbation, ζ, as wavelengths cross the cosmological horizon [5][6][7][8] where they become constant [9,10].For many models this calculation can be accurately performed analytically, while for others, such as models with features, a numerical implementation is required [11][12][13][14].If additional fields are added, the problem becomes even more complex.ζ is no longer necessarily conserved after horizon crossing, and the evolution of all isocurvature modes needs to be accounted for -all the way from the initial vacuum state until such time as the system becomes adiabatic, or until the time at which we wish to know the statistics (see for example Ref. [15] for a discussion of adiabaticity).While analytic progress can be made in some circumstances using the In-In formalism and/or so called "super-horizon" techniques such as the δN formalism, in general for multiple field models numerical techniques become even more important.
The code documented here accounts for all tree-level effects present in multi-field inflation.This includes the super-horizon evolution of ζ, which can occur in models with multiple light fields, as well as the effect of features in the multidimensional potential, and the effect of quasi-light or heavy fields orthogonal to the inflaton (which are important if the inflationary trajectory is not straight).As discussed above the code utilises the transport approach to inflationary correlation functions [16][17][18][19][20][21].This approach can be viewed as a differential version of the integral expressions of the In-In formalism, and evolves correlations of inflationary perturbations from their vacuum state on sub-horizon scales until we wish to evaluate the statistics.We note for clarity that in its original form the "moment transport" method was restricted to super-horizon scales, but it was later shown how it could be extended to include sub-horizon scales in Ref. [21].A recent paper studies this extension further and develops it into a working algorithm [1] with many additional details provided.The present document details the code PyTransport which is discussed in that paper.
At the background level an inflationary cosmology is completely determined by the evolution of the fields, φ i , and their associated velocities (the rate at which the fields change with cosmic time), φi , as a function of the number of e-folds (the logarithm of the scale factor, N = ln(a)) which occurs.At the perturbed level, the key objects are the correlations of the perturbations in these fields, and correlations of other perturbative quantities.Here we have used the label i to run over the number of fields present.The code numerically solves the equations of motion for the background fields, and the equations of motion for the evolution of correlations of the field and field velocity perturbations defined on flat hyper-surfaces.Ultimately the quantities probably of most interest for observations are the statistics of the curvature perturbation ζ -in particular the power spectrum and bispectrum -which the code calculates from the field space correlations.Defining the array X = {Q I , P J } of the covariant field space perturbation and its momentum 2 , where the components are labelled X a and a now runs over the total number of fields and field velocities, we recall the following definitions for later clarity: where P ζ and B ζ are power spectrum and bispectrum of ζ respectively, and Σ and B are the equivalent functions for the correlations and cross correlations in field space, and f NL is the reduced bispectrum.Σ and B together with the background values of the fields and field velocities are the objects directly evolved by the code using the equations detailed in section 5 of Ref. [1] with initial conditions detailed in section 6 of Ref. [1] 3 .As discussed in that paper, B and Σ can then readily be converted to give P ζ and B ζ through the use of the "N " tensors with components N a and N ab , described in section 7 of Ref [1] and updated to the case of a curved field space metric in Ref. [2] (see also [23]).The equations of motions for the correlations are given in Eqs.(5.5) and (5.16) of Ref. [1], the initial conditions in Eqs.(6.2) and (6.7), and the conversion to ζ in Eq. (7.4).
It is worth briefly commenting on how our code compares with existing ones.A number of publicly available codes exist to calculate the power spectrum from canonical multi-field inflation.For example Pyflation [24] and MultiModeCode [25] employ a method originally used by Salopek and Bond [26] in which the mode functions of the QFT of inflationary perturbations are evolved.These codes are in written Python4 and Fortran respectively.Moreover, a Mathematica code which implements the transport method for curved field space metric models, but is restricted to the power-spectrum, is also available [27].At the level of the three-point function a number of authors have undertaken numerical work directly utilising the In-In formalism, for example in Refs.[11][12][13][14]28].The only publicly released code we are aware of, however, is BINGO [13] which is restricted to single field models.No general multi-field codes have been undertaken until now.
A further advantage of PyTransport (and CppTransport) over previous codes is that it leaves little for the user to calculate analytically.It needs the user only to provide the inflationary potential.Then all the derivatives are automatically calculated using symbolic Python (SymPy) and written automatically into the C ++ code which is then compiled.Compared to Fortran or pure C ++ implementations PyTransport has the advantage of easy access to the extensive and easy to use modules available to Python, and compared to a pure Python or Mathematica implementation we have the advantage of speed.

III. CODE OVERVIEW
The code structure should become familiar though the extended example we provide, but here we give a brief summary.
The code is distributed in a folder called PyTransportDist/5 , which also contains a copy of this document (possibly updated compared with the arXiv version) in the PyTransportDist/docs/ folder.The base code for PyTransport is written in C ++ and has a simple object orientated structure.This code can be found in the PyTrans-portDist/PyTransport/CppTrans folder and we provide a few more details about its structure and functionality in appendix 1.The C ++ code is deliberately as simple as possible to ensure transparency and adaptability.The idea of the PyTransport package as a whole is that after a potential and a field space metric (if the metric is non-Euclidean) are provided by the user the C ++ code is automatically edited and complied into a Python module by supporting Python functions (called from the PyTransportDist/PyTransport/PyTransSetup.py file which is described in full in appendix 2), meaning a lot of work is done for the user.The end result is a Python module consisting of a set of Python functions for a specific inflationary model, called the PyTrans*** module.The functions of this module provide key routines for inflationary cosmology (including calculating the evolution of the two and three point correlations).The asterisks, ***, indicate we can label the module with a tag telling us what model it corresponds to, and we can therefore install multiple modules if we want to work with many models simultaneously.The key functions available to these modules are defined in the file PyTransportDist/PyTransport/PyTrans/PyTrans.cpp (which is a C ++ file defining the Python module), these functions are detailed in appendix 3. The scripts that edit the C ++ code and compile the module are discussed further below in the setup section, and by default they place the compiled module in the local folder PyTransportDist/PyTransport/PyTrans/lib/python/ to avoid access issues if, for example, you do not have root privileges.Other useful Python functions that perform common tasks, such as producing a power spectrum by looping over calls to the compiled module, can be found in PyTransportDist/PyTransport/PyTransScripts.py, and we describe them below, and in full in detail in appendix 4. Python treats functions written in Python inside a file, such as PyTransScripts.pyand PyTransSetup.py, in the same way as a compiled module.So there are effectively three modules within PyTransport, one to setup a compiled module for the potential we want to study (Py-TransSetup), the compiled module itself (PyTrans***) (or multiple complied modules labeled with different tags) and a module with various functions automating common tasks that use the functions of the compiled module (Py-TransScripts).Also in the PyTransportDist/ folder is an example folder PyTransportDist/Examples containing the examples discussed below.There are no dependancies external to the folders provided except for a working Python installation (with appropriate packages downloaded), and a C ++ compiler -this is deliberate to make the code as easy as possible to use.An MPI installation such as openMPI is also needed if the module is required to be used across multiple cores.
We note that all the C ++ code is written by the transport team except for an included Runge-Kutta-Fehlberg (rkf45) integrator routine written by John Burkardt and distributed under a GNU LGPL license detailed here 6 .We choose this lightweight integrator over other options, such as using integrators included with the BOOST library, in order that it could easily be included with the distribution with no external dependancies being introduced.In our (limited) testing it functions well for all the models we have looked at.There are no dependancies external to the folders provided except for a working Python installation (with appropriate packages downloaded), and a C ++ compiler -this is deliberate to make using the code as easy as possible to use.An MPI installation such as openMPI is also needed if the module is required to be used across multiple cores.

IV. SETUP A. Prerequisites
So what is needed?The idea is as little as possible beyond Python: Python: A working Python installation is needed, in development we used Python 2.7 which we recommend, but have subsequently attempted to ensure compatibility with versions of Python 3.For convenience we recommend downloading a complete Python distribution, for example Enthought Canopy or Continuum Anaconda7 , which come with all the core packages used by the code as well as interactive development environments.Python packages currently used by PyTransport or by provided examples include Numpy, Matplotlib, SciPy, Gravipy (needed only for models with an non-trivial field space metric), SymPy, Distutils, Math and Sys, as standard, and Mpi4Py and Mayavi are used for MPI and 3D bispsectra plots respectively.Of these only Mpi4Py and Mayavi may need to be downloaded separately from the distributions mentioned.One way to install a package such as Mpi4Py is to type "pip install Mpi4Py" in the terminal.If using Canopy, Anaconda or similar, they come with their own package managers which are an even easier way to install packages.There are many easily found resources on the internet to help with such installations if there is a system specific snag.Note that you should not attempt to install Mpi4Py without installing MPI first (which we deal with next).We also note that although apparently possible we have not easily been able to install Mayavi with Python 3.5, and recommend searching for online resources to help with this.

MPI:
As computing the bispectrum can be computationally expensive, distributed computing can be helpful (even if only across the multiple cores of modern PCs).In some of the scripts in the PyTransScripts module, we use the Mpi4Py module to implement this.Mpi4Py needs a working MPI installation such as openMPI installed on your computer.Note that Mpi4py or openMPI are not needed for PyTransport in general, and if you do not have these installed you simply cannot run the scripts that use MPI, but can run the code in serial instead.A nice guide to installing openMPI is at this link8 .
C ++ complier: Python needs to be able to find a C ++ complier in order to compile the PyTrans module(s).This is bundled with most Linux distributions.If not present on a Mac system, downloading Xcode from the app store is the easiest way to install one (or Xcode command line tools can be downloaded separately, as explained here9 ).

B. Getting going
Once you have Python running and a C ++ compiler, to get started take the PyTransport/ folder from the PyTrans-portDist/10 folder and place it anywhere convenient in your computer's file system.It is essential that you don't change the structure of the sub-directories within PyTransport/ , but you can place this folder wherever you want.That's more or less all you have to do.You can do this by copying the entire folder PyTransportDist/ (which also contains examples and this guide) to a convenient location, but equally well you could run examples from anywhere else on your computer.In each example, we will see one needs to add the path of the PyTransport/ folder to the paths which Python includes when looking for code, so that Python can find the setup file PyTransSetup.py(or this could be done permanently).Now you can get started, no other installation is required which is not handled for you by provided Python scripts.
Lets say you want to analyse a canonical inflationary model (we will give an example of a model with a curved field space metric in Section V B) defined by the potential The first step is to create (by compiling the C ++ code) a Python module for this potential.This is achieved by writing the potential in a Python file in SymPy (symbolic Python notation) and calling the appropriate functions from the PyTransSetup module.First define two SymPy arrays one for the fields and the other for parameters which define the model (the parameters you might wish to change the value of).These have length nF and nP respectively.Then define the symbolic potential V using these arrays.The potential must be written into the C ++ code by calling the function potential(V,nF,nP) which is in the file PyTransSetup.py.If we wanted to work with a non-Euclidean field space metric it would also be specified at this stage, but if one is not specified explicitly the code assumes the model is canonical.If working with a version of Python 2, the next step is to call the function compileName("Quad") (where "Quad" can be replaced by any name the user likes (the *** from above)).If using Python 3, users must use the function compileName3("Quad") to achieve the same thing.Multiple modules for different potentials with different names can therefore be created and used simultaneously.In the previous version of the code the tolerances used by the numerical integrator had to be specified at this setup stage, but these can now be set at the point we want to calculate a given evolution (and hence more easily adjusted).Below is a screen shot of the procedure just described with copious comments which should make the procedure clear: This example is contained in the Examples/DoubleQuad/ folder which accompanies the code, with this script in the file DQuadSetup.py file.In this script we have used the two functions from the PyTransSetup module mentioned above.Appendix 2 contains a summary of all the functions available in the setup module.The complied Python module can now be used.To do so we need to point Python to the path of the new module (and the scripts module if we wish to call the provided Python scripts).This can be done automatically by calling the function in the setup module pathSet().We recommend using the PyTransQuad module in a separate file from the one used to set it up, and if working in an integrated development environment to restart the Python kernel (this is to ensure the most recent version is always imported).Below is a screen shot of the start of a file in which we use the module we set up in the previous paragraph.In the screen shot we first calculate the value of the potential and the first derivative of the potential for a particular choice of field values and parameters.Then we use these to set up an array containing field values and the associated fields velocity (using the slow roll equation): This screen shot is of the start of the file Examples/DoubleQuad/SimpleExample.py.Of course we will usually want to use the module for more sophisticated tasks.Appendix 3 contains a summary of all the functions available within the PyTrans*** module.We will see the use of a number of the more sophisticated functions in the Examples section.

V. EXAMPLES A. Double quadratic
First lets continue with the double quadratic example using more of the functions available from the compiled module.In the screen shot below we use the background evolution function to calculate a fidicual background trajectory in field space using the array we set up in the last part of the example as initial conditions.Here and for all the functions and output of the PyTransport package e-folds (N ) are used as the time variable.The function used to calculate the background evolution is the PyTrans.backEvolvefunction.Details of the format it outputs the background evolution in can be found in Appendix 3. Essentially it provides information about the fields and their rate of change (in cosmic time) at every e-fold value given by the array t.In addition to this array and the initial conditions, one must provide the parameter value used to define the potential, and the absolute and relative tolerances to be used by the integrator (this argument has been added in the second version of the code).The final argument of this function indicates whether the evolution should terminate at the end of inflation (and has also been added in the second version of this code).If set to true data is only returned up to the first value of t that is past the end of inflation if that point is reached (note it does not find the exact end of inflation, and the last value returned will only be close to the end of inflation if the array t finely samples the evolution).If set to false, the code will attempt to generate output for all the times contained in t.If the last entry is long past the end of inflation the code could take a long time to run (or crash) as typically the number of oscillations in field space grows exponentially with N .The output is plotted shown in Fig. 1.
Then we run the two point evolution function to calculate the evolution of Σ and the power spectrum of ζ for a k mode which crossed the horizon 15 e-folds into this fiducial run using the PyTrans.sigEvolvefunction.We plot the correlations and cross correlations of the fields in Fig. 2. We repeat for a neighbouring k to give us a crude estimate of the spectral index, n s .Finally, we use the PyTrans.alphaEvolvefunction to calculate the evolution of the field space three-point function and the bispectrum of ζ for a set of three ks.We plot the three-point correlations and cross correlations of the fields in Fig. 2 and also the evolution of the reduced bispectrum, f NL .In the plots it can clearly be seen that the heavier field drops out of the dynamics at around 40 e-folds.At this point the system becomes adiabatic and ζ and its statistics become constant.A screen shot of the code which does all this from the SimpleExample.pyfile is below: There are few things to note from this script.
First we note that it is important that we fix our final time sensibly.If we were to fix it after the end of inflation, when the lighter field oscillates indefinitely about its minimum, then the code would become very slow, especially when it evaluates the correlations (since the the field correlations also oscillate a lot in response).This can be achieved by using a time array which runs up to the last output time of the background evolution (if, as in this example, we asked the background evolution to terminate at the end of inflation).This final time may not be what we require (if for example the moments of ζ becomes constant before the end of inflation), and in this case we could use an earlier time.
Next we note the use of a function within PyTransScripts.pywhich finds the k value which corresponds to an exit time of 15 e-folds after the start of the fiducial run, the PyTransScripts.kexitNfunction.This function uses the background trajectory and Python spline routines to find this k.We also note that it is essential that we run the two and three point correlation evolution from a time the k mode of interest is deep inside the horizon.In the script, we calculate this time as well as the field and field velocity values at this time (which are then fed into the two and three point evolution routines) using another function from the scripts module, the PyTransScripts.ICsBE function.
To use this function we need to specify the number of e-folds we require before horizon crossing, and here we specify 6.0.The function returns a time roughly 6 e-folds before the horizon crossing time and a numpy array containing the fields' values and velocities at that time.It is important to point out that this function is only approximate in the sense that it simply finds the first value of N in the background array (back) which is before the specified number of e-folds before horizon crossing, and returns this value and the corresponding field and field velocities at this time.It therefore requires the background fiducial evolution which is fed into it to be finely sampled (10 points per e-fold is a rough guide) to be accurate.
Finally, we note that within PyTransScripts.pywe also include the related functions kexitPhi which finds the k value which crosses the horizon at a particular field value, and and ICsBM which finds initial conditions a fixed time before the "massless condition" which is discussed further in section V C.These functions and all others in this module are detailed in appendix 4. Calculating the value of n s in the manner presented here is clearly a bad way of doing things, since it involves using only two points in the power spectrum to calculate a derivative, and the step between them is arbitrarily chosen.An alternative is to calculate the power spectrum around a given exit time using a number of points, and to fit a spline to it and differentiate.We now generate the power spectrum for this model with the following script which fits a spline to the entire spectrum, differentiates and produces n s at every value of k over roughly 30 e-folds, the results are plotted in Fig. 3 (the start of the file is identical to that above): In this script we used the PyTranScripts.pSpectrafunction to generate the power spectrum over a range of ks.This function essentially just loops over calls to the compiled function which evolves the two-point function.
Next we wish to calculate the bispectrum.Here we first we calculate the bispectrum in the equilateral triangle configuration as a function of the k values we calculated the power spectrum for above.Then we generate (and plot using a separate plots file) a slice through the bispectrum for a given k t as a function of the α, β variables defined such that k 1 = k t /2 − βk 2 /2, k 2 = k t /4 * (1 + α + β) and k3 = k t /4 * (1 − α + β).We use two separate scripts for each of these tasks (and the plots file) which are pasted below, both use MPI to speed things up.They should be called using the command "/usr/local/bin/mpiexec', '-n', '10', 'python', 'MpiEqBi.py'" (for the equilateral run, where the first part should be replaced with your location of mpiexec if different, and 10 replaced by the number of processes one desires to call): The results are in Figs. 4 and 5 respectively.Note that in all these scripts we use more functions available from the PyTransScripts module whose function should be self evident and which are described in full in appendix 4. If one didn't wish to use MPI the only change needed would be to call the function alpBetSpec rather than alpBetSpecMpi from the scripts module, and to remove MPI related lines and reference to rank etc.When using MPI we recommend calling more processes that the system has cores.This is because we have not implemented sophisticated load sharing, and since some ranges will be faster to evaluate than others, if the number of processes is larger than cores, the cores that have capacity first will end up running more processes, sharing the load in a simple way.In this example we also save the data at the end for future use.This can be done simply for any numpy array in various ways and then read back into Python easily.For reasons of simplicity and flexibility we leave data management up to the user, but note that Python is a powerful tool for this purpose.
If we want to generate the full bispectrum we would simply loop over the alpBetSpec or alpBetSpecMpi function for many k t .

B. Double quadratic with a field-space metric
To demonstrate the evaluation of a non-canonical inflation model with a curved field-space metric we extend the previous example to include a non-trivial metric.We retain the potential from Eq. (4.1), but now consider the model to have the field metric (5.1) Previously we needed to create a Python module for a model with the double quadratic potential (that was canonical), now a Python module must be created for a model with this potential and the field metric.This is achieved by writing out the setup script we had before but including a G matrix which encodes Eq. (5.1) in SymPy notation.This is then included as a final argument in the potential(V,nF,nP,Simple,G) function which sets up the files needed to compile this model into a python module (this final argument is optional, if not included then model must be canonical).The fourth argument is also optional and we didn't discuss it when dealing with the previous canonical example.It tells the routines which use sympy whether or not to attempt to simplify the expressions using SymPy's simplification routines.By default it is set to False for reasons discussed in Section VI A (this ability to switch simplification off has also been added to the new version of the code).Finally, the function compileName('DQuadNC','True'), or compileName3('DQuadNC','True') if the user is using Python 3, is used to create a python module PyTransDQuadNC, which can be used to study this model.The final argument in the compileName function is optional, and care must be taken to include it and set it to True if the model has a non-trivial field space metric.If set to true the then the code that deals with the compilation is pointed to different C ++ files which contain the more complex equations needed deal with models of inflation with a non-trivial field space.Below is a screenshot of setup file for this model.This example can be found in the Examples/DoubleQuadraticNC/ folder with some simple scripts similar to those above for the canonical case.The complied model can be used in exactly the same way as for canonical example.

C. Heavy field examples
In the previous examples both fields which played a role in the dynamics were light (at least at the start of the evolution).Interesting dynamics can also occur when the field orthogonal to the direction of travel in field space is heavy, if the field trajectory curves.In this kind of example it is imperative that the initial conditions for the evolution of the two and three point functions are set when the k 2 /a 2 term in the equation of motion for the scalar field perturbations dominates over the mass squared of the heavy field.This is a requirement for our initial conditions to be accurate as discussed in Ref [1].There is a script in PyTransScripts.pywhich will achieve this, the ICsBM function.This finds initial conditions a user specified number of e-folds before the massless condition where k 2 /a 2 = M 2 (where M is the largest eigenvalue of the mass matrix).There is also the function ICs which evaluates initial conditions using ICsBE (the before horizon exit function) and ICsBM, and takes the earliest one.The power spectrum and bipsectrum routines use this latter function.
One example is the potential: where which is from Ref. [29] and is in the Examples/LH/ folder with some simple scripts to those discussed above for the double quadratic potential for users to play with.This example is also discussed at length in Ref. [1].

D. Further examples
Also in the Examples/ folder is another light field example with more interesting dynamics than the double quadratic example which we refer to as the axion quartic model.This example is again accompanied with scripts and plots for users to explore and was also discussed in Ref. [1].It is in the QuadAx folder, and has the potential Finally in the Examples/ folder and discussed in Ref.
[1] is a single field example with a step in the potential: which is in the folder SingleField/ and was discussed in Refs.[11,12], as well as in Ref. [1] .

VI. THINGS THAT CAN GO WRONG
While using the PyTransport package some issues have presented themselves which it might be useful for new users to know about.Many more have been ironed out, but it is neither practical nor desirable to make the code fully immune to misuse.Below we detail a few common problems we have faced.

A. Potential computation fails and problems with simplification of expressions
This is the most severe bug/problem we have found, but also the least common.The function PyTransScripts.potential()takes a potential V written in SymPy format (and a field space metric, if the model is not canonical) and calculates and then simplifies related functions such as the derivatives.To do so it uses the function sympy.simplify.As discussed at this reference11 , however, this is not always the best way to simplify an expression and can take some time to complete the simplification for complicated functions.For example the simplification process takes a relatively long time (tens of seconds) for the heavy field model above.
We found a more serious problem occurred when looking at the example: which represents a semi-circular valley in field space.For this example, there appears to be a simpy.simplifybug which made it crash with the potential written in the form given above, when the powers where written as doublesi.e. as 2.0 -rather than simply as 2. The problem appears when taking cross derivatives.Uncommon errors such as this are something users might need to watch out for.
We have also noticed that it is helpful to write powers as integers (such as 2, rather than 2.0) in general, as this can result in much shorter expressions after simplification.If analysing a complex model, or if an error is encountered, the user can inspect the potential.hfile or the fieldmetric.hfile to see what expressions Python has generated.The location of these files is discussed in appendix 1.
Because of these issues simplification is by default switched off, but can be useful to speed up the code.

B. Make sure latest module version is imported
A more common problem is that if we wish to update a complied PyTrans*** module, for example after altering the tolerances, then after recompiling the module we need to ensure the new module is imported.The only reliable way to do this seems to be to restart or open a new Python session and then use the import command.If working in the Canopy editor, for example, this can be achieved by selecting "restart kernel" from the "run" menu.

C. Selecting the absolute and relative tolerances
The evolution of the three point function can be a numerically intensive task, requiring high numerical accuracy.The question arises how low (the lower the higher the accuracy) do we need to set numerical tolerances.This question can't be answered absolutely, and must be dealt with on a model by model basis.Models with finer features in the potential, or in which the excitation of the two and three point function occurs on sub-horizon scales will require lower tolerances (high accuracy).Models which produce a small signature may also need higher accuracy to resolve the true answer from noise than models which produce a large bispectrum.
Convergence is the key criterion in selecting tolerances.If one calculates the evolution of the three-point function of ζ for a example run with a given potential, reduce the tolerances and run again without the answer changing in any significant manner, then the tolerances are likely to be sufficient.As a rule of thumb 10 −8 for the absolute tolerance and 10 −8 for the relative tolerance is usually sufficient.However, some examples do require lower values, while for reasonably accurate results for simple models one can sometimes get away with higher values.As the values are lowered, the code takes longer to run and eventually will fail (described below).Therefore, there is significant benefit for picking a required accuracy which is sufficient for the task, but not one which is too stringent.Those with experience of solving ODEs numerically will be familiar with this problem, which is of course more general than the specific integration at hand.One can determine the accuracy needed for a given model by examining the output for various choices for a small number of runs and checking for convergence.

D. Integration stalling or failing
The consequences of picking a tolerance that is too demanding can be that the integrator will not finish.But this can also be a consequence of setting silly initial conditions or parameter choices.It is always a good idea once you begin with a new model of inflation to build up gradually.First integrate the background.Even here it is possible for the code to take a long time.For example, if the final time is after the end of inflation the code will try and track all the oscillations of the field.As these increase exponentially with e-fold it can be very time consuming.Once the background is giving sensible output, move onto integrating the correlations, initially just for a single k value of the power spectrum or single triangle of the bispectrum.If everything looks good then run over many values to calculate the power or bi-spectrum.If you feel you are waiting too long try just asking the code to evolve a short e-folding time (.1 say).Typically a single triangle of the bispectrum evaluated before the end of inflation will take from between a fraction of a second to half a minute to run with 4 − 6 efolds of sub-horizon evolution, depending on the required accuracy.Heavy field models however typically have many more sub-horizon efolds given that the massless condition is already met only deep inside the horizon.If the run time seems to be much longer than you expect, double check you are working with the correct potential for your initial conditions and parameter choices.

E. Integration failing
If the code can't reach the accuracy demanded by the user the rkf45 routine will stop running and issue an error message that the required accuracy couldn't be reached.

F. Not enough efolds before horizon exit
A requirement that needs to be met in order to get accurate power spectra and bispectrum is that all the k values involved in a given correlation must be sufficiently deep inside the horizon initially for the initial conditions to be accurate.Since the functions which find the initial conditions for the two and three point evolutions take a background trajectory as their input, this trajectory must have enough e-folds prior to the exit of the k values of interest so that the correct initial conditions can be found.
Moreover, we must choose how many e-folds these k values stay sub-horizon for.As described above we can either measure backwards from horizon crossing itself or from the massless condition (only suitable for models with heavy fields) or pick the earliest of the two conditions.Normally 4-5 e-fold from these points is about right.But as with the setting of tolerances, the only way to tell for sure is to demand convergence.Too little run in time will lead to spurious oscillations in the spectra.

G. Not enough efolds after horizon exit
For single field models or effective single field models, such as models additional heavy fields, the correlations of ζ become conserved after horizon crossing.Likewise if there are multiple additional light fields which then decay, leading to adiabatic evolution, ζ also becomes conserved.In the former case conservation only occurs a few e-folds after horizon crossing, and in the latter the process also takes a few e-folds to complete.We must ensure, therefore, that the statistics are evaluated at a time when ζ has become conserved if that is what is intended, and so our evaluation time must be sufficiently far after horizon crossing or field decay.

H. Computer crashes and data loss
The PyTrans*** compiled module returns data in memory to Python.Moreover, the functions provided to calculate power spectra and bispectra which loop over the complied module also return arrays of data in memory to the calling process.This is true even for the MPI functions.If anything goes wrong before these functions complete the data is lost.For example if the computer crashes, or if using distributed computing any one of process crashes.These modules can be easily edited by users to instead write data to disk periodically, particularly advisable for distributed computing.We have not included this functionality to avoid overly complicated functions, and in any case, we anticipate that different users will have different needs, but it is an issue to look out for.

VII. SUMMARY
We have presented the PyTransport package for the calculations of inflationary spectra.This package compliments and extends currently available tools, and is also complimentary to two related packages mTransport [27] and Cpp-Transport [1,3] all described at this website 12 .In its most recent version it has also now been updated to include functionality to deal with models that have a curved field space metric [2].
Through use of a detailed example we have shown how PyTransport can be used in practice.We have also summarised the structure of the code, with some more details provided in the appendices.
• PyTrans.backEvolve(Narray,fields-dotfields, params, tols, exit) takes a numpy array of times in efolds "Narray" at which we want to know the background value of the fields and field velocities.This must start with the initial time in e-folds (initial N) we wish to evolve the system from, and finish with the final value of N. It takes a numpy array of length 2 nF "fields-dotfields" which contains the background field and velocity values at the initial time, as well as a numpy array of length nP "params" containing parameter values.It also takes an array "tols" containing the relative and absolute tolerance and a boolean which if True will exit the numerical evolution when inflation ends ( = 1) or if False continues until the desired number of e-folds has elapsed.It returns a two dimensional numpy array.This array contains the fields, and field velocities at the times contained in Narray.The format is that the array has 1 + 2nF columns, with the zeroth column the times (Narray), the next columns are the field values and field velocity values at those times.
• PyTrans.sigEvolve(Narray,k, fields-dotfields, params, tols, full) takes a numpy array of times in efolds "Narray" at which we want to know the value of the two point function of inflationary perturbations.This must start with the initial N we wish to evolve the system from, and finish with the final N. It also takes a Fourier mode value "k", and the initial conditions of the background cosmology (field and field velocity values) at the initial time as a numpy array of length 2 nF, the parameters of the system as a numpy array of length 2 nP, an array "tols" containing the relative and absolute tolerance, and an integer "full" set to False or True (if another value it defaults to True).The initial time and the initial field and field velocity array are used to calculate the initial conditions for the evolution of the two-point function.The function returns a two dimensional numpy array the zeroth column of which contains the times (Narray).If full=False the next column contains the power spectrum of ζ at these times, and this is the final column (there are therefore only 2 columns in total).If full = True then there is in addition 2 nF + 2 nF * 2 nF columns.The first 2 nF of these columns contain the fields and field velocities at the time steps requested.Then the final 2 nF * 2 nF contain the elements of matrix Σ ab r (the power spectrum and cross correlations of the fields and field velocities).There are therefore [1 + 1 +2 nF +2 nF *2 nF] columns in total.The convention is that the element Σ ab r corresponds to the [ 1 + 2nF + a + 2nF x (b-1)]th column of the array (recall the columns start at the zeroth column -a and b run from 1 to 2 nF).
• PyTrans.alpEvolve(Narray,k1, k2, k3, fields-dotfields, params, tols, full) takes a numpy array of times in efolds (N) at which we want to know the value of the three point function of inflationary perturbations.This must start with the initial N we wish to evolve the system from, and finish with the final N. It also takes three Fourier mode values (k1, k2, k3), the initial conditions of the background cosmology (fields and field velocities) at the initial time as a numpy array of length 2 nF, the parameters of the system as a numpy array of length 2 nP, an array "tols" containing the relative and absolute tolerance, and an integer "full" set to False or True (if another value it defaults to 1).The function returns a two dimensional numpy array, the first column of which contains the times (Narray).If full=False the next four columns contains the power spectrum of zeta for each of the three k values input, and the value of the bispectrum of zeta for a triangle with sides of length of these k values.A total of [1+ 4] columns.If full=True, there are an additional 2*nF + 6 *2 nF* 2 nF + 2 nF * 2 nF * 2 nF columns.The first 2nF of these columns contain the fields, and field velocities at the time steps requested (the background cosmology).The next 2nF * 2nF of these contain the real parts of Σ ab (k 1 ) in the same numbering convention as above.Then the real part of Σ ab (k 2 ) and then the real parts of Σ ab (k 3 ), the following 3 * 2nF * 2nF columns are the imaginary parts of the Σ ab (k 1 ), Σ ab (k 2 ) and Σ ab (k 3 ) matrices.So for example if one wanted access to the Σ ab i (k 2 ) element that would be the  The functions which are part of the PyTransScripts modules are detailed below.The code which provides these scripts in the PyTransScripts.pyfile should be clear to those familiar with Python.The scripts are simply an indication of what is possible, and it is intended that users will modify them for their own purposes, or write their own, as well as using the ones provided.
• PyTransScripts.ICsBE(NBExit, k, back, params, PyT) takes in the number of e-folds before horizon exit of a scale k at which initial conditions for the evolution of correlations are to be set, the scale k itself, a numpy array "back" containing the background cosmology (returned from PyTrans.backEvolve), the parameters of the model as a numpy array, and the PyTrans module being used.It returns a double, and an numpy array.The former is the starting number of efolds which are at least NBExit before exit, as measured with respect to the background trajectory, "back".The array contains the initial value of the fields and field velocities at this starting time.This script simply runs through the elements of back to find the first one before the exit time minus NBExit.As such it requires back to have enough entries for the result to be useful (roughly 10 per e-fold is fine), as discussed in the main text.
• PyTransScripts.ICsBM(NBMassless,k, back, params, PyT) works like PyTransScripts.ICsBM, but instead of calculating initial conditions before exit time, it evaluates the time when k 2 /a 2 = M 2 where M is the largest eigenvalue of the mass matrix of the potential (we call this the massless condition), and returns conditions more than NBExit before that time.This is useful for example with a heavy field, for which we need to ensure the approximation we use to fix initial conditions is accurate (which requires k 2 /a 2 M 2 ).
• PyTransScripts.ICsBM(NB,k, back, params, PyT) takes the same arguments as the two previous functions and calls each in turn, then outputs the number of e-folds and the fields and field velocities at the earliest time of the two.This is so that the start time can be set to either NB before exit or NB before the massless condition, whichever is earlier.
• PyTransScripts.kexitN(Nexit,back, params, PyT) takes an exit time in e-folds, a background evolution which runs though this time (returned from PyTrans.backEvolve), a set of parameters and the PyTrans module being used.It returns as a double the k mode that exited at the time Nexit.This routine uses NumPy spline routines to find the value of k at the exit time.
• PyTransScripts.kexitPhi(PhiExit,n, back, params, PyT) takes an exit value of one of the fields, and a number indicating one of the fields (in range from 1 to nF), a background evolution which runs through this field value, a set of parameters and the PyTrans module being used.It returns the k mode that exited at the that field value.This routine uses NumPy spline routines to find the right k.
• PyTransScripts.pSpectra(kA,back, params, NB, tols, PyT): takes a numpy array specifying a range of k "kA", a background evolution (output from the backEvolve function), a set of parameters, a double indicating the number of e-folds (before massless or before exit whichever is earlier) of sub-horizon evolution, an array containing the relative and absolute tolerance, and the PyTrans module being used.It returns two numpy arrays.The first has the corresponding values of P ζ (to the input array of k) at the end of the evolution, and the second the times taken to perform the integration for each element.
• PyTransScripts.pSpecMpi(kA,back, params, NB, tols, PyT): does the same as the function above but spreads the calculation across as many process as are active using Mpi4Py.The script which contains this function should be called using the the command "mpiexec -n N python Script.py",where N is the number of processes to be opened.The length of kA must be divisible by N. We recommend calling at least twice as many processes as cores are available so that cores that run processes which finish first don't simply remain idle.The function returns the two numpy arrays to the process with rank 0. Empty arrays are returned to the other processes.
• PyTransScripts.eqSpectra(kA,back, params, NB, tols, PyT): takes in the same information as the function above, and returns three arrays.The first has the corresponding values of P ζ , the second the corresponding values of B ζ in the equilateral configuration (k1 = k2 = k3 = k) at the end of the evolution, and the final the times taken to perform the integration for each element.
• PyTransScripts.eqSpecMpi(kA,back, params, NB, tols, PyT): does the same as the function above but spreads the calculation across as many process as are active using Mpi4Py.The script which contains this function should be called using the the command "mpiexec -n N python Script.py",where N is the number of processes to be opened.The length of kA must be divisible by N. We recommend calling at least twice as many processes as cores are available so that cores that run processes which finish first don't simply remain idle.The function returns the three numpy arrays to the process with rank 0. Empty arrays are returned to the other processes.
• PyTransScripts.alpBetSpectra(kt,alpha,beta, back, params, NB, nsnaps, tols, PyT): takes in a value of kt, and two numpy arrays defining a range of values of alpha and beta, as well as a background evolution, set of parameters and a number of e-folds before horizon exit or the massless condition is met and an array containing the relative and absolute tolerance.It also takes an integer nsnaps, and the PyTrans module being used.nsnaps tells the code at how many different times at which to provide output.The function returns six arrays.The first output array is a three dimensional numpy array containing B ζ for k values corresponding to the values of alpha, beta and kt input, and at nsnaps different times between the start and finish time.The second array is also three dimensional and corresponds to P ζ (k 1 ) for these values and times, the next P ζ (k 2 ) and the next P ζ (k 3 ).If nsnaps is 0 or 1, the third dimension of these arrays is only 1 element long, and the values returned correspond to those at the end of the evolution.If it is greater than one the output is given at evenly spaced times until the end of the evolution.This allows us to see how a slice through bispectrum evolves with time if we wish.The fifth array is two dimensional and corresponds to the times taken to perform the integrations associated with every combination of alpha and beta.The last array is the times at which the nsnaps are taken.
• PyTransScripts.alpBetMpi(kt,alpha,beta, back, params, NB, nsnaps, tols, PyT): does the same as the function above but spreads the calculation across as many process as are active using Mpi4Py.The script which contains this function should be called using the the command "mpiexec -n N python Script.py",where N is the number of processes to be opened.The length of alpha must be divisible by N. We recommend calling at least twice as many processes as cores are available so that cores that run processes which finish first don't simply remain idle.The function returns the five numpy arrays to the process with rank 0. Empty arrays are returned to the other processes.

FIG. 3 :
FIG.3:The power spectrum and ns in the double quadratic model.