CMB constraints on cosmic strings and superstrings

We present the first complete Markov chain Monte Carlo analysis of cosmological models with evolving cosmic (super)string networks, using the unconnected segment model in the unequal-time correlator formalism. For ordinary cosmic string networks, we derive joint constraints on Lambda cold dark matter (CDM) and string network parameters, namely the string tension Gmu, the loop-chopping efficiency c_r and the string wiggliness \alpha. For cosmic superstrings, we obtain joint constraints on the fundamental string tension Gmu_F, the string coupling g_s, the self-interaction coefficient c_s, and the volume of compact extra dimensions w. This constitutes the most comprehensive CMB analysis of LambdaCDM cosmology + strings to date. For ordinary cosmic string networks our updated constraint on the string tension is, in relativistic units, Gmu<1.1x10^-7, while for cosmic superstrings our constraint on the fundamental string tension is Gmu_F<2.8x10^-8, both obtained using Planck2015 temperature and polarisation data.


I. INTRODUCTION
Cosmic strings are line-like concentrations of energy that can arise as topological defects in theories of the early Universe [1][2][3][4][5]. In particular, they form naturally in models of hybrid inflation [6][7][8][9][10][11][12] in which the inflationary phase ends with a second-order phasetransition [7,[13][14][15]. Although they were originally considered as an alternative candidate for providing the seeds for structure formation in the Universe [16][17][18][19], it is now understood that they cannot give rise to the observed acoustic peak structure in the power spectrum [20][21][22][23][24], but can play a subdominant role. There are a wide range of potential observational signatures of cosmic strings, for example line-like discontinuities in the cosmic microwave background (CMB) temperature anisotropy via the Kaiser-Stebbins effect [25,26]. Thus, strings provide a powerful tool for testing theories of the early Universe. Observations have strongly constrained the contribution of cosmic strings to the total CMB anisotropy [20,[27][28][29][30][31][32]. Current data place a 2σ upper bound on the string tension of Gµ < 1.3 × 10 −7 for Nambu-Goto strings [33] or Gµ < 2.7 × 10 −7 for Abelian-Higgs strings [34], which corresponds to ∼1% of the total temperature anisotropy at = 10. G is the gravitational constant, µ is the tension of the string and c = 1 in relativistic units. Although this may seem insignificant, there is still constraining power left in the data since strings generate specific signatures in the primordial B-mode polarisation spectrum [27,[35][36][37][38][39][40], which can now be analysed with the Planck2015 polarisation [41] and joint BICEP2 data [42].
Going beyond the simplest cosmic string models, complex networks of multiple types of interacting superstrings, each with a different tension, can also be considered. Notably, interacting networks of fundamental F-strings, one dimensional D-branes (Dstrings) and bound (FD) states between F-and D-strings, collectively referred to as cosmic superstrings, arise naturally in string theoretic inflation [7,43,44]. These networks are notably different to their simpler, single-type string counterparts since the different string types have intercommutation probabilities that are not necessarily unity [44][45][46][47][48][49][50]. The interactions among different string types are also much more complex, as colliding strings can zip together or unzip, producing heavier or lighter FD-string states carrying different charges. These features affect CMB signatures allowing us to obtain constraints on string theory parameters such as the string coupling g s and the fundamental string tension µ F [51,52].
In this paper we use the Planck2015 public data [41] to perform the first full Markov chain Monte Carlo (MCMC) analysis of ΛCDM models with cosmic string or superstring networks. For "ordinary" cosmic string networks we work in the unconnected segment model (USM) framework and utilise our analytic method [53] for fast computation of the string unequal-time correlator (UETC). This is used as a source to compute CMB anisotropies and hence obtain joint constraints on ΛCDM and the string network parameters, including the tension Gµ, the loop chopping efficiency c r and the wiggliness parameter α. In the case of cosmic superstring networks we extend our method to deal with multiple network components. The UETC approach is efficient, meaning we can compute the superstring spectrum in much less time than previous codes and obtain joint constraints on the fundamental string tension Gµ F , the string coupling constant g s , the self-interaction coefficient c s , and the parameter w of [52], quantifying the volume of compact extra dimensions.
In Sec. II we describe the UETC formalism applied to evolving Nambu-Goto string networks. In Sec. III we summarise our modelling of cosmic superstrings and the adaptation of our UETC method to these multi-string component networks. In Sec. IV we present the results of our MCMC analysis for cosmic string and superstring networks using Planck2015 CMB data. Our constraints on string network parameters and possible future directions are discussed in Sec. V.

II. UNEQUAL-TIME CORRELATOR
Unlike passive inflationary perturbations which are set as initial conditions, metric perturbations from cosmic string networks are actively sourced at all times. To compute the string spectra the components of the string network's energy momentum tensor must be used as sources in the linearised Einstein-Boltzmann equations. The relevant quantity to calculate is the unequal-time correlator (UETC), whose dominant eigenmodes, found by diagonalising, can be used as source functions, each individual mode being coherent [19]. The UETC Θ µν (k, τ )Θ * αβ (k, τ ) ≡ C µν,αβ (k, τ, τ ) (2.1) determines all the two-point correlation functions such as the CMB temperature C and matter power spectra P (k), defined as in [54]. Θ µν (k, τ ) is the string energymomentum tensor defined below.

A. String Energy-Momentum Tensor
Nambu-Goto strings are one-dimensional defects in the zero-width limit. They provide a good description for long cosmic strings, whose correlation length is many orders of magnitude larger than their width, at least away from string intersections. A string moving in spacetime spans a two-dimensional surface, the worldsheet x µ (σ a ), where the indices µ = 0, 1, 2, 3 label spacetime coordinates and a = 0, 1 are the indices of coordinates on the worldsheet [55,56]. The worldsheet action is reparametrisation invariant and a gauge can be chosen by imposing two conditions on the spacetime coordinates x µ as functions of σ a . In an FRW background, a useful choice of gauge is such that σ 0 = τ , the conformal time, and x ·ẋ = 0, wherė ≡ ∂/∂τ and ≡ ∂/∂σ, relabelling σ 1 , which in this gauge is a spacelike worldsheet coordinate, as σ. In this gauge the Nambu-Goto string energy-momentum tensor is Here, U is the string energy per unit length and T is the string tension. For Nambu-Goto strings on arbitrarily small scales, Lorentz invariance requires that T = U = µ. However, if we coarse-grain the string, then the integrated effect of small-scale structure is to make the effective tension smaller than the energy density. We can then include the effect of small-scale wiggles on the string via a "string wiggliness" parameter α, such that The Fourier transform of the 00-component of the energy-momentum tensor of a representative string segment in a network is where v and ξ are the string network velocity and comoving correlation length, defined in Sec. II B below, and x 0 is the position of the string endpoint. The string segment is parametrised by with the string orientations and velocity orientationŝ X is transverse toX such thatX ·Ẋ = 0. Note that the position of the string endpoint appears only through a phase in the cosine factor in equation (2.4), which we will denote as χ ≡ k · x 0 . The other components of the string energy-momentum tensor are given by with i, j = 1, 2, 3. Choosing coordinates so that k lies along thek 3 axis, the scalar, vector and tensor anisotropic stresses are given by The velocity one-scale model (VOS) equations dictate the values of the string network correlation length L, and the average velocity v, of string segments in the network [57]. The correlation length L is the average length of string segments which, for scaling networks (that have a random walk structure), is also equal to the average string separation. The network velocity v, is the root-mean-square (RMS) velocity of these correlation-length-sized string segments averaged over all (shorter) length scales. The macroscopic evolution equations for these network parameters can be derived from the Nambu-Goto action by applying a statistical averaging procedure over the string worldsheet [58][59][60]. Expressed in terms of the physical time t they read:L where a(t) is the scale factor,ȧ(t)/a(t) is the Hubble function and from now on˙≡ d/dt, unlike in equation (2.2). The loop chopping efficiency parameter c r , quantifies the energy loss due to loop production andk provides a phenomenological description of the smallscale structure on the string, which, for relativistic strings, is given bỹ (2.14) The correlation length can be written in comoving units as ξτ = L/a. The VOS equations in comoving units are where now ≡ d/dτ unlike in equation (2.2). For fixed expansion rate the scaling solutions, found by the requirement ξ = 0 and v = 0, read where β is the physical time FRW expansion exponent a(t) ∝ t β and is equal to 1/2 and 2/3 in the radiation and matter eras respectively. Note in the scaling solutions of (2.18) the implicit velocity dependence of k through equation (2.14).  both epochs [60]. Here, we also adopt this approach: at any particular τ , the values of ξ and v found using the VOS equations (2.15-2.16) are used for calculating the UETC, keeping c r constant throughout and explicitly accounting for the velocity dependence (2.14) ofk. In earlier versions of CMBACT the wiggliness α, was also an evolving parameter, but it is now kept constant in CMBACT4, which is the approach we take here. The evolution of the network parameters can be seen for a range of c r in Fig. 1 showing that a wide range of correlation lengths and velocities are available. Detailed comparison of the VOS model with Nambu-Goto simulations of ordinary string networks (i.e. single string type with unit intercommuting probability [63]) determine the loop chopping efficiency to c r = 0.23 ± 0.04 [60], corresponding to the black dot-dot-dashed curves in Fig. 1. Models of cosmic superstrings generally have suppressed intercommutation probabilities [45][46][47][48], which effectively reduces c r and so they correspond to the purple region in the figure. Such networks have relativistic RMS velocities v ∼ 1/ √ 2 and correlation lengths much smaller than the horizon, corresponding to a much higher string number density compared to ordinary string networks. However, they also have smaller string tension so their overall effect on the CMB can be small, consistent with the data.
It should be noted that the RMS network velocity used in the VOS model arises from a worldsheet average and is thus integrated over all (short) length scales. Therefore, it provides an accurate measure of the energy stored in a wiggly string segment, but does not explicitly correspond to (and in fact is expected to be larger than) the coherent velocity on correlation-length scales. Indeed, numerical simulations of Nambu-Goto strings reveal a network veloc-ity distribution with larger velocities at short scales, implying that the RMS velocity is dominated by relativistic speeds at short distances. On length scales of order the correlation length, coherent velocities as low as v coh 0.2 have been reported [64][65][66][67]. Other network velocity measures (again containing information from a range of length scales) in both Nambu-Goto and Abelian-Higgs string simulations also tend to be lower than the VOS RMS velocity, with velocities in the Abelian-Higgs model v AH 0.5, significanlty slower than in Nambu-Goto simulations [68][69][70]. For further discussion about the impact of string velocities on the UETC and the string power spectrum see the end of Sec. II F.

C. Unconnected Segment Model
Simulations of evolving string networks are numerically very expensive. Strings decay as 1/(ξτ ) 3 , eventually reaching a scaling solution (ξ = constant) with a number density of tens to hundreds of strings per horizon volume. At early times, the box contains a huge number of strings whose dynamics and interactions have to be tracked at each time step. The unconnected segment model (USM) [21,61] dramatically reduces the required computational resources by approximating the string network as a collection of correlation-length-sized segments, with the time evolution of the correlation length and segment velocity described by the VOS equations. Moreover, the model consolidates these string segments by collecting all strings that decay between any two times, and so fewer strings need to be tracked. The number of strings that decay between any two conformal times in a volume V , is where n(τ ) is the number density of strings at conformal time τ , given by n(τ ) = C(τ )/(ξτ ) 3 . In CMBACT, the factor C(τ ) is chosen so as to keep the number of strings at any time proportional to 1/(ξτ ) 3 . The energy-momentum tensor for the string network is then given by the sum over the total number of consolidated string segments K , with a factor accounting for string decay The string decay factor T off (τ, τ i , L f ) is a function interpolating between 1 and 0 and is responsible for turning off the contribution of the i th consolidated segment after the time it has decayed. Its steepness is controlled by a string decay parameter 0 < L f ≤ 1, as follows: Thus, in the limit L f → 1 the string decay factor T off (τ, τ i , L f ) approaches a Heaviside function, sharply switching off the contribution of the i th consolidated segment to the network energy-momentum tensor for times τ > τ i .
Since the number of consolidated segments also sets the number of decay epochs, a finite number of consolidated segments leads to discrete steps in the number density of strings. The string decay parameter L f was introduced to allow a fraction of the consolidated strings to decay before the end of their respective decay epoch, thus making the number density evolution smoother. The function C(τ ) was also introduced to ensure that the number of strings at any conformal time τ is kept proportional to (ξτ ) −3 . However, one consequence of L f < 1 is that it is possible that L f τ i+1 < τ i , meaning strings can start to decay earlier than their respective epoch and the number density is systematically lower.
In the CMBACT4 implementation we have found that changing the number of consolidated segments from 200 to 10000 has very little impact on the string spectra, as shown in Fig. 2. However, the amplitude of the C is dependent on the value of L f . The change is scale dependent, but can be as much as 30%, for example near the peak of the scalar temperature signal. Previous analyses which have used the results from CMBACT have overlooked this dependence. Although not entirely degenerate with the amplitude of C , which scales proportional to (Gµ) 2 , it will clearly have some affect on the inferred values of Gµ from the USM. We compare this to our approach in the following section.

Infinite Consolidated String Segments
We are able to accommodate a large number of segments analytically. As discussed in [53], the scaling factor, that weights the UETC taking into account string decay, has a particularly simple form when the number of consolidated string segments tends to infinity, L f → 1 and C(τ ) → 1. This is An analytic expression for the scaling factor can also be found for arbitrary L f using the form of T off quoted in equation (2.21). However, it seems natural to consider only the case L f = 1 when the number of consolidated string segments is very large. In the infinite limit the segments will decay at an infinite number of epochs which are infinitesimally separated, a   tensor C modes respectively. The first column contains the temperature (TT) C , the second column has the EE mode contribution, BB modes are in the third and the TE cross-correlation in the final column. We also plot the corresponding spectra derived from our analytic USM method, shown in green dot-dot-dashed lines.
continuous limit in which the string number density is smooth. We have shown that the number density scales according to (ξτ ) −3 with our approach. While infinite consolidated segments may seem unphysical, it is just a limit used to obtain the correct scaling relation. We obtain very similar results to CMBACT4 when using between 200 to 10000 segments with L f = 1. The question of whether the observed resulting modification of scaling from early string decay obtained when L f < 1 is physical or not requires investigation.
Since we takE C(τ ) = 1 we avoid considering different scaling behaviour. Ultimately, the USM is a simplified model which aims to match the UETC from simulations by adjusting the network parameters. Overall it has been shown to match Nambu-Goto simulations well [71]. However, due to the correlation between the inferred values for Gµ for a given L f , this issue should be considered more closely.
Since the number density scales according to (ξτ ) −3 using our approach, we believe this to be reasonable and will adopt this for the comparison to data.

D. Analytic Calculation of the Unequal-Time Correlator
The UETC can be computed analytically [53] by integrating over all string configurations (orientations and positions) in the network. For the two point correlator between Θ(τ 1 , k 1 , χ 1 ) and Θ(τ 2 , k 2 , χ 2 ) translational invariance implies k 1 = −k 2 = k and so χ 1 = −χ 2 = χ. Considering that, due to equations (2.4) and (2.8), Θ(τ, k, χ) is a symmetric function of k the integral is Without loss of generality k can be chosen to lie along the k 3 -axis, such that k = kk 3 . Θ here represents each of Θ 00 , Θ S , Θ V and Θ T of equations (2.9-2.11). The φ, ψ and χ integrals can be done analytically in this case leaving only the θ integral in terms of Bessel functions. The UETC can then be written as the sum over six integral identities respectively. This extends the corresponding result of [53] in that ξ and v are now functions of τ instead of being kept constant. This means that the expressions of the amplitudes A i , presented in Table II, are now time-dependent. The integral identities (shown in Table I) remain the same. It should be noted that I 1 (x, ) and I 4 (x, ) diverge but the combination I 1,4 (x − , ) − I 1,4 (x + , ) is regular and, in the limit where x 1,2 x 2,1 , has an analytic approximation given by In the small x 1,2 , limit the UETC can be written as (2.28) and at equal times, when x 1 = x 2 = x and = 0, the equal-time correlator is given by The form of B and C are similar to [53] but again depend on the values of v and ξ at τ 1 and τ 2 . These coefficients have also been included in Table II. Thanks to these analytic approximations, computational times can be greatly reduced compared to the case where the integral identities I i are used for computation over the whole range of kτ 1 , kτ 2 . The regions where these approximations are valid are shown in Fig. 3, only the white region is computationally intensive. It should be noted that, because ξ is a function of time, the shape of the approximated regions in Fig. 3 changes for different values of k and so we must consider a large number of k-modes when computing the UETC. This is in contrast to [53], where the approximation of constant ξ and v meant that the UETC was only a function of the combinations kτ 1 and kτ 2 .

Negative values of the UETC
It has been noted in [72] that there are negative regions in the string UETC calculated analytically through our formalism, which do not appear in the Gaussian model for the string UETC used in [72]. These can be seen in Fig. 4.
There are two distinct types of regions with negative values of our UETC. First, regions with small kτ 1 and large kτ 2 (and vice versa), corresponding to the top left and bottom right corners of Fig. 3 or Fig. 4: in these regions the UETC should be zero, but small negative (and positive) values can arise from the finite order truncation of the Bessel series expansions of I 1 (x ± , ρ) and I 4 (x ± , ρ) in Eq. (2.25). These values are spurious and can be thought of as noise arising from the truncation. The order of truncation must then be chosen such that this noise is at a tolerable 1, red when | log x1 − log x2| < and blue when |x1 − x2| 1. In the code the x1,2 1 region is set for x1,2 < 0.2, = 0.001 for x1 ≈ x2 and |x1 − x2| > 10 for x1,2 x2,1. level.
Second, in the regions off the diagonal with large kτ 1 ≈ kτ 2 (corresponding to the top right corner of Fig. 3 or Fig. 4) there is a ringing pattern with successive positive and negative peaks that decay as we move away from the diagonal. These oscillatory patterns are a consequence of causality [21,73,74], built into the USM: as the correlator must vanish at superhorizon scales (in fact in the USM it vanishes at scales larger than the correlation length, which is smaller than the horizon), this introduces a sharp edge in physical space that becomes oscillatory in Fourier space. This oscillatory pattern therefore has a clear physical origin, but in the USM it is somewhat artificially enhanced due to the fact that the model assumes all string segments have the same length. If segments are instead given a length distribution peaking at the network correlation length, the sharp edge is smoothed and the oscillatory pattern gets suppressed. Further, considering a segment velocity distribution peaking near the network RMS velocity again suppresses these oscillations. The Gaussian model assumes a wide Gaussian distribution of string lengths (but also assigns non-zero values to the correlator at superhorizon scales) so this causal oscillatory feature is absent from the UETC in that model.
The suppression of oscillations in the UETC can be seen in Fig. 5 where the blue solid line shows the profile of the UETC across the diagonal as calculated using the velocity and correlation lengths from VOS. In red dot-dot-dash is the same profile when a Gaussian distributed sample of velocities and correlation lengths, peaking on the VOS values, are chosen. The oscillatory features are mostly washed out but the first trough remains a prominent feature. The off-diagonal dip in the correlation functions that we find after considering a range of segment lengths and velocities has also been observed in Abelian-Higgs simulations [69].
It may also be related to the velocity anti-correlation observed in Nambu-Goto simulations on correlationlength scales and can be attributed to string intercommutations [66].

E. Eigenmode Decomposition
The UETC is generally rescaled by a factor of √ τ 1 τ 2 , which, for ξ and v constant, makes it a function of kτ 1 and kτ 2 only. This is not true in the present case because now we are tracking the time-dependence of ξ and v, so the UETC depends separately on k, τ 1 and τ 2 . However, it is still useful to introduce this rescaling in order to facilitate direct comparison of the UETC with previous results. This rescaled UETC can then be discretised onto a logarithmic grid in kτ 1 and kτ 2 with n × n grid points and then diagonalised giving the eigenvectors and eigenvalues [19] Due to the explicit dependence on k, this diagonalisation procedure has to be repeated for a large number of k-modes, and the eigenvalues are k-dependent. This significantly increases the computation time compared to [53]. The extra factor (k 2 τ 1 τ 2 ) γ is used for more efficient reconstruction of the UETC when the eigenmodes are truncated below n. The choice γ = 0.25 gives the best reconstruction on scales that give the dominant contribution to the CMB anisotropies.
There is no correlation between the scalar, vector and tensor modes and so the vector and tensor UETC can be diagonalised independently. However, the density Θ 00 , and scalar anisotropic stress Θ S , are correlated. The diagonalisation is done over a 2n × 2n grid constructed from is the symmetric combination of the cross-correlation between Θ 00 and Θ S . After diagonalisation, the first half of the eigenvectors refer to the density and the second to the anisotropic stress. The diagonalisation creates orthogonal eigenvectors which are then used as source terms in the CAMB [75] linear Einstein-Boltzmann code. The C are calculated using each individual eigenvector u i (kτ )/( √ τ (kτ ) γ ), as a source function C i , which can be summed to give the total power spectra (2.32) By ordering λ i from largest to smallest, the required accuracy in the C can be achieved by including relatively few eigenmodes. This can be seen in the middle row of Fig. 6 where there is only about 10% difference between using all 512 eigenmodes of a 512 × 512 grid  compared to only using 32 eigenmodes when fixing the value of Gµ. Also, it can be seen in the top row of Fig. 6 that reducing the grid resolution reduces the amplitude of the C . A grid resolution of 128 × 128 is about 5% lower, on average, than using the 512 × 512 grid but convergence times decrease drastically. It should be noted that there is negligible difference between using a 512×512 and a 1024×1024 grid meaning that the former is reliably giving the full C contribution. The bottom row shows what happens when using more k values in the calculation. Wiggly features arise from using too few k values and can be removed at the expense of a much longer calculation. Using these findings we can choose the optimal UETC parameters to give good quality C in a reasonable amount of time. The resulting spectra obtained from our analytical method are shown in Fig. 2 in green dot-dot-dashed curves and agree well with USM string realisations, especially in the limit of large numbers of simulated segments.

F. Comparison of the String Power Spectrum
In Fig. 7 we compare our temperature power spectrum (scaled by Gµ in the upper subplot and normalised at = 10) to that of CMBACT4 [61], Nambu-Goto simulations [71], and Abelian-Higgs simulations [69]. Both CMBACT4 and our method use the same velocity dependent one-scale model parameters, but CMBACT4 uses L f = 0.5. The Nambu-Goto simulations are performed in an expanding background from recombination to today, including Λ domination. Large loops are kept in the simulation and contribute to the total energy-momentum tensor of the network, but these simulations cannot resolve smallscale physics near the string width and do not include the effects of radiation backreaction. In contrast, the Abelian-Higgs simulations can resolve small-scale structure and radiative effects [76]. These, however, have smaller dynamical range and cannot easily evolve through the radiation-matter transition (so the UETC is instead interpolated), but see recent progress in [76] where the authors simulate through the transition.
Overall, when normalised at = 10, the four spectra agree reasonably well. The USM variants (CMBACT4 and our approach) both predict slightly more power at the peak than either of the simulations. The Nambu-Goto simulations predict more power on very small scales, around twice as much as the Abelian-Higgs model. It is well known that Nambu-Goto calculations yield higher string densities than field theoretic ones, which will increase their overall normalisation. The resulting constraints on Gµ are therefore around a factor of 50% lower [33] as can be inferred from the upper subplot in Fig. 7. The USM variants are closer to the Nambu-Goto simulations in this respect [71]. Within this paper we will not consider using the analytic USM to mimic the Abelian-Higgs spectra. As we have shown, there is some additional uncertainty in the USM, as the normalisation depends somewhat on the choice of L f .
In summary, given the large differences in modelling between the various approaches we find this comparison encouraging, although more work is needed to further delineate the differences. In particular, as discussed at the end of Sec. II B, the VOS RMS velocity is defined through a worldsheet integral over all scales and receives a large contribution from relativistic wiggles on the string. On the other hand, the USM assumes straight segments moving at a given speed and the small-scale structure on the segments is captured via a "renormalisation" of their tension. This implies that the speed to be associated to the USM segments must be lower than the VOS RMS velocity, and should correspond to the network velocity at correlation length scales. Numerical simulations show this to be significantly lower than the RMS speed. This issue has not been examined before, partly because the calculated string spectra from different approaches can differ by up to a factor of two, and partly because it can be offset by choosing a lower value for the USM parameter L f (see below). As quantitative agreement between the different approaches is now being established, it is important to fully understand this issue. To this end it will be important to extract the network velocity distribution as a function of length scale in both Nambu-Goto and Abelian-Higgs simulations. Fig. 7 in purple dot-dot-dash is the C obtained when v = 0.4. As can be seen, the peak of the velocity fixed C has a very similar amplitude to the Nambu-Goto simulation C in dotted red, although the simulations still have larger power at both lower and higher . This supports the idea that the discrepancy in the amplitude of string spectra could be related to different predictions/assumptions on string velocity in the different approaches (cf. discussion at the end of Sec. II B). Note that the parameter L f in the USM is somewhat degenerate with the string velocity -for fixed v a lower L f reduces the density of   strings by increasing the string decay rate, thus reducing the C amplitude and matching simulations better than using L f = 1. In the absence of a more complete quantitative understanding of the string velocity distribution -input required from string evolution simulations -our string spectra obtained from the USM have a larger amplitude (see the solid blue line in the upper subplot of Fig. 7). This leads to slightly tighter constraints on cosmic strings than in numerical simulations. Marginalising over the network parameters c r and α, partly takes care of the differences between L f = 0.5 and L f = 1 in the USM since high c r reduces the velocity (as seen from equation (2.16) and pictorially in Fig. 1).

III. COSMIC SUPERSTRINGS
A cosmic superstring network can be modelled as a collection of string segments of different types, each string type having its own tension and self-intercommuting probability [44-50, 52, 77-79]. Strings of different types interact with each other via "zipping" or "unzipping" leading to heavier or lighter strings respectively, that are connected to the original strings at trilinear Y-shaped junctions [80]. The fundamental building blocks for these networks are light (fundamental) F-strings and heavier (Dirichlet) D-strings, with a tension hierarchy controlled by the fundamental string coupling [80][81][82]. Heavier strings arise as bound states between p F-strings and q Dstrings, where p,q are coprime. Given the fundamental string tension, the corresponding tensions of these heavier (p, q)-strings are controlled mainly by p,q and the value of the string coupling. These networks generally behave very differently than their ordinary cosmic string counterparts. They are typically characterised by small intercommutation probabilities, thus leading to higher string number densities [44,45,49,52]. The complex interactions present imply that several string types with different tensions and correlation lengths can simultaneously contribute to the string network CMB spectra.
In scaling superstring networks, the string number density is dominated by the lightest F-strings, followed by D-strings and the first bound state, i.e. (1,1)strings. Heavier bound states are suppressed, so the number of string types considered in the model can be truncated at a finite number. Following [52] we shall describe the network by keeping seven distinct types of strings: where the last column describes the (p, q) charges of the corresponding string type.
The large-scale dynamics is then modelled by seven copies of the VOS equations, appropriately extended to account for transfer of energy among the different string types through zipping and unzipping interactions [44,50]. In each copy of the VOS equations describing a single string, say of type i, the self interaction coefficient c r in equation (2.15) is replaced by the corresponding self-interaction coefficient c i , and new cross-interaction terms with coefficients d k ij are added to describe zipping and unzipping. The coefficients c i , d k ij are controlled by the corresponding microphysical intercommuting probabilities P ij [52], which can be estimated [46,48] from the corresponding string theoretic amplitudes (and field theory approximations in the case of non-perturbative interactions between heavy strings [47]). They can be expressed as a product of two pieces: one that is dependent on the volume of the compact extra dimensions V ij (w, g s ), and a quantum interaction piece F ij (v, θ, g s ). Physically, one can think of V ij as arising from string position fluctuations around the minimum of a localising potential well, giving rise to an effective volume seen by each type of string. The heavier the string the smaller the fluctuations are and so the smaller the value of V ij [46]. The parameter w corresponds to the effective volume in the compact extra dimensions seen by F-strings. g s is the fundamental string coupling and v and θ are the relative velocity and angle of the incoming strings. For a pair of strings colliding at an angle θ, and relative speed v, the intercommuting probability is (3.2) Details of how F ij and V ij are calculated can be found in [52]. Since the network contains a large number of individual strings with a range of velocities and orientations, the coefficients c i and d k ij are determined by the integral of P ij over a Gaussian velocity distribution centred on the scaling network velocities of each string type and over all angles. This gives the average intercommuting probabilities P ij (w, g s ) ≡ P ij . Numerical simulations of single-type Nambu-Goto strings with small intercommuting probability [49] suggest that the self-interaction coefficients c i scale as: where c s is the standard self-interaction coefficient in three dimensions corresponding to the value c r in Sec. II B. This choice of c s implies a convenient normalisation of the coefficients c i so that one recovers the ordinary cosmic string value c r when P ii = 1. This facilitates direct comparison with ordinary cosmic strings.
For cross-interactions between two strings of types i and j (i = j), producing a segment of type k, there is an additional factor arising from the kinematic constraints of Y-junction formation [83,84] that we denote as S k ij (i = j). This also arises as an integral over relative velocities and string orientations [52,85]: where S is a normalisation factor [52], Θ(−f→ µ (v, θ)) imposes the kinematic constraints [84] and σ 2 v is the variance of the velocity distribution peaked on the relative scaling velocitiesv ij = (v 2 i + v 2 j ) 1/2 between strings of type i and j. The cross-interaction coefficients are then given by ij . The overall normalisation κ is the analogue of c s , but for cross-interactions. There is no obvious choice for this phenomenological parameter, but it may be expected to be of order unity by analogy to the ordinary self-interacting string result for c r , obtained by numerical simulations. Strictly speaking it should be treated as an extra parameter for the model but, given the large computational resources required in our MCMC analysis, we will set it to unity in this work. Our analysis will still indirectly capture the effects of changing this parameter as it is somewhat degenerate with w. To see this, note that d ij is also proportional to P 1/3 ij which depends weakly on w through the volume factor V ij (w, g s ).
The leading effect of w is to change the relative amplitude between self-interactions (FF interactions having the strongest w dependence) and cross-interactions of heavy strings, thus mimicking somewhat the effect of varying κ relative to c s . As computational power improves and our methodology is refined, κ should be re-introduced as an additional MCMC parameter.
The modified VOS equations [50,52], in comoving units, are where i ab is the average length of segments of type i formed by the zipping/unzipping of string types a and b at conformal time τ , and µ i is the tension of the i th string type. All string tensions can be expressed in terms of the fundamental string tension µ F , and in flat spacetime [80][81][82] are given by: where p i and q i are the charges of string type i as listed in (3.1). The coefficients b i ab appearing in the velocity evolution equations (3.7) are related to energy conservation and allow for the energy saved from zipping interactions to be redistributed as kinetic energy of the new segment (b i ab = d i ab ) [50] or radiated away (b i ab = 0) as in [44]. A more realistic model should have a specific radiation mechanism so that 0 < b i ab < d i ab such that some of the energy is redistributed whilst the rest is radiated away. However, for cosmic superstring networks (for which d ij are much smaller than unity) this term has negligible impact on the string scaling densities and velocities [52,85], so here we take b i ab = 0.
Once the velocities and correlation lengths of all string types in the network are obtained by solving (3.6-3.7), their unequal-time correlators can be calculated independently as laid out in Sec. II. Although N > 3 string types are needed in order to accurately construct the abundances of the dominant three lighter strings (in this case seven string types are used (3.1)), the resulting scaling densities of the higher charged states with N > 3 are strongly suppressed compared to the lighter F-, D-and FDstrings [44,50,85]. This allows us to only consider these first three states in the computation of CMB signatures through our UETC analytic method. The evolution of the network parameters for the three lightest strings can be seen in Fig. 8 for c s = 0.23, w = 1 and g s = 0.3.
Once the UETC of each of the three lighter strings are calculated they can simply be summed to give the total string UETC, since the individual segments are uncorrelated in the USM. This can then be diagonalised and the eigenvectors and eigenmodes used as sources for finding the contribution from cosmic superstrings to the CMB anisotropy. We have checked that our analytic UETC method reproduces the results of Fig. 4 in [52], including the shift in the location of the peak as we vary g s . We have found a slightly lower amplitude in the B-mode spectrum that can be attributed to the extra factor of 2 in the vector modes that was present in CMBACT3 (which [52] was based on) and has been corrected in CMBACT4 [86].

IV. STRING CONSTRAINTS
We obtain joint constraints on cosmic string network and ΛCDM parameters using a modified version of COSMOMC. To reduce computational time in our analysis we have tested two methods for deriving string network constraints. In the first method, the string C are pre-calculated for a range of c r = [0.1, 1] and α = [1,10] at the Planck best fit values for the cosmological parameters, i.e. Ω b h 2 , Ω c h 2 and H 0 . These C are read into COSMOMC, interpolated at the MCMC c r and α values and then scaled by (Gµ) 2 . This is an extremely efficient way for obtaining network constraints since only the ΛCDM C need to be calculated, while the interpolation takes very little time. We have checked that the difference in the resulting string C when calculated at the upper and lower 3σ bounds in Ω b h 2 , Ω c h 2 and H 0 is ∼ 0.5% in the temperature, E-and B-modes and no more than ∼ 10% in the TE cross-correlation. This uncertainty in the string C is 1% of the total C . The C for different c r and α are plotted in Fig. 9. The different bands of colour indicate the value of c r , solid red being the lowest (c r = 0.1) then progressing through long-dashed yellow, short-dashed green, dot-dashed blue and dotdash-dotted purple in steps of 0.2, up to c r = 0.9. The upper (patterned) and lower (dot-patterned) edges of the bands indicate α = 10 and α = 1 respectively. From this it can be seen that the effect of α is to change the amplitude of the C , with lower α also flattening the small features (as best seen in the second column and to a lesser extent in the third column of Fig. 9). Increasing c r reduces the amplitude of the C and, as best seen in the third column of Fig. 9, shifts the main peak towards slightly smaller . In the second method, which is computationally expensive, we simply calculate the string and ΛCDM C for each (network) parameter value and compare to CMB data.
The same process of pre-calculating string spectra can be done for cosmic superstring networks in the parameter ranges c s = [0.1, 1], g s = [0.01, 0.9] and w = [0.001, 1]. The superstring C can be seen in Fig. 10, where same colours and patterns are used for the steps in c s as in Fig. 9. The bands indicate values of w, with w = 10 −3 corresponding to the solidpatterned lines and w = 1 to the dotted version of the same pattern. The rows indicate varying values of g s , with g s = 0.01, g s = 0.1 and g s = 0.9 for the top, middle and bottom rows respectively. The first point to notice is that the C amplitudes at low g s are much greater than those at large g s . For large c s values there is less difference between the greatest and smallest values of w, especially at low g s , i.e. the purple dotdash-dotted lines in the top row of Fig. 10 overlap, but are well separated in the bottom row. This is because for large c s the cross-interaction terms d k ij (which are less dependent on w than the self-interaction terms c i ) play a more important role in setting the scaling string number densities. For small values of c s , the c i coefficients become smaller (while d k ij are unaffected) leading to small correlation lengths and so large string number densities. The C amplitudes are then affected more strongly by c i , giving rise to a stronger dependence on w.
The datasets used in the MCMC analysis come from the Planck2015 mission [41], in particular: Planck2015 TT+lowP: This contains the 100-GHz, 143-GHz, and 217-GHz binned half-mission TT crossspectra for = 30 − 2508 with CMB-cleaned 353-GHz map, CO emission maps, and Planck catalogues for the masks and 545-GHz maps for the dust residual contamination template. It also uses the joint temperature, E and B cross-spectra for = 2 − 29 with E and B maps from the 70-GHz LFI full mission data and foreground contamination determined by 30-GHz LFI and 353-GHz HFI maps. Planck2015 TT+Pol+lowP: This contains the same data as Planck2015 TT+lowP but also uses the TE and EE cross-spectra for = 30 − 1996. Planck2015 TT+Pol+lowP+BKPlanck: This again contains all of the data used in Planck2015 TT+Pol+lowP but includes also the cross-frequency spectra between BICEP2/Keck maps at 150 GHz and Planck maps at 353 GHz including the B-mode spectra at multipoles ∼ 50 − 250.
We first consider our interpolation method, where the C are pre-calculated on a grid in c r and α (or in the case of cosmic superstring networks c s , g s and w), and then a spline interpolation used between grid values. The results obtained from this method are very quick and accurate due to our ability to use all 512 eigenmodes of the 512 × 512 grid for the UETC. The constraints on network parameters derived from this method are shown in Fig. 11. Gµ is implemented into the MCMC analysis through a logarithmic prior of [−10, −5] such that Gµ = 10 [− 10,−5] .
There is no significant difference in our constraints when using Planck2015 TT+lowP, or including EE and TE or both EE and TE and BB results. The upper 2σ value for the tension is Gµ < 1.1 × 10 −7 for Planck2015 TT and is similarly Gµ < 9.6 × 10 −8 and Gµ < 8.9 × 10 −8 for Planck2015 TT+Pol+lowP and Planck2015 TT+Pol+lowP+BKPlanck. These agree well with the Gµ < 1.8 × 10 −7 and Gµ < 1.3 × 10 −7 from the Planck cosmological parameters analysis [33]. The slightly tighter constraints obtained here are due to the amplitude of the C not scaling with the value of L f , i.e. the C are larger when L f = 1 as assumed here, while previous results were obtained from CMBACT with L f = 0.5. There is little difference between using the Planck temperature data alone and including polarisation data as expected from [33]. As can be seen in the other two columns of Fig. 11, c r and α are not constrained. There is a slight preference for higher values of c r and lower values of α since both of these lead to smaller C . Features, such as the position of the main peak or the pronounced lower peak make very little difference to the overall constraints. There is a very slight correlation between Gµ and c r and anticorrelation between Gµ and α, as expected from the C seen in Fig. 9. A combination of high α and low c r is mildly disfavoured. Further, by comparing the constraints on Gµ and c r to their affect on the C in Fig. 9 there is a larger difference between changes at small c r than changes at large c r . For this reason we expect to see greater correlation between Gµ and c r on a logarithmic scale from values c r 1 to c r ≈ 0.1 than implied over our prior range.
Considering our direct calculation method, where the string spectra are calculated every time along with the C from ΛCDM, the constraints are slightly weaker. This is because there is a pay-off between the resolution of the UETC and number of eigenmodes used in the reconstruction and the time spent computing the spectra. To efficiently calculate the constraints a grid resolution of 128 × 128 with 64 eigenmodes has been used. As can be seen in Fig. 6 we expect a reduction in power of about 10 − 20% which means the value of Gµ is allowed to be higher than when the high resolution, full reconstruction interpolation method is used. For Planck2015 TT+lowP this is Gµ < 4.3 × 10 −7 . The constraints on c r and α also show a slight preference for lower c r and larger α, as in our interpolation method.
For cosmic superstrings, Gµ F , g s and w are marginalised over logarithmic priors, and c s over a flat prior. Again all 512 eigenmodes of the 512 × 512 grid for the UETC are used. The likelihood contours obtained from our interpolation method can              be found in Fig. 12. It can be seen that w and c s are almost flat (columns 3 and 4), again with larger values of c s favoured as this leads to smaller amplitude C . As the string density grows with decreasing g s , the constraints on g s favour larger values, as seen in the second column. Note, however, that the model is not reliable for large values of g s as the perturbative expansion starts to break down and the string interaction amplitudes used in c i and d k ij have large uncertainties. Finally, the first column shows our constraints on the fundamental string tension Gµ F , which is much smaller than for ordinary cosmic strings. We find Gµ F < 2.8×10 −8 for Planck2015 TT+lowP when marginalising over g s , c s and w, and the same constraint for Planck2015 TT+Pol+lowP and Planck2015 TT+Pol+lowP+BKPlanck. Fig. 12 we show the constraints when using the direct calculation method, where the string spectra are calculated at every step in the Markov chain. This is a much more intensive computation and so a lower resolution grid and fewer eigenmodes in the reconstruction had to be used. As for cosmic strings the optimal balance between computing time and accuracy suggested using a 128 × 128 grid with 64 eigenmodes. The constraints are thus slightly weaker, with the main result Gµ F < 4.2 × 10 −8 . The results from our two methods are in good agreement, justifying the use of our interpolation method, and showing that varying ΛCDM parameters within Planck priors has little effect on the string constraints.

V. CONCLUSIONS
Currently, there are two main approaches to the detection of cosmic strings. Firstly, since they actively generate scalar, vector and tensor perturbations they lead to signatures in the temperature, polarisation, and non-Gaussian spectra of the CMB. Secondly, a cosmic string network will emit gravitational waves, primarily from loop decay. This leads to a stochas- tic background which can be constrained using pulsar timing, laser interferometry experiments such as LIGO and eLISA, and also the CMB [87]. A transient gravitational wave signal is also expected from cusps and kinks in the network [88]. The latter class of tests has the potential to provide even stronger constraints on the string tension Gµ, but there are large uncertainties in the loop size, which is fixed by gravitational back-reaction. Model dependence on gravitational waves from cosmic strings further makes it difficult to determine signatures, for example, whilst Nambu-Goto strings decay into loops, Abelian-Higgs strings primarily decay into particles [88][89][90]. It is therefore important to use a variety of complementary observational probes.
The first class of tests also suffer from uncertainties, but these are less significant. The string UETC can be obtained from simulations and used as source functions in CMB codes, but simulations are numerically expensive and suffer from issues in dynamical range. An alternative approach is to model the string network as an ensemble of segments using the USM. Crucially, although the USM provides a simplified picture of the network, it is able to match simulations by adjusting the free parameters of the model, namely the correlation length, RMS velocity and string wig-gliness.
In this paper we have significantly improved and extended our previous work on string power spectra from the USM: 1. We have analytically solved the UETC for an evolving string network, whereas our previous work was restricted to constant network parameters. The UETC itself can be computed in under a minute. For the CMB power spectrum, although the time taken is increased due to tracking a larger number of Fourier modes, on a 3.1 GHz Intel Xeon CPU with 8 threads, our code runs in ∼ 60 minutes. For comparison, around 2000 network realisations are required for CMBACT4 to achieve the same accuracy and since this code is serial, the computation time is ∼ 30 hours.
2. We have extended the formalism to cosmic superstring networks with multiple string types and different network parameters. Here the UETC can be computed for each string type and added, since the segments are assumed to be uncorrelated. The UETC calculation is much quicker than the CMB line-of-sight integration, so the total computation time is not significantly increased over the single string case. 3. For the first time we have been able to marginalise over the string network parameters when fitting to Planck2015 and joint Planck-BICEP2 data. The data is consistent with no strings for both the single and multi-string case. Since other network parameters are unconstrained when the tension is very small, it is only possible to present joint constraints on these with Gµ. In the superstring case, for example, the constraint on the string coupling g s is degenerate with Gµ F .
There are several possibilities to explore in future work. Firstly, there are various ways in which the USM could be improved. Superstring networks contain Y-type junctions, but in the present formulation these only impact the evolution of the network parameters. Since junctions are relatively rare in the limit of large and small coupling, the USM is expected to provide a sufficient description. However, in some regimes the energy density of the network may not be dominated by a single string type, and junctions may become important. In this case the USM could be modified to include a correlation between segments. A further improvement is the inclusion of loops. The decay of string segments in the USM should mimic the energy loss in loops, but it is possible these may lead to additional interesting signatures.
Given that Planck has largely exhausted the available signal in the temperature data, future string constraints from the CMB will be driven by polarisation and non-Gaussianity. The non-Gaussian signal from post-recombination simulations has been used to obtain constraints on Gµ [32], and attempts have been made to compute the bi-spectrum analytically using a Gaussian model for the string correlators [91]. It is also possible to compute the non-Gaussian signal using the USM which will, by design, include physics from recombination and along the line-of-sight. This has already been demonstrated for the CMB bi-spectrum [92] by performing many realisations of the network. It is possible to employ a similar analytic method used in this work to compute the string bi-and tri-spectrum, which we would expect to be significantly faster [93].
The detection of gravitational waves by LIGO is particularly exciting for strings, and the next generation of ground and space based experiments can potentially provide much stronger limits than those from the CMB. However, these limits strongly depend on modelling, for example, the loop, kink and cusp distribution. Further work is needed to understand these and until then the CMB will continue to be an important tool in the search for strings.