|
| 1 | +# evolve the RGEs of the standard model from electroweak scale up |
| 2 | +# by dpgeorge |
| 3 | + |
| 4 | +import math |
| 5 | + |
| 6 | +class RungeKutta(object): |
| 7 | + def __init__(self, functions, initConditions, t0, dh, save=True): |
| 8 | + self.Trajectory, self.save = [[t0] + initConditions], save |
| 9 | + self.functions = [lambda *args: 1.0] + list(functions) |
| 10 | + self.N, self.dh = len(self.functions), dh |
| 11 | + self.coeff = [1.0/6.0, 2.0/6.0, 2.0/6.0, 1.0/6.0] |
| 12 | + self.InArgCoeff = [0.0, 0.5, 0.5, 1.0] |
| 13 | + |
| 14 | + def iterate(self): |
| 15 | + step = self.Trajectory[-1][:] |
| 16 | + istep, iac = step[:], self.InArgCoeff |
| 17 | + k, ktmp = self.N * [0.0], self.N * [0.0] |
| 18 | + for ic, c in enumerate(self.coeff): |
| 19 | + for if_, f in enumerate(self.functions): |
| 20 | + arguments = [ (x + k[i]*iac[ic]) for i, x in enumerate(istep)] |
| 21 | + try: |
| 22 | + feval = f(*arguments) |
| 23 | + except OverflowError: |
| 24 | + return False |
| 25 | + if abs(feval) > 1e2: # stop integrating |
| 26 | + return False |
| 27 | + ktmp[if_] = self.dh * feval |
| 28 | + k = ktmp[:] |
| 29 | + step = [s + c*k[ik] for ik,s in enumerate(step)] |
| 30 | + if self.save: |
| 31 | + self.Trajectory += [step] |
| 32 | + else: |
| 33 | + self.Trajectory = [step] |
| 34 | + return True |
| 35 | + |
| 36 | + def solve(self, finishtime): |
| 37 | + while self.Trajectory[-1][0] < finishtime: |
| 38 | + if not self.iterate(): |
| 39 | + break |
| 40 | + |
| 41 | + def solveNSteps(self, nSteps): |
| 42 | + for i in range(nSteps): |
| 43 | + if not self.iterate(): |
| 44 | + break |
| 45 | + |
| 46 | + def series(self): |
| 47 | + return zip(*self.Trajectory) |
| 48 | + |
| 49 | +# 1-loop RGES for the main parameters of the SM |
| 50 | +# couplings are: g1, g2, g3 of U(1), SU(2), SU(3); yt (top Yukawa), lambda (Higgs quartic) |
| 51 | +# see arxiv.org/abs/0812.4950, eqs 10-15 |
| 52 | +sysSM = ( |
| 53 | + lambda *a: 41.0 / 96.0 / math.pi**2 * a[1]**3, # g1 |
| 54 | + lambda *a: -19.0 / 96.0 / math.pi**2 * a[2]**3, # g2 |
| 55 | + lambda *a: -42.0 / 96.0 / math.pi**2 * a[3]**3, # g3 |
| 56 | + lambda *a: 1.0 / 16.0 / math.pi**2 * (9.0 / 2.0 * a[4]**3 - 8.0 * a[3]**2 * a[4] - 9.0 / 4.0 * a[2]**2 * a[4] - 17.0 / 12.0 * a[1]**2 * a[4]), # yt |
| 57 | + lambda *a: 1.0 / 16.0 / math.pi**2 * (24.0 * a[5]**2 + 12.0 * a[4]**2 * a[5] - 9.0 * a[5] * (a[2]**2 + 1.0 / 3.0 * a[1]**2) - 6.0 * a[4]**4 + 9.0 / 8.0 * a[2]**4 + 3.0 / 8.0 * a[1]**4 + 3.0 / 4.0 * a[2]**2 * a[1]**2), # lambda |
| 58 | +) |
| 59 | + |
| 60 | +def drange(start, stop, step): |
| 61 | + r = start |
| 62 | + while r < stop: |
| 63 | + yield r |
| 64 | + r += step |
| 65 | + |
| 66 | +def phaseDiagram(system, trajStart, trajPlot, h=0.1, tend=1.0, range=1.0): |
| 67 | + tstart = 0.0 |
| 68 | + for i in drange(0, range, 0.1 * range): |
| 69 | + for j in drange(0, range, 0.1 * range): |
| 70 | + rk = RungeKutta(system, trajStart(i, j), tstart, h) |
| 71 | + rk.solve(tend) |
| 72 | + # draw the line |
| 73 | + for tr in rk.Trajectory: |
| 74 | + x, y = trajPlot(tr) |
| 75 | + print(x, y) |
| 76 | + print() |
| 77 | + # draw the arrow |
| 78 | + continue |
| 79 | + l = (len(rk.Trajectory) - 1) / 3 |
| 80 | + if l > 0 and 2 * l < len(rk.Trajectory): |
| 81 | + p1 = rk.Trajectory[l] |
| 82 | + p2 = rk.Trajectory[2 * l] |
| 83 | + x1, y1 = trajPlot(p1) |
| 84 | + x2, y2 = trajPlot(p2) |
| 85 | + dx = -0.5 * (y2 - y1) # orthogonal to line |
| 86 | + dy = 0.5 * (x2 - x1) # orthogonal to line |
| 87 | + #l = math.sqrt(dx*dx + dy*dy) |
| 88 | + #if abs(l) > 1e-3: |
| 89 | + # l = 0.1 / l |
| 90 | + # dx *= l |
| 91 | + # dy *= l |
| 92 | + print(x1 + dx, y1 + dy) |
| 93 | + print(x2, y2) |
| 94 | + print(x1 - dx, y1 - dy) |
| 95 | + print() |
| 96 | + |
| 97 | +def singleTraj(system, trajStart, h=0.02, tend=1.0): |
| 98 | + tstart = 0.0 |
| 99 | + |
| 100 | + # compute the trajectory |
| 101 | + |
| 102 | + rk = RungeKutta(system, trajStart, tstart, h) |
| 103 | + rk.solve(tend) |
| 104 | + |
| 105 | + # print out trajectory |
| 106 | + |
| 107 | + for i in range(len(rk.Trajectory)): |
| 108 | + tr = rk.Trajectory[i] |
| 109 | + print(' '.join(["{:.5f}".format(t) for t in tr])) |
| 110 | + |
| 111 | +#phaseDiagram(sysSM, (lambda i, j: [0.354, 0.654, 1.278, 0.8 + 0.2 * i, 0.1 + 0.1 * j]), (lambda a: (a[4], a[5])), h=0.1, tend=math.log(10**17)) |
| 112 | + |
| 113 | +# initial conditions at M_Z |
| 114 | +singleTraj(sysSM, [0.354, 0.654, 1.278, 0.983, 0.131], h=0.1, tend=math.log(10**17)) # true values |
0 commit comments