In this case, the only branch cuts and branch points will come from the square root. The cuts of $\sqrt{f(z)}$ occurs along the half line $\text{Im}(f(z)) = 0 \,\wedge\, \text{Re}(f(z)) \leq 0$. The branch points lie at $f(z) = 0$ or $f(z) = \tilde\infty$.
Your example:
With[{z = x + I y},
expr = (Tanh[z] - Tanh[2 z])^2 + (Tanh[z] Tanh[2 z] + 1)^2 - 1 - 2 ((Tanh[2 z])^2) ((Tanh[z])^2);
branchCutRegion[x_, y_, __] = Re[expr] <= 0;
];
bpvals = Union[{x, y} /. Solve[(expr == 0 || 1/Together[TrigToExp[expr]] == 0) && -10 < x < 10 && -10 < y < 10, {x, y}]];
Here we needed to help Solve find the branch points corresponding to $\tilde\infty$.
We can visualize the cut by plotting the constraint on the imaginary part, restricted to the region defined by the constraint on the real part. Here I've overlaid the branch points:
ContourPlot[Im[expr] == 0, {x, -10, 10}, {y, -10, 10},
RegionFunction -> branchCutRegion, PlotPoints -> 100,
Epilog -> {Red, Point[bpvals]}
]

For fun we can add a plot of the expression under the cuts. Here I'll use domain coloring. Here, the complex argument varies with hue and the absolute value varies with saturation and brightness -- the darker the pixel, the larger the absolute value. I've also binned the absolute value to show some contours.
binnedabs = Compile[{{z, _Complex}},
Module[{f, abs, rnd, sgn, val},
f = (Tanh[z] - Tanh[2 z])^2 + (Tanh[z] Tanh[2 z] + 1)^2 - 1 - 2 Tanh[2 z]^2 Tanh[z]^2;
abs = Abs[f];
rnd = Round[abs, .2];
val = If[rnd == 0, f, rnd Sign[f]];
{
Divide[Mod[Arg[val], 2π], 2π],
Power[1 + 0.3*Log[Abs[val] + 1], -1],
Power[1 + 0.5*Log[Abs[val] + 1], -1]
}
],
CompilationTarget -> "C",
Parallelization -> True,
RuntimeAttributes -> {Listable},
RuntimeOptions -> "Speed"
];
lattice = Array[List, {2048, 2048}, {{-10., 10.}, {-10., 10.}}].{I, 1};
raster = Raster[binnedabs[lattice], {{-10, -10}, {10, 10}}, ColorFunction -> Hue];
cutplot = ContourPlot[Im[expr] == 0, {x, -10, 10}, {y, -10, 10},
RegionFunction -> branchCutRegion, PlotPoints -> 100, ContourStyle -> Black];
Show[
cutplot,
ImageSize -> 800,
Prolog -> raster,
Epilog -> {EdgeForm[Black], GrayLevel[.8], Disk[#, Scaled[.0045]] & /@ bpvals}
]

As of version 12 we can use ComplexPlot to visualize the domain coloring:
exprz = (Tanh[z] - Tanh[2 z])^2 + (Tanh[z] Tanh[2 z] + 1)^2 - 1 - 2 ((Tanh[2 z])^2) ((Tanh[z])^2);
exprxy = exprz /. z -> x + I y;
branchCutRegion[x_, y_, __] = Re[exprxy] <= 0;
bpvals = Union[{x, y} /. Solve[(expr == 0 || 1/Together[TrigToExp[expr]] == 0) && -10 < x < 10 && -10 < y < 10, {x, y}]];
domaincoloring = ComplexPlot[exprz, {z, -10 - 10 I, 10 + 10 I},
ColorFunction -> "CyclicLogAbsArg", ImageSize -> 800];
cutplot = ContourPlot[Im[exprxy] == 0, {x, -10, 10}, {y, -10, 10},
RegionFunction -> branchCutRegion, PlotPoints -> 100, ContourStyle -> Black];
Show[
domaincoloring,
cutplot,
Epilog -> {EdgeForm[Black], GrayLevel[.8], Disk[#, Scaled[.0045]] & /@ bpvals}
]

To achieve the same image from my original answer, you can use
domaincoloring = ComplexPlot[exprz, {z, -10 - 10 I, 10 + 10 I},
ColorFunction -> {Hue[Divide[Mod[#8, 2π], 2π],
Power[1 + 0.3*Log[#7 + 1], -1],
Power[1 + 0.5*Log[#7 + 1], -1]] &, None},
ColorFunctionScaling -> False,
Exclusions -> None,
ImageSize -> 800
];
x + I*2 yin your formula. Do you mean for that to be2x + I*2 y? $\endgroup$